R help - May 2012 - nls and if statements

Hi All,

I have a situation where I want an 'if' variable to be parameterized.
It's
entirely possible that the way I'm trying to do this is wrong, especially
given the error message I get that indicates I can't do this using an
'if'
statement. 

Essentially, I have data where I think a relationship enters when a variable
(here Pwd) is below some value (z). I don't know that value, so I want to
fit it. The data is Pw, Tsoil, and Cfl. 

nlstest=nls(if(Pw>z){Cfl~K*exp(-b*Pw)+c
	}else{
	Cfl~K*(exp(-a*Tsoil))*exp(-b*Pw)+c
	},start=c(K=5.5, a=0.1, b=0.1, c=2, z=5))

Which returns the error:
Error in inherits(object, "formula") : object 'z' not found

Is there a better way to try and allow a conditional variable to, well,
vary?

--
View this message in context:
http://r.789695.n4.nabble.com/nls-and-if-statements-tp4630391.html
Sent from the R help mailing list archive at Nabble.com.

Hello,

I don't know if this is what you want, but your formula is equivalent to

nlstest=nls(Cfl~K*ifelse(Pw > z, 1, exp(-a*Tsoil))*exp(-b*Pw)+c
        ,start=c(K=5.5, a=0.1, b=0.1, c=2, z=5)) 

Hope this helps,

Rui Barradas


--
View this message in context:
http://r.789695.n4.nabble.com/nls-and-if-statements-tp4630391p4630398.html
Sent from the R help mailing list archive at Nabble.com.

On May 17, 2012, at 10:06 AM, DWatts wrote:
> Hi All,
>
> I have a situation where I want an 'if' variable to be  
> parameterized. It's
> entirely possible that the way I'm trying to do this is wrong,  
> especially
> given the error message I get that indicates I can't do this using  
> an 'if'
> statement.
>
> Essentially, I have data where I think a relationship enters when a  
> variable
> (here Pwd) is below some value (z). I don't know that value, so I  
> want to
> fit it. The data is Pw, Tsoil, and Cfl.
>
> nlstest=nls(if(Pw>z){Cfl~K*exp(-b*Pw)+c
> 	}else{
> 	Cfl~K*(exp(-a*Tsoil))*exp(-b*Pw)+c
> 	},start=c(K=5.5, a=0.1, b=0.1, c=2, z=5))
>
Perhaps recasting as Boolean equivalent:

nlstest== <- nls({Cfl~ (Pw>z)*K*exp(-b*Pw)+c) +
                        !(Pw>z)*(K*(exp(-a*Tsoil))*exp(-b*Pw)+c)},
              start=c(K=5.5, a=0.1, b=0.1, c=2, z=5))
> and provide commented, minimal, self-contained, reproducible code.
Untested in absence of data.

-- 

David Winsemius, MD
West Hartford, CT

Your issue is that nls returns "singular gradient", but I believe the
real trouble is that
if or ifelse are not in the derivatives table. My own (experimental, but in
R-forge)
package nlmrt has nlxb to mirror nls but be much more aggressive in finding a
solution,
and it gives an error msg that ifelse is not in the derivatives table. It may be
that nls
has an alternative for getting numerical derivatives (perhaps Doug Bates or
Saikat DebRoy
could comment).

It may be better to convert to a function and use Nelder-Mead in optim, or
better, nmk
from dfoptim which is more up to date. I'm uniquely positioned to say this,
as I wrote the
version of Nelder-Mead in optim and I think it's getting long in the tooth,
though it
still works remarkably well. One of the routines in minqa may be better for
speed, but
possibly less stable.

Code:


wfn<-function(par, Cfl, Tsoil, Pw){
   K<-par[1]
   a<-par[2]
   b<-par[3]
   c<-par[4]
   z<-par[5]
   m<-length(Pw)
   res<-rep(NA,m)
   for (i in 1:m){
     tmp<-1
     if (Pw[i] <= z) tmp<-exp(-a*Tsoil[i])
     ff<-K*tmp*exp(-b*Pw[i])+c
     res[i]<-Cfl[i]-ff
   }
   as.numeric(crossprod(res))
}
a1<-optim(st,wfn, control=list(trace=1), Cfl=Cfl, Tsoil=Tsoil, Pw=Pw)

(which converges pretty fast)



Best, JN