[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