[R] error using nls with logistic derivative

[Francisco Mora Ardila]

It seems to work!

Thanks and sorry for the personal message

Francisco

On Tue, 17 Apr 2012 09:09:24 +0200, peter dalgaard wrote
> On Apr 17, 2012, at 06:23 , Francisco Mora Ardila wrote:
> 
> > Hi
> > 
> > I´m trying to fit a nonlinear model to a derivative of the logistic function
> > 
> > y = a/(1+exp((b-x)/c)) (this is the parametrization for the SSlogis function with 
nls)
> > 
> > The derivative calculated with D function is:
> > 
> >> logis<- expression(a/(1+exp((b-x)/c)))
> >> D(logis, "x")
> > a * (exp((b - x)/c) * (1/c))/(1 + exp((b - x)/c))^2
> > 
> > So I enter this expression in the nls function:
> > 
> > ratelogis <- nls(Y ~ a*(exp((b-X)/c)*(1/c))/(1 + exp((b-X)/c))^2,
> > start=list(a = 21.16322, b = 8.83669, c = 2.957765),
> > )
> > 
> > The data is:
> > 
> >> Y
> > [1]  5.5199668  1.5234525  3.3557000  6.7211704  7.4237955  1.9703127
> > [7]  4.3939336 -1.4380091  3.2650180  3.5760906  0.2947972  1.0569417
> >> X
> > [1]   1   0   0   4   3   5  12  10  12 100 100 100
> > 
> > The problem is that I got the next error:
> > 
> > Error en nls(Y ~ a * (exp((b - X)/c) * (1/c))/(1 + exp((b - X)/c))^2,  : 
> >  step factor 0.000488281 reduced below 'minFactor' of 0.000976563
> > 
> > I trien to change the minFactor using the control argument inside nls
> > 
> > control=nls.control(maxiter=50, tol=1e-5, minFactor = 1/2048
> > 
> > but got a new error message:
> > 
> > 
> > Error en nls(Y ~ a * (exp((b - X)/c) * (1/c))/(1 + exp((b - X)/c))^2,  : 
> >  step factor 0.000244141 reduced below 'minFactor' of 0.000488281
> > 
> > So it seems that as I modify minFactor, the step factor reduces also and I can never 
> > reach a solution.
> > 
> > Does anybody Know what am I doing wrong? Is there a problem with the formula? How 
can I 
> > solve it? I tried some suggestions on R-help related topics but did not work.
> 
> (Please don't send private messages. I don't do free consulting outside of the 
> mailing lists.)
> 
> With nls(), there's really no substitute for good data and good starting 
> values. Did you actually try plotting those data? At best, they are extremely 
> noisy. How did you obtain the starting values? They seem remarkably accurate 
> for something that fits data so poorly.
> 
> What is happening is that the algorithm shoots off into a region of the 
> parameter space where it can't make any progress:
> 
> > ratelogis <- nls(Y ~ a*(exp((b-X)/c)*(1/c))/(1 + exp((b-X)/c))^2,+ start=list(a = 
21.16322, b = 8.83669, c = 2.957765),+ trace=TRUE)
> 151.098 :  21.163220  8.836690  2.957765 
> 127.1149 :  103.49326  11.43274  20.29663 
> 111.344 :  833.02386 -45.86244 140.32985 
> 111.3375 :  978.97214 -76.20571 159.90833 
> 111.3374 :  1097.1376 -101.6771  174.2037 
> 111.3228 :  1179.8451 -119.7416  183.3794
> 
> Notice that the "b" parameter which is supposed to be the maximum point starts 
> off well to the right of the position of the largest Y values, then shoots 
> into large negative values.
> 
> With a little eyeballing, I can do better:
> 
> > ratelogis <- nls(Y ~ a*(exp((b-X)/c)*(1/c))/(1 + exp((b-X)/c))^2,
> + start=list(a = 40, b = 3.5, c = 10),trace=T)
> 139.5173 :  40.0  3.5 10.0 
> 97.28614 :  26.513424  4.052639  2.267709 
> 81.53127 :  32.384910  3.473411  2.307504 
> 65.11387 :  39.53542  3.01097  2.07678 
> 50.66328 :  36.223529  2.590102  1.133965 
> 50.50921 :  35.984466  2.565698  1.067731 
> 50.50162 :  36.149993  2.573420  1.081673 
> 50.50145 :  36.129962  2.572195  1.079504 
> 50.50144 :  36.133656  2.572402  1.079862 
> 50.50144 :  36.133061  2.572368  1.079803
> 
> However, the fit isn't exactly stellar. Try this:
> 
> X <- c(1, 0, 0, 4, 3, 5, 12, 10, 12, 100, 100, 100)
> Y <- c(5.5199668, 1.5234525, 3.3557, 6.7211704, 7.4237955, 1.9703127, 
>    4.3939336, -1.4380091, 3.265018, 3.5760906, 0.2947972, 1.0569417)
> plot(X,Y)
> ratelogis <- nls(Y ~ a*(exp((b-X)/c)*(1/c))/(1 + exp((b-X)/c))^2,
>  start=list(a = 40, b = 3.5, c = 10),trace=T)
> x0 <-  seq(0,100,,1000)
> lines(x0,predict(ratelogis,newdata=data.frame(X=x0)))
> lines(x0,evalq(a*(exp((b-X)/c)*(1/c))/(1 + exp((b-X)/c))^2,envir=list(a = 
> 21.16322, b = 8.83669, c = 2.957765, X=x0)), lty="dashed")
> 
> -- 
> 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


----------------------
Francisco Mora Ardila
Estudiante de Doctorado
Centro de Investigaciones en Ecosistemas
Universidad Nacional Autónoma de México