Skip to content
Prev 12918 / 20628 Next
Original Raw

Variance components analysis using a GLMM, how to insert a variance-covariance matrix in the model ?

Viechtbauer Wolfgang (STAT) 
Thanks, Jesse. I was about to (cautiously) suggest giving metafor a try.

Indeed, if you have a vector of estimates and a corresponding (approximately) known var-cov matrix, then one can tackle this from a meta-analytic perspective (which, in the end, is nothing else than two-stage multilevel modeling, maybe with a twist here and there). You may find this illustration useful:

http://www.metafor-project.org/doku.php/tips:two_stage_analysis

It compares the two-stage approach to a single mixed-effects model for longitudinal data. If one can use a single mixed-effects model for the analyses, then this is usually better, but this may not be possible in all cases and there may be circumstances where the two-stage approach has advantages.

In the present case, things are even simpler, since there is no additional grouping variable. 

You (C?lia) wrote that the response variable is a proportion, but you also wrote logit(phi_t), which sounds like a logit-transformed proportion. So, is Rcov the (estimated) var-cov matrix of the raw or logit transformed proportions?

At any rate, you would then want to use rma.mv() from the metafor package, which allows for an entire var-cov matrix of the estimates as input. So:

library(metafor)

load(url("http://mammal-research.org/data/example.RData"))

### this assumes that Rcov is the var-cov matrix of the raw proportions
dat <- data.frame(yi = marmot$estimates$estimate, time = as.numeric(as.character(marmot$estimates$time)) - 1990)

res <- rma.mv(yi, V=marmot$vc, mods = ~ time, data=dat)
res

var.components.reml(dat$yi, design = model.matrix(~ dat$time), vcv=marmot$vc)

But you probably want something like:

dat$id <- 1:nrow(dat)
res <- rma.mv(yi, V=marmot$vc, mods = ~ time, random = ~ 1 | id, data=dat)
res

plot(dat$time, dat$yi, pch=19, xaxt="n", xlab="Year", ylab="Estimate")
axis(side=1, at=dat$time, labels=dat$time+1990)
abline(res)
lines(dat$time, predict(res)$ci.lb, lty="dotted")
lines(dat$time, predict(res)$ci.ub, lty="dotted")

for a random-effects model. The fit looks pretty decent.

An issue here is that the analysis above is based on a model that assumes that the sampling distributions of the estimates can be approximated with normal distributions. With raw proportions, that may not be sensible, especially when some of the proportions are so close to 1 as in the present case.

Some form of a generalized mixed-effects model may be an alternative that is more suitable, but you have proportions and I see no information about the denominator on which those proportions are based. Also, I don't know how MARK computes that var-cov matrix (I downloaded the "gentle introduction" book from the MARK website, but when I realized it has more than 1000 pages, I quickly abandoned the idea of figuring this out). It may already be based on some normal approximation in the first place. Also, I don't see an obvious way of incorporating a known var-cov matrix into a GLMM, since the mean and variance structures of such models are typically intertwined.

Best,
Wolfgang

--   
Wolfgang Viechtbauer, Ph.D., Statistician   
Department of Psychiatry and Psychology   
School for Mental Health and Neuroscience   
Faculty of Health, Medicine, and Life Sciences   
Maastricht University, P.O. Box 616 (VIJV1)   
6200 MD Maastricht, The Netherlands   
+31 (43) 388-4170 | http://www.wvbauer.com

Thread (4 messages)

REZOUKI CELIA p1314815 Variance components analysis using a GLMM, how to insert a variance-covariance matrix in the model ? Jan 28 ONKELINX, Thierry Variance components analysis using a GLMM, how to insert a variance-covariance matrix in the model ? Jan 28 Jesse Whittington Variance components analysis using a GLMM, how to insert a variance-covariance matrix in the model ? Jan 28 Viechtbauer Wolfgang (STAT) Variance components analysis using a GLMM, how to insert a variance-covariance matrix in the model ? Jan 28