[R] nls, how to determine function?

Tue Aug 9 21:02:53 CEST 2011

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.