1

I am attempting to calculate pseudo R squared's for a multilevel mixed-effect meta-regression that includes random slopes in the metafor package (i.e., rma.mv() object) using a similar approach to the following citation:

https://besjournals.onlinelibrary.wiley.com/doi/10.1111/2041-210X.12225

Digging into the code of this citation that uses the lme4/nlme packages to demonstrate how to perform the calculations, it seems like a metafor equivalent of a VarCorr() output of a lme4/nlme fit model is needed.

I am struggling to figure out how to obtain this summary of the variances (or SD) at each level of the random effects for a rma.mv() object. Any insight here is appreciated!

This is the code from the citation used to calculate the pseudo R squared's:


# Fit Model
orangemod.rs <- 
    lmer(circumference ~ ageYears + (ageYears | Tree), data = Orange)


# First we need the design matrix (X), the number of observations (n)
# and the fixed effect estimates (Beta)

  X <- model.matrix(orangemod.rs)
  n <- nrow(X)
  Beta <- fixef(orangemod.rs)

# First the fixed effects variance (eqn 27 of Nakagawa & Schielzeth): 

  Sf <- var(X %*% Beta)

# Second, the list of covariance matrices of the random effects.
# Here the list has only a single element because there is only
# one level of random effects.

  (Sigma.list <- VarCorr(orangemod.rs))

# Use equation 11 in the paper to estimate the random effects variance 
# component.
# Using sapply ensures that the same code will work when there are multiple
# random effects (i.e. where length(Sigma.list) > 1)

  Sl <- 
    sum(
      sapply(Sigma.list,
        function(Sigma)
        {
          Z <-X[,rownames(Sigma)]
          sum(diag(Z %*% Sigma %*% t(Z)))/n
        }))

# As this model is an LMM, the additive dispersion variance, Se, is
# equivalent to the residual variance. The residual standard deviation
# is stored as an attribute of Sigma.list:
  
  Se <- attr(Sigma.list, "sc")^2

# Finally, the distribution-specific variance, Sd, is zero for LMMs:

  Sd <- 0

# Use eqns 29 and 30 from Nakagawa & Schielzeth to estimate marginal and
# conditional R-squared:

  total.var <- Sf + Sl + Se + Sd
  (Rsq.m <- Sf / total.var) 
  (Rsq.c <- (Sf + Sl) / total.var)

All of these steps can be easily reproduced for the rma.mv() object besides the step with VarCorr() as this function doesn't accept rma.mv() objects:

# Second, the list of covariance matrices of the random effects.
# Here the list has only a single element because there is only
# one level of random effects.

  (Sigma.list <- VarCorr(orangemod.rs))

# Groups   Name        Std.Dev. Corr  
# Tree     (Intercept)  1.8966        
#          ageYears     8.7757  -1.000
# Residual             10.0062

Here is a simple model that illustrates the point and highlights two of the components of the rma.mv() model that I think might have the necessary information but I'm not exactly sure how to get it into the final VarCorr() format.

library(metafor)
library(tidyverse)

data("dat.konstantopoulos2011")
df=dat.konstantopoulos2011%>%
  as.data.frame()

#fit random slope model
m1=rma.mv(yi ~ year, vi, random = list(~year|study, ~year|district, ~1|school), data = df,
          cvvc = TRUE)

lme4::VarCorr(m1)
nlme::VarCorr(m1)

#potentially useful components
m1$V
m1$vvc


What I'm looking to recreate

# Groups   Name        Std.Dev. Corr  
# study    (Intercept)  [VALUE]        
#          ageYears     [VALUE]  [VALUE]
# district (Intercept)
#          ageYears              [VALUE]
# school   (Intercept)
# Residual             [VALUE]
zr2015
  • 13
  • 3

1 Answers1

0

The variance-covariance matrix of the random effects corresponding to ~year|study can be extracted with m1$G. If you only need the variances, then m1$tau2 will do. The var-cov matrix for ~year|district is m1$H and the variances are in m1$gamma2. The variance corresponding to ~1|school is m1$sigma2.

Note that if you want random slopes, then the model you are fitting doesn't provide that. You would need to use:

rma.mv(yi ~ year, vi, 
       random = list(~year|study, ~year|district, ~1|school),
       struct=c("GEN","GEN"), data = df)

See: https://wviechtb.github.io/metafor/reference/rma.mv.html#specifying-random-effects

Wolfgang
  • 2,810
  • 2
  • 15
  • 29
  • I appreciate your response and apologies for the question on GitHub. So it looks like here if I get the diagonals of G / H, that will give me the standard deviations as in the VarCorr() table for a lmer() fit? I have edited my message to show the output I'm looking to recreate. Further, it sounds like what you're saying is that sigma2 will give me the standard deviation for the intercept of school. How can I get the standard deviation of the remaining residual error? Again, thank you for your time and expertise I sincerely appreciate it! – zr2015 Aug 17 '23 at 10:28
  • Yes, the square-root of the diagonals are the standard deviations. You can get the correlation from `cov2cor(m1$G)` and `cov2cor(m1$H)`. The sampling variances you feed into the model (via `vi`) are the variances of the 'residuals' (which are heteroscedastic, so there isn't just one variance component for 'residuals'). – Wolfgang Aug 17 '23 at 14:02