[R-sig-ME] getVarCov for lmer()

[David Afshartous]

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.