nls, how to determine function?
Hi Petr, thanks for your help on this. I will most definitely get this book as it appears to be a good one. I am happy with how nls() and the self starting function appears to be fitting the data set. What I am wondering about is how to write out this function outside of R. For example, if I sat down with a paper and pencil and wrote this function out in terms of f(x) = ?. In this example (http://www.graphpad.com/curvefit/id203.htm) there is a formula contains for what they are calling the Boltzmann Sigmoid function. In order to have an answer, is this something I need to generate on my own? Or, is there a common form that is being applied by the SSlogis and the SSfpl functions? In Pinheiro and Bates (pg 518) there is a model f(x) = A + B/1+ exp[(scal-x)/expL]. Is this the form of the equation? Or, is there a way to extract the local min and max values from nls () results? Do I have to break apart the individual components of the regression into piecewise functions? Thanks again for your help. Katrina
On Tue, Aug 9, 2011 at 12:38 AM, Petr PIKAL <petr.pikal at precheza.cz> wrote:
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.