Odp: nls, how to determine function?
Hi
Hi R help, I am trying to determine how nls() generates a function based on the self-starting SSlogis and what the formula for the function would be. I've scoured the help site, and other literature to try and figure this out but I still am unsure if I am correct in what I am coming up with.
Thanks for providing data and your code
**************************************************************************
dat <- c(75.44855206,NA,NA,NA,82.70745342,82.5335019,88.56617647,80. 00128866,94.15418227,86.63987539,93.91052952,74.10612245,86.62289562,90. 47961047,NA,NA,82.45320197,72.14371257,NA,71.44104803,72.59742896,68. 36363636,NA,NA,61,NA,NA,71.26502909,NA,85.93333333,84.34248284,79. 00522193,79.64223058,97.2074017,88.43700548,96.40413877,95.13511869,92. 57379057,93.97498475,NA,97.55995131,89.53321146,97.21728545,93.21980198, 77.54054054,95.85392575,86.25684723,97.55325624,80.03950617,NA,91. 34023128,92.42906574,88.59433962,65.77272727,89.63772455,NA,NA,NA,NA,74. 86344239,83.57594937,70.22516556,65.30543319,NA,NA,67.84852294,60. 90909091,54.79303797,NA,52.18735363,33.47003155,NA,41.34693878,24. 5047043,NA,NA,NA,NA,9.944444444,13.6875,NA,11.90267176,84.14285714,3. 781456954,NA,1.432926829,4.26557377,1.823529412,0.444620253,4.
711155378,NA,6.320284698,0.581632653,0.144578313,3.666666667,0,0,0,0,0,NA,
0.032947462,0,0,10.54545455,0,NA,0.561007958,0.75,NA,0.048780488,0. 74137931,NA,2.023339318,0,0,0,NA,NA,0.156950673,NA,0.283769634,32.
81818182,NA,NA,0,NA,0,0,0,NA,0.212454212,3.120181406,NA,0.011811024,NA,0,
0.120430108,5.928571429,1.75,0.679292929,0.97,NA,0,NA,NA,1,0.38547486,NA,
1.460732984,0.007795889,0.05465288,0.004341534)
dat.df.1 <- data.frame(dat)
unnecessary
dat.df.2 <- data.frame(x=x.seq, dat.df=dat.df.1)
some correction dat.df.2 <- data.frame(x=seq_along(dat), dat=dat)
fit.dat <-nls(dat ~ SSlogis(x, Asym, xmid,scal), data = dat.df.2, start =list(Asym=90, xmid = 75, scal = -6)) plot(dat.df.2, axes=FALSE, ann=FALSE, ylim=c(0,100)) lines(dat.df.2$x[complete.cases(dat.df.2)], predict(fit.dat),
ylim=c(0,100))
summary(fit.dat)
**************************************************************************
Formula: dat ~ SSlogis(x, Asym, xmid, scal)
Parameters:
Estimate Std. Error t value Pr(>|t|)
Asym 85.651 1.716 49.900 < 2e-16 ***
xmid 72.214 1.036 69.697 < 2e-16 ***
scal -6.150 0.850 -7.236 7.9e-11 ***
---
Signif. codes: 0 ?***? 0.001 ?**? 0.01 ?*? 0.05 ?.? 0.1 ? ? 1
Residual standard error: 10.33 on 105 degrees of freedom
Number of iterations to convergence: 10
Achieved convergence tolerance: 4.405e-06
(45 observations deleted due to missingness)
**************************************************************************
From r-help, SSlogis parameters asym, xmid and scal are defined as:
Asym: a numeric parameter representing the asymptote. xmid: a numeric parameter representing the x value at the inflection point of the curve. The value of SSlogis will be Asym/2 at xmid. scal: a numeric scale parameter on the input axis. and it states that the value of SSlogis "is a numeric vector of the same length as input. It is the value of the expression sym/(1+exp((xmid-input)/scal)). If all of the arguments Asym, xmid, and scal are names of objects the gradient matrix with respect to these names is attached as an attribute named gradient." However, how do I get the actual function for the curve that is generated? I don't think it can just be: y= asym/((1+e^((xmid-x)/scal)))?
Yes. I think that the best source of information about nonlinear regression is book by Bates, Pinheiro - Mixed effect models with S and S+. There you can find how to determine starting parameters, how to construct and use your own function together with selfstart feature.
Also, how do you determine the starting parameters to input in for asym, xmin, and scal? Perhaps I need to start at the beginning and define my own function, and not rely on SSlogis to provide it? What I want to be able to do is determine a local maximum for my curve (the x value at which this curve inflects (the upper inflection)), and the x value for the local minimum (the lower inflection curve), and the x value counts in between these values. I think in order to do this I need to differentiate the function.
Maybe I do not understand well but looking at the picture it seems to me that logistic model is fitting your data quite well. You can use also four parameter logistic model.
fit.dat.2 <-nls(dat ~ SSfpl(x, A, B, xmid,scal), data = dat.df.2, start
=list(A=85.65, B=0, xmid = 72, scal = -6))
summary(fit.dat.2)
Formula: dat ~ SSfpl(x, A, B, xmid, scal)
Parameters:
Estimate Std. Error t value Pr(>|t|)
A 1.6729 1.5927 1.050 0.296
B 85.5555 1.7065 50.134 < 2e-16 ***
xmid 71.7628 1.0762 66.679 < 2e-16 ***
scal -5.8051 0.9162 -6.336 6.13e-09 ***
---
Signif. codes: 0 ?***? 0.001 ?**? 0.01 ?*? 0.05 ?.? 0.1 ? ? 1
Residual standard error: 10.32 on 104 degrees of freedom
Number of iterations to convergence: 9
Achieved convergence tolerance: 7.629e-06
(45 observations deleted due to missingness)
As you can see parameter A is insignificant so simple logistic can be used
too. In that case upper asymptote is 85.6, lower asymptote is zero,
inflection point is 72 (x where y is in the middle between both
asymptotes) and scal is rate at which the curve is falling (growing).
There is however some wave in the beginning of your data
fit <-loess(dat ~ x, data = dat.df.2, span=0.3)
lines(dat.df.2$x[complete.cases(dat.df.2)], predict(fit), col=3)
So it is up to you to decide if you are satisfied with getting asymptotic
values from logistic model or you want to set something more elaborated.
Regards
Petr
Any insight on this would be greatly appreciated. Sincerely, Katrina
______________________________________________ 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.