[R-sig-ME] algal nonlinear mixed effects problem

[Reinhold Kliegl]

Thanks, Ben, for the illustration of how to use derive(), deparse(),
and eval() for setting up such functions. The following is based on
lme4_0.999999-0.

I modified Ben's attempt at a fit in two ways to get convergence.
(1) The start argument threw an Error; I provide the values in a
vector rather than a list. This might be related to differences
between nlmer implementations.
(2) The start value for asymp.R2 had to be moved away from 0.

The result suggests that there is no reliable evidence for the
variance component associated with "Individual asymp.R1"

A part of the protocol is listed below.
Happy New Year to everyone.

Reinhold Kliegl

Here is the protocol:
> # 4. Attempt the fit
> nlmerfit2 <- nlmer(
+   X ~ fpl(Day, asymp.L, asymp.R1, asymp.R2, asymp.R3, xmid, scale) ~
+          asymp.R1|Individual,
+      start =  list(nlpars=c(asymp.L=0.7,
+      asymp.R1=0.6,asymp.R2=0,asymp.R3=0,xmid=5,scale=1)),data=d)
Fehler: length(start$fixef) > 0 is not TRUE
> # ERROR: length(start$fixef) > 0 is not TRUE
>
> nlmerfit2a <- nlmer(
+   X ~ fpl(Day, asymp.L, asymp.R1, asymp.R2, asymp.R3, xmid, scale) ~
+          asymp.R1|Individual,
+      start =  c(asymp.L=0.6,
asymp.R1=0,asymp.R2=0,asymp.R3=0,xmid=5,scale=1),data=d)
Fehler in mer_finalize(ans) : Downdated X'X is not positive definite, 5.
> # Fehler in mer_finalize(ans) : Downdated X'X is not positive definite, 5.
>
> nlmerfit2b <- nlmer(
+   X ~ fpl(Day, asymp.L, asymp.R1, asymp.R2, asymp.R3, xmid, scale) ~
+          asymp.R1|Individual,
+      start =  c(asymp.L=0.6,
asymp.R1=.1,asymp.R2=0,asymp.R3=0,xmid=5,scale=1),data=d)
> nlmerfit2b
Nonlinear mixed model fit by the Laplace approximation
Formula: X ~ fpl(Day, asymp.L, asymp.R1, asymp.R2, asymp.R3, xmid,
scale) ~      asymp.R1 | Individual
   Data: d
   AIC   BIC   logLik deviance
 16.11 32.02 -0.05555   0.1111
Random effects:
 Groups     Name     Variance   Std.Dev.
 Individual asymp.R1 8.9101e-22 2.985e-11
 Residual            2.0575e-03 4.536e-02
Number of obs: 54, groups: Individual, 9

Fixed effects:
           Estimate Std. Error t value
asymp.L  -2.173e+02  2.187e+05  -0.001
asymp.R1  6.140e+01  5.002e+04   0.001
asymp.R2 -1.917e-02  4.370e+00  -0.004
asymp.R3 -1.207e-01  2.752e+01  -0.004
xmid      6.152e+03  6.167e+06   0.001
scale     4.807e+03  3.610e+06   0.001

Correlation of Fixed Effects:
         asym.L asy.R1 asy.R2 asy.R3 xmid
asymp.R1 -0.357
asymp.R2 -0.677 -0.446
asymp.R3 -0.677 -0.446  1.000
xmid     -1.000  0.357  0.677  0.677
scale    -0.597  0.962 -0.185 -0.185  0.597


On Mon, Dec 31, 2012 at 10:09 PM, Ben Bolker <bbolker at gmail.com> wrote:
> Ben Bolker <bbolker at ...> writes:
>
>>
>>
>>   Inspired by a recent question of Diego Pujoni's on the list
>> ( http://article.gmane.org/gmane.comp.lang.r.lme4.devel/9363/match= ),
>> I started working on the problem, with two realizations:
>>
>>  1. the problem really needs a nonlinear fit: see
>> http://rpubs.com/bbolker/3363
>>
>>  2.  I have no idea how to do an nlmer fit with a non-trivial
>> fixed-effects model (e.g. in nlme we could do
>>
>> nlmefit1 <- nlme(model = X ~ SSfpl(Day, asymp.L, asymp.R,
>>                                    xmid, scale),
>>    fixed = list(asymp.R ~ Group, xmid + scale + asymp.L ~ 1),
>>    random = asymp.R ~ 1 | Individual, ...)
>>
>>   This was raised earlier on stack overflow:
>>
>>
> http://stackoverflow.com/questions/11056625/how-to-add-fixed-effect-to-four-parameter-logistic-model-in-nlmer
>>
>
>   Reinhold Kliegl e-mailed off-list to say that he was able to do
> it / had to do it by hacking up an objective function that accounts
> for the grouping (i.e., more or less doing it by hand).  I have
> implemented this and updated http://rpubs.com/bbolker/3363 , which
> also shows a solution with AD Model Builder; it turns out that at
> least the development version of nlmer is a little too fragile to
> fit the model I wanted to fit (I haven't tested with the stable
> version).
>
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models