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