[R-sig-ME] non-linear mixed effects: polynomial

[Ben Bolker]

On 09/20/2011 10:46 AM, Douglas Bates wrote:
> On Tue, Sep 20, 2011 at 8:17 AM, Ben Bolker <bbolker at gmail.com> wrote:
>> On 09/20/2011 04:22 AM, Gabriel Yvon-Durocher wrote:
>>> Thanks for the reply Ben. I have 16 subjects, however the number of
>>> observations per subject is not even, which is why I opted for nlmer
>>> rather than nlme. I presumed nlmer (as is the case for lmer) would be
>>> better at handling unbalanced designs.
>>
>>  No, nlme is also intended for unbalanced designs.
>>
>>  nlmer is likely to be (a) faster and (b) better for crossed designs,
>> but otherwise I would stick (for the time being) to nlme.
> 
> I don't think you need either nlme or nlmer though, do you?  A
> polynomial model is linear in the coefficients so this could be fit by
> lmer. 


D'oh.

 In an ideal world I would recommend

 rate~poly(Temp,3)

for the fixed effect specification (unless you really need to be able to
interpret the slope, quadratic, cubic factors at Temp==0 (!!)), but I
think this could be hard to manage with the random effect specification.
 I might do something like:

polymat <- model.matrix(~poly(Temp,3),data=mydat)
colnames(polymat) <- c("int","linT","quadT","cubicT")

newdat <- data.frame(log10.rate=mydat$log10.rate,polymat[,-1])

lmer(log10.rate~linT+quadT+cubicT+(1|Subject)+(0+linT|Subject)+ ....


 Ben's caution about over-specification of the random-effects
> variance-covariance model still applies to the linear mixed-effects
> model, however.  You are trying to estimate 15 variance-covariance
> parameters, for which you would need a great deal of data.
> 
>>   A couple of other things occur to me:
>>
>> (1) what you got was a warning about false convergence (the PORT library
>> algorithm built into lme4 is pretty finicky ...) If the result seems
>> sensible and accords with your individual-level fits, I would probably
>> proceed (with caution)
>>
>> (2) the random-effects specification (1+a+b+c+d|Subject) is likely to be
>> quite tricky (I guess you have lots of data per individual?) because it
>> will fit a full, unstructured 5 x 5 variance-covariance matrix for the
>> random effects across individuals.  It's a pain, but maybe try
>>
>> (1|Subject)+(0+a|Subject)+(0+b|Subject)+(0+c|Subject)+(0+d|Subject)
>>
>>  (i.e. independent random effects) for comparison?  (Admittedly the
>> among-subject correlations among parameters could be quite interesting,
>> although probably more interpretable if you centered your Temperature
>> variable (another thing worth trying to improve convergence) or used
>> orthogonal polynomials (ditto).
>>
>>
>>>
>>> I will check out the suggested books. Thanks.
>>>
>>> Gabriel Yvon-Durocher
>>>
>>>> Gabriel Yvon-Durocher <g.yvon-durocher at ...> writes:
>>>>
>>>>>
>>>>> Dear all,
>>>>>
>>>>> I am trying to fit a 3rd order polynomial to some rate and
>>>>> temperature
>>>> measurements, but am running into
>>>>> problems. I would like to treat the parameter a,b,c,d as random
>>>>> effects across
>>>> subjects, and use the model
>>>>> to estimate these parms for each experimental subject.
>>>>>
>>>>> Here is my model:
>>>>>
>>>>> resp.fun<-function(Temp,a,b,c,d)
>>>>> -a*(Temp^3)+b*(Temp^2)-c*(Temp)+d
>>>>>
>>>>> then I build a gradient function:
>>>>>
>>>>> gr.model<- deriv(body(resp.fun), namevec = c('a','b','c','d'),
>>>> function.arg = resp.fun)
>>>>>
>>>>> and then run the mixed effects model:
>>>>>
>>>>> model<-nlmer(log10.rate ~ gr.model(Temp,a,b,c,d) ~
>>>>> 1+a+b+c+d|Subject, start = c(a=4.897e-05, b=3.198e-03,
>>>>> c=3.569e-02, d=1), data = lab.analysis, na.action=na.omit)
>>>>>
>>>>> unfortunately I get the error message
>>>>>
>>>>> Warning message: In mer_finalize(ans) : false convergence (8)
>>>>>
>>>>> Any ideas on what I may be doing wrong? I have manually fit this
>>>>> polynomial to the data for each subject and the fits are good,
>>>>> but I want to run the mixed effects analysis to get the best
>>>>> possible parameter estimates given the experimental design.
>>>>>
>>>>
>>>> How many subjects?
>>>>
>>>> Can you do it in nlme instead, which is more mature/robust? If not,
>>>> you may end up needing to do this in AD Model Builder, or BUGS, or
>>>> rolling your own Laplace approximator (there is an example in
>>>> Madsen and Thyregod's book) ...
>>>>
>>>> Ben Bolker
>>>>
>>>> _______________________________________________
>>>> R-sig-mixed-models at r-project.org mailing list
>>>> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>>>
>>
>> _______________________________________________
>> R-sig-mixed-models at r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>>