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

Highland Statistics Ltd highstat at highstat.com
Wed Nov 16 21:07:31 CET 2011 





------------------------------

Message: 7
Date: Wed, 16 Nov 2011 14:20:18 -0500
From: Maria Susana Alvarado Barrientos<masaba79 at gmail.com>
To: "r-sig-mixed-models at r-project.org"
	<r-sig-mixed-models at r-project.org>
Subject: [R-sig-ME] Variance structure and random effects using lme,
	what's the difference?
Message-ID:
	<CAB62b8bhh9hYD__ifZDSwmktA-X8XWGWFJ6DWwt_rjTyhE=GZg at mail.gmail.com>
Content-Type: text/plain; charset="windows-1252"

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





Maria,

This is what I would do.
1. Fit an lme with random intercept fTree....just by design of the study. Some people will not even bother to test whether the random intercept is needed...which is a valid point.
But you could do this.
2. Check the lme for heterogeneity. And if there is heterogeneity trouble...choose something way more simplistic than your varIdent by tree. You were estimating a lot (!!) of variances
(hence the error). It would also have been better/easier if the heterogeneity issue would have been part of your underlying questions.

Good luck,

Alain











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.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-0001.png>

------------------------------

_______________________________________________
R-sig-mixed-models mailing list
R-sig-mixed-models at r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models


End of R-sig-mixed-models Digest, Vol 59, Issue 22

More information about the R-sig-mixed-models mailing list