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

[Ben Bolker]

On 11-09-20 08:05 PM, Gabriel Yvon-Durocher wrote:
> Ben and Doug,
> 
> Thank you very much for your suggestions. 
> 
> I managed to get nlme to work really nicely (see below code). The model
> diagnostics, and examination of the fitted values suggest it fits the
> data well.
> 
> Ben, I tried your code for the lmer but it gave very odd parameter
> estimates. These are key, because I use these in the derivative function
> of the polynomial to estimate Q10 at any temperature.
> 

  Glad it worked.

  My orthogonal-polynomial stuff has the advantages of being a bit more
numerically stable, and of allowing independent tests of the linear,
quadratic, cubic components of the model, BUT at some cost to
interpretability.  The parameters will not be substitutable into your
derivs function unless your derivs function is itself defined in terms
of the orthogonal parameterization.

  If you like the results you got, then it doesn't seem worth the
trouble redoing it that way.

   Ben


> Thanks again for your help.
> 
> successful code for reference:
> 
> # groupedData object
> structure<-groupedData(log10.rate~round.temp|Subject, outer=~Stream,
> data=lab.analysis)
> 
> # respiration function - 3rd order polynomial
> resp.fun<-function(Temp,a,b,c,d) a*Temp^3-b*Temp^2+c*Temp-d
> 
> # differentiate polynomial function
> gr.model<- deriv(body(resp.fun), namevec = c('a','b','c','d'),
> function.arg = resp.fun)
> 
> control<-nlmeControl(maxIter=500, msMaxIter=500, tolerance=1e-6)
> 
> # build model using nlme
> model<-nlme(log10.rate ~ gr.model(Temp,a,b,c,d),
> fixed = a + b + c + d ~ 1,
> random = a + b + c + d ~ 1,
> start =c(0.0000008,0.0000655,0.0238379,1.9547086),
> data = structure,
> verbose=FALSE,
> control=control)
> 
> 
> 
> On 20 Sep 2011, at 17:15, Ben Bolker wrote:
> 
>> 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
>>> <mailto: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
>>>>>> <mailto: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
>>>> <mailto:R-sig-mixed-models at r-project.org> mailing list
>>>> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>>>>
>>
>