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]]
>
> ______________________________________________
> [email protected] 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: [email protected] Priv: [email protected]
______________________________________________
[email protected] 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.
