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

[Reinhold Kliegl]

I agree that the estimates are not acceptable; I was just trying to
show that one can get rid of the nlmer() error messages. I also put up
data to show that, in principle, one can get interpretable results
with nlmer() for "well-behaved" data, but apparently the program is
not quite ready for prime time.

As far as the present set of data is concerned, first I would probably
remove the asymp.R1 variance component, but I am not sure in which
direction to take the model from there. I doubt that the logistic
function is a good reference for all three groups; I only see the
logistic function justified for the t2-group. I know about "borrowing
strength", but in this case, I would probably rather borrow strength
to fit the t2-group to the other two. I did not study your ADMB
results in detail, but if your happy with them, this probably is a
defensible solution.


On Tue, Jan 1, 2013 at 8:30 PM, Ben Bolker <bbolker at gmail.com> wrote:
> 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
>