[R] SSbiexp

[peter dalgaard]

On Apr 23, 2013, at 14:05 , catalin roibu wrote:

> Hello all!
> I have a problem to use a biexponential regression model:
> I use this code:
> n3<-nls(proc~SSbiexp(cls,a,b,c,d),data=bline) and this is the error message:
> 
> Error in nls(y ~ cbind(exp(-exp(lrc1) * x), exp(-exp(lrc2) * x)), data =
> xy,  :
>  singular gradient
> 

Well, the basic issue is that the model doesn't fit well.

First, try 

plot(proc~cls, bline, log="y")

and notice that this is rather close to a straight line, i.e. a monoexponential curve fits the original data rather well, which indicates that it can be difficult to find sufficient information to fit a biexponential.

Next, notice that the model is partially linear, so we can try plotting the residual SS as function of the log rate-constants, i.e.

attach(bline) # (a bit sloppy, but gets the job done...)
f <- function(a,b) {
 x1 <- exp(-exp(a)*cls)
 x2 <- exp(-exp(b)*cls)
 r <- lm(proc ~ x1+x2-1)$residual
 sum(r^2) 
}
grd <- seq(-4,4,,101)
grdval <- outer(grd,grd,Vectorize(f))
image(log(grdval))

There's some spiking along the diagonal, but that's not of interest here. The interesting bit is that there's no definite minimum SS. Apparently, you can chose one of the rate constants as large as you please and fix the other one at around 0.4. Some further digging shows that the coefficient to x1 for the minimum cell above is -2.64e+08   while x1 itself is 

> exp(-exp(4)*cls) 
[1] 1.180541e-06 1.393678e-12 1.645294e-18 1.942338e-24 2.706993e-36
[6] 3.772675e-48 5.257894e-60 7.327809e-72 1.021260e-83 1.423308e-95

Notice that this multiplied by -2.64e08 is essentially zero except for the first element. I.e., the effect of the second exponential is essentially to take the first observation out of the data. However, that amounts removing one observation using two parameters, which is probably what causes the singularity.

-pd




> My data is like this:
> 
> structure(list(proc = c(387.177462830167, 508.090511433077,
> 321.836817656365,
> 151.226805860727, 150.885315536942, 86.2895998400175, 56.3215592398958,
> 39.5044440630984, 25.5703078997907, 7.33494872198676), cls = c(0.25,
> 0.5, 0.75, 1, 1.5, 2, 2.5, 3, 3.5, 4)), .Names = c("proc", "cls"
> ), row.names = c("0.25", "0.5", "0.75", "1", "1.5", "2", "2.5",
> "3", "3.5", "4"), class = "data.frame")
> 
> Thank you!
> 
> -- 
> ---
> Catalin-Constantin ROIBU
> Lecturer PhD, Forestry engineer
> Forestry Faculty of Suceava
> Str. Universitatii no. 13, Suceava, 720229, Romania
> office phone     +4 0230 52 29 78, ext. 531
> mobile phone   +4 0745 53 18 01
>                       +4 0766 71 76 58
> FAX:                +4 0230 52 16 64
> silvic.usv.ro
> 
> 	[[alternative HTML version deleted]]
> 
> ______________________________________________
> R-help at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.

-- 
Peter Dalgaard, Professor
Center for Statistics, Copenhagen Business School
Solbjerg Plads 3, 2000 Frederiksberg, Denmark
Phone: (+45)38153501
Email: pd.mes at cbs.dk  Priv: PDalgd at gmail.com