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.

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@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@...> 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@r-project.org mailing list
>>>>> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>>>> 
>>> 
>>> _______________________________________________
>>> R-sig-mixed-models@r-project.org mailing list
>>> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>>> 
> 


	[[alternative HTML version deleted]]


