[R-sig-ME] getVarCov for lmer()
David Afshartous
dafshartous at med.miami.edu
Fri Feb 1 20:55:16 CET 2008
RE what I'd like to calculate, I'd like to calculate the marginal covariance
matrix, i.e., tr(Z_i) Psi Z_I + sigma^2 I, in the notation of your book.
I've done this for a series of models estimated with lme() using
getVarCov(model, type = "marginal"), but for the last model there is an
error message and I wanted to see if this is also the case for lmer().
Here is the problem on simulated data called dat.new (formed below;
longitudinal data from a crossover design; 20 patients, 4 measurements each,
both for drug and placebo, thus n=160 observations total):
# basic random intercept model, marginal covariance matrix is compound
# symmetric as expected:
> fm <- lme( dv ~ time.num*drug, random = ~ 1 | Patient.new, data=dat.new )
> getVarCov(fm, type = "marginal")
Patient.new 1
Marginal variance covariance matrix
1 2 3 4 5 6 7 8
1 18.8570 8.5199 8.5199 8.5199 8.5199 8.5199 8.5199 8.5199
2 8.5199 18.8570 8.5199 8.5199 8.5199 8.5199 8.5199 8.5199
3 8.5199 8.5199 18.8570 8.5199 8.5199 8.5199 8.5199 8.5199
4 8.5199 8.5199 8.5199 18.8570 8.5199 8.5199 8.5199 8.5199
5 8.5199 8.5199 8.5199 8.5199 18.8570 8.5199 8.5199 8.5199
6 8.5199 8.5199 8.5199 8.5199 8.5199 18.8570 8.5199 8.5199
7 8.5199 8.5199 8.5199 8.5199 8.5199 8.5199 18.8570 8.5199
8 8.5199 8.5199 8.5199 8.5199 8.5199 8.5199 8.5199 18.8570
Standard Deviations: 4.3424 4.3424 4.3424 4.3424 4.3424 4.3424 4.3424
4.3424
## now I make random effect variance stratified by treatment.
## as expected, the 8x8 matrix has the upper left 4x4 block
## and the lower right 4x4 block represent the compound symmetric
## covariance matrix for the different treatment conditions (recall
## that the same patient will experience both since it's a crossover).
## whereas the upper right and lower left 4x4 blocks represent
## covariance of measurements across treatment condition, e.g.,
## the covariance of measurement i on treatment to measurement j on control.
## such a covariance exists since although we have stratified the intercept
## random effect by treatment, a covariance of these random variables is
specified by default ( check via summary(fm1) ).
> fm1 <- lme( dv ~ time.num*drug, random = ~ 0 + Pind + Dind | Patient.new,
data=dat.new )
> getVarCov(fm1, type = "marginal")
Patient.new 1
Marginal variance covariance matrix
1 2 3 4 5 6 7 8
1 31.3030 29.3380 29.3380 29.3380 2.0209 2.0209 2.0209 2.0209
2 29.3380 31.3030 29.3380 29.3380 2.0209 2.0209 2.0209 2.0209
3 29.3380 29.3380 31.3030 29.3380 2.0209 2.0209 2.0209 2.0209
4 29.3380 29.3380 29.3380 31.3030 2.0209 2.0209 2.0209 2.0209
5 2.0209 2.0209 2.0209 2.0209 6.8509 4.8858 4.8858 4.8858
6 2.0209 2.0209 2.0209 2.0209 4.8858 6.8509 4.8858 4.8858
7 2.0209 2.0209 2.0209 2.0209 4.8858 4.8858 6.8509 4.8858
8 2.0209 2.0209 2.0209 2.0209 4.8858 4.8858 4.8858 6.8509
Standard Deviations: 5.5949 5.5949 5.5949 5.5949 2.6174 2.6174 2.6174
2.6174
## now, presumably the off-diagonal 4x4 blocks should go to zero if the
## covariance between the intercept random effects for treatment/control can
## be constrained to 0. This can be done via:
> fm1.no.cov <- lme( dv ~ time.num*drug, random = list(~ 0 + Pind |
Patient.new, ~ 0 + Dind | Patient.new), data=dat.new )
# But the following error message arises when trying to get the marginal
# covariance:
> getVarCov(fm1.no.cov, type = "marginal")
Error in getVarCov.lme(fm1.no.cov, type = "marginal") :
Not implemented for multiple levels of nesting
Does there exist a different model statement for fm1.no.cov such that this
problem doesn't arise? In any event, this is why I sought getVarCov in
lmer().
##################################
## code to generate simulated data above:
set.seed(500)
n.timepoints <- 4
n.subj.per.tx <- 20
sd.d <- 5;
sd.p <- 2;
sd.res <- 1.3
drug <- factor(rep(c("D", "P"), each = n.timepoints, times =
n.subj.per.tx))
drug.baseline <- rep( c(0,5), each=n.timepoints, times=n.subj.per.tx )
Patient.baseline <- rep( rnorm( n.subj.per.tx*2, sd=c(sd.d, sd.p) ),
each=n.timepoints )
time.baseline <- rep(1:n.timepoints,n.subj.per.tx*2)*as.numeric(drug=="D")
dv <- rnorm( n.subj.per.tx*n.timepoints*2,
mean=time.baseline+Patient.baseline+drug.baseline, sd=sd.res )
dat.new <- data.frame(drug, dv)
dat.new$Dind <- as.numeric(dat.new$drug == "D")
dat.new$Pind <- as.numeric(dat.new$drug == "P")
dat.new$time.num = rep(1:n.timepoints, n.subj.per.tx*2)
dat.new$Patient.new = rep(1:20, each=8)
On 2/1/08 1:29 PM, "Douglas Bates" <bates at stat.wisc.edu> wrote:
> On Feb 1, 2008 11:15 AM, David Afshartous <dafshartous at med.miami.edu> wrote:
>
>> All,
>
>> I've fitted some models using lme() and now I want to check some things by
>> fitting them in lmer() (latest version loaded from CRAN). Under the listed
>> methods for lmer() I don't seem to find functionality similar to that of
>> getVarCov() to obtain the variance-covariance matrix of the fitted model.
>> Does such functionality exist?
>
> I would need to go back and check what getVarCov produces. Can you
> tell me what you want to be able to calculate?
>
> Be aware that it is not always possible to reproduce all the
> capabilities of the nlme package in lme4 because nlme does not handle
> (except through egregious kludges) crossed or partially crossed random
> effects specifications. Some of the quantities that can be calculated
> for models with a single level of random effects or with multiple
> nested levels of random effects cannot even be defined for models with
> crossed random effects.
>
>> PS - running R v2.6.0.
>
>> PPS - if it must be computed manually, how does one access information
>> stored in slots (e.g., Z', the transpose of the model matrix for the random
>> effects)?
>
> If the fitted model is fm then Z' is available as
>
> fm at Zt
>
> It is stored as a compressed, column-oriented sparse matrix of class
> "dgCMatrix", defined in the Matrix package. Linear algebra operations
> with such objects should magically do the right thing. If you are
> working with very large models, of course, you will find that you need
> to be very careful about exactly what operations are done with Z. In
> many books on mixed models there are formulas for many calculations
> based on a generalized least squares representation of the model.
> Don't even consider using such formulas except for trivial examples.
>
> I've seen far too many papers where the authors write down an n by n
> matrix like V = sigma^2 I + Z Sigma Z' and then proceed to blithely
> describe calculations involving the inverse of V. I usually stop
> reading at that point because I know that the author has never tried
> to perform the calculation on a non-trivial example. When n is in the
> hundreds of thousands or in the millions you do not even construct n
> by n matrices, let alone try to decompose them. And anyone with
> experience in numerical linear algebra knows that you never, ever,
> ever expect to evaluate a formula that involves the inverse of a large
> matrix.
>
> This is all to say that you could write down an expression involving
> solve() applied to some matrix of size n by n derived from Z but you
> will be disappointed in the speed.
More information about the R-sig-mixed-models
mailing list