Skip to content
Prev 12775 / 20628 Next
Original Raw

lme with a known correlation matrix for the random effects

Viechtbauer Wolfgang (STAT) 
While this is not a meta-analysis, this model can be easily fitted with the metafor package. So-called phylogenetic meta-analyses require the specification of a known correlation matrix for a random effect term. The same principle can be used here:

library(metafor)

N <- 5
y <- c(7.38, 7.25, 7.34, 7.30, 7.06)
names(y) <- c("t4", "t2", "t5", "t1", "t3")
g <- e <- names(y)
data <- data.frame(y=y, g=g, e=e)

Psi <- matrix(c(1, 0.00, 0.00, 0.00, 0.00,
                0, 1.00, 0.28, 0.25, 0.78,
                0, 0.28, 1.00, 0.83, 0.26,
                0, 0.25, 0.83, 1.00, 0.23,
                0, 0.78, 0.26, 0.23, 1.00), nrow=N, ncol=N)

rownames(Psi) <-colnames(Psi) <- names(y)

res <- rma.mv(y, V=0, random = list(~ 1 | g, ~ 1 | e), R = list(g = Psi), data=data)
res

Note that you will get two warnings, one about non-positive sampling variances and the other about V being not positive definite. In meta-analyses, you typically have a known var-cov matrix of the sampling variances/covariances of the observed outcomes. Here, I am just setting that part of the model to 0, so it drops out. So you can ignore those warnings. The rest specifies exactly the requested model. Via the 'R' argument, one can specify a known correlation matrix for one (or more) terms. The default is REML estimation, but method="ML" will give you ML estimation.

The results look like this:

Multivariate Meta-Analysis Model (k = 5; method: REML)

Variance Components: 

            estim    sqrt  nlvls  fixed  factor    R
sigma^2.1  0.0135  0.1161      5     no       g  yes
sigma^2.2  0.0059  0.0771      5     no       e   no

Model Results:

estimate       se     zval     pval    ci.lb    ci.ub          
  7.2815   0.0792  91.8996   <.0001   7.1262   7.4368      *** 

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

So, sigma^2.1 is sigma^2_g and sigma^2.2 is sigma^2_e. To make sure that these are really good estimates of these variance components, one can profile the restricted log likelihood for each component with:

par(mfrow=c(1,2))
profile(res, sigma2=1)
profile(res, sigma2=2)

(this will generate a lot of the same warnings as mentioned above, since the model is being refitted repeatedly and each time the same issue comes up; you can ignore those again)

I hope this helps!

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 (2 messages)

Viechtbauer Wolfgang (STAT) lme with a known correlation matrix for the random effects Dec 16 Venelin Mitov lme with a known correlation matrix for the random effects Dec 17