Dear All,
I wish to model random effects that have known between-group covariance
structure using the lme() function from library nlme. However, I have yet
to get even a simple example to work. No doubt this is because I am
confusing my syntax, but I would appreciate any guidance as to how. I have
studied Pinheiro & Bates carefully (though it's always possible I've
missed something), the few posts mentioning pdSymm (some of which suggest
lme is suboptimal here anyway) and ?pdSymm (which has only a trivial
example, see later) but have not yet found a successful example of syntax
for this particular problem.
I am using the pdSymm class to specify a positive definite matrix
corresponding to the covariance structure of a random batch effect, and
passing this to lme() through the random= argument. To do this, I must
set the value= argument of pdSymm.
Consider the following simple and self-contained example:
library(nlme)
# make response and batch data
batch.names <- c("A", "B", "C")
data.df <- data.frame(
response = rnorm(100),
batch = factor(sample(batch.names, 100, replace=T))
)
# make covariance matrix for batch
batch.mat <- matrix(c(1,.5,.2, .5, 1, .3, .2, .3, 1), ncol=3)
colnames(batch.mat) <- batch.names
rownames(batch.mat) <- batch.names
# fit batch as a simple random intercept
lme(response ~ 1, data=data.df, random=~1|batch)
# ...works fine
# do the same using pdSymm notation
lme(response ~ 1, data=data.df,
random=list( batch=pdSymm(form=~1) )
)
# ...works fine also
# specify cov structure using value arg
lme(response ~ 1, data=data.df,
random=list( batch=pdSymm(
value=batch.mat,
form=~1,
nam=batch.names)
)
)
# throws error below
---snip---
Error in "Names<-.pdMat"(`*tmp*`, value = "(Intercept)")
:
Length of names should be 3
> traceback()
7: stop(paste("Length of names should be", length(dn)))
6: "Names<-.pdMat"(`*tmp*`, value = "(Intercept)")
5: "Names<-"(`*tmp*`, value = "(Intercept)")
4: "Names<-.reStruct"(`*tmp*`, value = list(batch =
"(Intercept)"))
3: "Names<-"(`*tmp*`, value = list(batch =
"(Intercept)"))
2: lme.formula(response ~ 1, data = data.df, random = list(batch =
pdSymm(value = batch.mat,
form = ~1, nam = batch.names)))
1: lme(response ~ 1, data = data.df, random = list(batch = pdSymm(value =
batch.mat,
form = ~1, nam = batch.names)))
---snip---
The length of batch.names is 3, so I find this error enigmatic. Note that
I had to specify all three of value, form and nam otherwise I got
missing args errors. Also note that doing
pdSymm(value=batch.mat, form=~1, nam=batch.names)
on the command line, like the similar invocation described on ?pdSymm,
works fine also. It's just lme() that doesn't like it.
Can anybody show me what I should be doing instead? Some successful code
will greatly clarify the issue. (My version details are below). Also, I
notice the pdMat scheme is absent from lme() in lme4. Is this
functionality deprecated in lme4 and excluded from lmer?
Many thanks,
William
Version details: running R 2.1.0 on windows XP, using nlme 3.1-57.
=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-Dr William Valdar
++44 (0)1865 287 717
Wellcome Trust Centre valdar at well.ox.ac.uk
for Human Genetics, Oxford www.well.ox.ac.uk/~valdar
William Valdar wrote:> Dear All, > > I wish to model random effects that have known between-group covariance > structure using the lme() function from library nlme. However, I have > yet to get even a simple example to work. No doubt this is because I am > confusing my syntax, but I would appreciate any guidance as to how. I > have studied Pinheiro & Bates carefully (though it's always possible > I've missed something), the few posts mentioning pdSymm (some of which > suggest lme is suboptimal here anyway) and ?pdSymm (which has only a > trivial example, see later) but have not yet found a successful example > of syntax for this particular problem. > > I am using the pdSymm class to specify a positive definite matrix > corresponding to the covariance structure of a random batch effect, and > passing this to lme() through the random= argument. To do this, I must > set the value= argument of pdSymm. > > Consider the following simple and self-contained example: > > library(nlme) > > # make response and batch data > batch.names <- c("A", "B", "C") > data.df <- data.frame( > response = rnorm(100), > batch = factor(sample(batch.names, 100, replace=T)) > ) > > # make covariance matrix for batch > batch.mat <- matrix(c(1,.5,.2, .5, 1, .3, .2, .3, 1), ncol=3) > colnames(batch.mat) <- batch.names > rownames(batch.mat) <- batch.names > > # fit batch as a simple random intercept > lme(response ~ 1, data=data.df, random=~1|batch) > > # ...works fine > > # do the same using pdSymm notation > lme(response ~ 1, data=data.df, > random=list( batch=pdSymm(form=~1) ) > ) > # ...works fine also > > # specify cov structure using value arg > lme(response ~ 1, data=data.df, > random=list( batch=pdSymm( > value=batch.mat, > form=~1, > nam=batch.names) > ) > ) > > # throws error below > > ---snip--- > Error in "Names<-.pdMat"(`*tmp*`, value = "(Intercept)") : > Length of names should be 3 > >> traceback() > > 7: stop(paste("Length of names should be", length(dn))) > 6: "Names<-.pdMat"(`*tmp*`, value = "(Intercept)") > 5: "Names<-"(`*tmp*`, value = "(Intercept)") > 4: "Names<-.reStruct"(`*tmp*`, value = list(batch = "(Intercept)")) > 3: "Names<-"(`*tmp*`, value = list(batch = "(Intercept)")) > 2: lme.formula(response ~ 1, data = data.df, random = list(batch > pdSymm(value = batch.mat, > form = ~1, nam = batch.names))) > 1: lme(response ~ 1, data = data.df, random = list(batch = pdSymm(value > = batch.mat, > form = ~1, nam = batch.names))) > ---snip--- > > The length of batch.names is 3, so I find this error enigmatic. Note > that I had to specify all three of value, form and nam otherwise I got > missing args errors. Also note that doing > > pdSymm(value=batch.mat, form=~1, nam=batch.names) > > on the command line, like the similar invocation described on ?pdSymm, > works fine also. It's just lme() that doesn't like it. > > Can anybody show me what I should be doing instead? Some successful code > will greatly clarify the issue. (My version details are below).I'm afraid that I don't understand what you are trying to do. With a formula of ~ 1 the pdSymm generator creates a 1x1 variance-covariance matrix, which you are initializing to a 3x3 matrix. What is batch.mat supposed to represent?> Also, I > notice the pdMat scheme is absent from lme() in lme4. Is this > functionality deprecated in lme4 and excluded from lmer?Yes, although I think you mean lmer in lme4. Because the lmer function allows multiple nested or non-nested grouping factors, the need for the pdMat classes is eliminated (or greatly reduced) and the code can be simplified considerably. There is an article in the 2005/1 issue of R News describing the use of lmer.
Maybe Matching Threads
- Function nlme::lme in Ubuntu (but not Win or OS X): "Non-positive definite approximate variance-covariance"
- lme fitted correlation of random effects: where is it?
- pdMat class in LME to mimic SAS proc mixed group option? Group-specific random slopes
- lme or nlme
- Multiple random effects inlme?