[R-sig-ME] Variance structure and random effects using lme, what's the difference?

[Maria Susana Alvarado Barrientos]

Dear all,

I have a question about variance structure and random effects using lme.

To give a little context, I want to develop a function to predict a
parameter, Cs, related to hourly changes in the radial profile of sap
velocity, with a few meteorological variables (VPD and S; and possibly HH,
the hour of day, as there is hysteresis in the relation of Cs and VPD, and
Cs and Sin). The data is nested: Cs for each hour of many days for 18 trees
from 2 different stands. The data are unbalanced as I have different number
of Cs data points for each tree.

I thought of fitting a mixed effects model using Tree (as a factor) as
random effects (random slope and random intercept), to be able to make
generalizations.

Here is an example of the dataset

> head(CSData)

         Cs SerialTime   YY MM HH        Sin  VPD PP Tree Stand

1 15.979194   733801.6 2009  1 14 937.166670 0.97  0    1     1

2 10.828222   733801.6 2009  1 15 909.500000 0.83  0    1     1

3  8.431478   733801.7 2009  1 16 782.350000 0.65  0    1     1

4  3.902576   733801.7 2009  1 17 335.866670 0.43  0    1     1

5  4.468554   733802.3 2009  1  7  -8.227667 0.82  0    1     1

6 12.060395   733802.3 2009  1  8  -3.882500 0.92  0    1     1

> tail(CSData)

             Cs SerialTime   YY MM HH     Sin    VPD PP Tree Stand

12726 1.8884800   734285.5 2010  5 11 605.983 -0.188  0   21     2

12727 0.6314967   734285.5 2010  5 12 468.333 -0.130  0   21     2

12728 0.3366606   734286.4 2010  5  9 267.700 -0.206  0   21     2

12729 0.1030525   734286.5 2010  5 12 153.480 -0.190  0   21     2

12730 4.5476530   734287.5 2010  5 12 494.350 -0.024  0   21     2

12731 4.2261614   734289.6 2010  5 15 346.517 -0.117  0   21     2

* *

I’m learning from Zuur et al. 2009 book, so I fitted this model with gls
first:


>full=formula(Cs ~ VPD*Sin + Sin*HH + VPD*HH) #full fixed effects

>M.gls=gls(full,data=CSData,method="REML")


After checking the diagnostics plots (attaching two: Normalized residuals
vs VPD per tree, and Normalized residuals vs Tree ID) I have decided that
this model including a variance structure for Tree and VPD is promising to
deal with heterogeneity in my data:


>full = formula(Cs ~ VPD*Sin + Sin*HH + VPD*HH)

>vf2 = varComb(varIdent(form=~1|fTree), varExp(form=~VPD))

>M.gls4 = gls(full, data=CSData, weights=vf2, method="REML")

Everything seems fine until I add the random effects, because I get the
error below:


>M1.lme = lme(full, data=CSData, random=~1|fTree, weights=vf2, method="REML")

Error in lme.formula(full, data = CSData, random = ~1 | fTree, weights
= vf2,  :

  nlminb problem, convergence error code = 1

  message = iteration limit reached without convergence (9)


I don't know what this error means, but then again I am not sure if the
problem is conceptual (the model parametrization might be wrong). Then I
would like to hear your thoughts on what am I really doing by adding both a
variance structure to allow heterogeneity per tree and at the same time
adding the trees identity as random effect?


I'd appreciate very much your help with comments and suggestions,


Sincerely,


M. Susana Alvarado Barrientos

PhD candidate - NRESS program

Dept. Natural Resources & the Environment

University of New Hampshire
-------------- next part --------------
A non-text attachment was scrubbed...
Name: Rplot.png
Type: image/png
Size: 11714 bytes
Desc: not available
URL: <https://stat.ethz.ch/pipermail/r-sig-mixed-models/attachments/20111116/9a89b9b6/attachment-0006.png>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: NormResVSvpd.png
Type: image/png
Size: 12731 bytes
Desc: not available
URL: <https://stat.ethz.ch/pipermail/r-sig-mixed-models/attachments/20111116/9a89b9b6/attachment-0007.png>