lme4::VarCorr() equivalent in rma.mv object with random slopes

[zacrobinson5]

I am attempting to calculate pseudo R squared's for a multilevel mixed-effect meta-regression that includes random slopes (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) #neither function works

#potentially useful components
m1$V
m1$vvc

I also posted about this issue on SO here.

[wviechtb]

Dear Zac,

The GitHub issues page is for discussing the development of the package (including potential bug reports), not for getting help with using the package. General questions about the use of the metafor package should not be asked here. For this, we have the r-sig-meta-analysis mailing list (https://stat.ethz.ch/mailman/listinfo/r-sig-meta-analysis) or places like Stack Overflow. I will respond there.

Best, Wolfgang