[R] error using nls with logistic derivative

peter dalgaard pdalgd at gmail.com
Tue Apr 17 09:09:24 CEST 2012


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



More information about the R-help mailing list