[R-sig-ME] Specifying correlation structure in lme
Ben Bolker
bbolker at gmail.com
Tue May 8 19:22:34 CEST 2012
Angelina Mukherjee <angelina.mukherjee88 at ...> writes:
> I'm trying to use AR1 as a correlation structure in modeling data. My
> variables, Patient and Region, are factors with Region nested in Patient.
> Expression is my Response variable.
>
> *Patient <- c(rep('1', 10), rep('2', 12))*
>
> *Region <- c(rep('1', 3), rep('2', 3), rep('3', 4), rep('4', 4), rep('5',
> 4), rep('6', 4))*
> library(nlme)
Are you really using character variables? I don't think that works
(although other aspects of your coding appear to have been mangled
in transition, so maybe that's another one).
I can replicate your error as follows:
df <- data.frame(Patient=factor(rep(1:2,c(10,12))),
Region=rep(1:6,c(3,3,4,4,4,4)),
Expression=rnorm(22))
library(nlme)
m0 <- lme(Expression ~ Patient + Region,
random = ~1 |Patient/Region,
correlation=corAR1(form=~1 | Patient/Region), data=df)
lme(Expression ~ Patient + Region,
random = ~1 | Patient/Region,
correlation=corAR1(form=~Region | Patient/Region), data=df)
You're basically stuck here, the constraint on gls models is
as suggested by the error message -- covariate values must be
unique within grouping factors.
By the way, there are several other potential issues here (some
of which I think I've pointed out previously)
* I don't remember whether you said you really only have two
patients; in that case treating patient as a random effect is
unlikely to work very well.
* In any case, including a categorical predictor (such as Patient)
as both a random and a fixed effect probably leads to an
unidentifiable model (intervals(m0) indicates an non-positive-definite
variance-covariance matrix, a sign of trouble)
* ideally I would probably prefer to use a random effect of
interaction(Patient,Region), and put correlations on Region
within Patient, but that is really beyond the scope of lme.
Accepting a cruder model, or finding a local statistical
consultant, may be your choices ...
More information about the R-sig-mixed-models
mailing list