[R-sig-ME] compound symmetry correlation structure in nlme

[Ben Bolker]

-----BEGIN PGP SIGNED MESSAGE-----
Hash: SHA1

On 06/16/2011 03:32 PM, Rian Dickson wrote:
> hello all,
> 
> I am analyzing some data with repeated measures on individuals, and
> want to test different correlation structures. In Zuur et al (2007)
> they state "the inclusion of a random intercept in a GLM is imposing
> the compound symmetrical correlation structure, just as it did in the
> linear mixed model" (Chapter 13, p.323), which I understand as
> meaning that if a random intercept is included in a model, then the
> compound symmetry correlation structure is implied.	
> 
> When I run the following two models (which are identical, except that
> in the second one I explicitly specify the correlation structure), I
> get the same estimates for the random and fixed effects, but the AIC
> and BIC values differ, because the second model has one more
> parameter (for the correlation structure).
> 
> global.ri.lme <- lme(Log.minH ~ site + f.year + cohort + emerg +
> emergSQ + pri + priSQ + start.t + startSQ + sea + tidem + tided, data
> = SUSCforage, random = ~1|scoterID, method = "REML")
> 
> global.ri.cs1.lme <- lme(Log.minH ~ site + f.year + cohort + emerg +
> emergSQ + pri + priSQ + start.t + startSQ + sea + tidem + tided, data
> = SUSCforage, random = ~1|scoterID, correlation = corCompSymm(form =
> ~1|scoterID), method = "REML")
> 
> So, if I want to obtain an accurate AIC value, which should I use?
> 

   Off the top of my head, I would say that since AIC values are only
defined in a relative sense, you should just be consistent throughout
your analysis -- always use one or the other.  The only time when
absolute AIC values matter is in the case of AICc.  Then, arguably,
nuisance variables such as scale parameters that can be computed without
doing any estimation should not be counted.

  Pinheiro and Bates 2000 make the point (p. 227-228, paperback edition)
that corCompSymm is not quite equivalent to the intercept-only model:
"The two correlation models are not equivalent, however, because the
intraclass correlation in the linear mixed-effects model can only take
values between 0 and 1, while (5.23) allows \rho to take negative values
(to have a positive-definite compound symmetry correlation structure, it
is only required that $\rho > -1/[max_{i\le M} (n_i) -1]$).

  Gang Chen kindly pointed this out to me when I asked a question
recently as to why the results of paired t tests and REML fits of
intercept-only models with two factors per block did not agree:
<https://stat.ethz.ch/pipermail/r-sig-mixed-models/2011q2/006091.html>

  Bottom line: I'm not sure of the answer, but I would hope that you can
get away with simply being consistent.  There are a variety of other
issues surrounding AIC and random effects (see
<http://glmm.wikidot.com>, as usual ...)  I also hope that the decision
to add or subtract 1 from the number of parameters does not
qualitatively affect your conclusions -- if it does then they're pretty
fragile ...

  good luck,

  Ben Bolker
-----BEGIN PGP SIGNATURE-----
Version: GnuPG v1.4.10 (GNU/Linux)
Comment: Using GnuPG with Mozilla - http://enigmail.mozdev.org/

iEYEARECAAYFAk36gvIACgkQc5UpGjwzenMj7gCff2evcCMcb7xt8pstm0pN3jPD
aSMAoJI0i2ZJ/IjCdUby6FCnnAgf+TXf
=528J
-----END PGP SIGNATURE-----