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

[Ben Bolker]

On 13-01-01 05:17 AM, Reinhold Kliegl wrote:
> 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"

  The results are crazy, though -- take a look at the fixed effect
coefficients, or the predicted values ...

> 
> 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