英文:
How to extract metafor::rma.mv() equivalent of lme4::VarCorr() for multilevel mixed effect meta-regression with random slopes
问题
我正在尝试计算多层混合效应元回归(包括随机斜率)的伪R平方,使用metafor
包中的类似方法(即rma.mv()
对象),该方法与以下引文类似:
https://besjournals.onlinelibrary.wiley.com/doi/10.1111/2041-210X.12225
我正在研究这篇引文中的代码,该代码使用lme4
/nlme
包演示了如何进行计算,似乎需要一个metafor
等效的lme4
/nlme
拟合模型的VarCorr()
输出摘要。
我正在努力弄清楚如何获得rma.mv()
对象中每个随机效应层级的方差(或标准差)的摘要。在这方面,任何见解都将不胜感激!
这是引文中用于计算伪R平方的代码:
# 拟合模型
orangemod.rs <-
lmer(circumference ~ ageYears + (ageYears | Tree), data = Orange)
# 首先我们需要设计矩阵(X)、观测数量(n)和固定效应估计值(Beta)
X <- model.matrix(orangemod.rs)
n <- nrow(X)
Beta <- fixef(orangemod.rs)
# 首先计算固定效应的方差(Nakagawa & Schielzeth的方程27):
Sf <- var(X %*% Beta)
# 其次,计算随机效应的协方差矩阵列表。
# 这里列表只有一个元素,因为只有一个随机效应层级。
(Sigma.list <- VarCorr(orangemod.rs))
# 使用论文中的方程11估计随机效应的方差分量。
# 使用sapply确保在存在多个随机效应时(即length(Sigma.list) > 1)也能正常工作。
Sl <-
sum(
sapply(Sigma.list,
function(Sigma)
{
Z <-X[,rownames(Sigma)]
sum(diag(Z %*% Sigma %*% t(Z)))/n
}))
# 由于该模型是LMM,加性离散方差Se等于残差方差。残差标准差存储为Sigma.list的属性:
Se <- attr(Sigma.list, "sc")^2
# 最后,LMM的分布特定方差Sd为零:
Sd <- 0
# 使用Nakagawa & Schielzeth的方程29和30估计边际和条件R平方:
total.var <- Sf + Sl + Se + Sd
(Rsq.m <- Sf / total.var)
(Rsq.c <- (Sf + Sl) / total.var)
除了使用VarCorr()函数的步骤之外,所有这些步骤都可以轻松地用于rma.mv()
对象:
# 其次,计算随机效应的协方差矩阵列表。
# 这里列表只有一个元素,因为只有一个随机效应层级。
(Sigma.list <- VarCorr(orangemod.rs))
# Groups Name Std.Dev. Corr
# Tree (Intercept) 1.8966
# ageYears 8.7757 -1.000
# Residual 10.0062
这是一个简单的模型,用于说明这一点,并突出了rma.mv()模型的两个可能有必要信息的组成部分,但我不确定如何将其转换为最终的VarCorr()格式。
library(metafor)
library(tidyverse)
data("dat.konstantopoulos2011")
df=dat.konstantopoulos2011%>%
as.data.frame()
# 拟合随机斜率模型
m1=rma.mv(yi ~ year, vi, random = list(~year|study, ~year|district, ~1|school), data = df,
cvvc = TRUE)
lme4::VarCorr(m1)
nlme::VarCorr(m1)
# 可能有用的组成部分
m1$V
m1$vvc
我希望重新创建的内容
# Groups Name Std.Dev. Corr
# study (Intercept) [VALUE]
# ageYears [VALUE] [VALUE]
# district (Intercept)
# ageYears [VALUE]
# school (Intercept)
# Residual [VALUE]
英文:
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]
答案1
得分: 0
~year|study
对应的随机效应的方差-协方差矩阵可以通过m1$G
提取。如果只需要方差,可以使用m1$tau2
。~year|district
的方差-协方差矩阵是m1$H
,方差在m1$gamma2
中。~1|school
对应的方差是m1$sigma2
。
请注意,如果您想要随机斜率,那么您正在拟合的模型不提供该功能。您需要使用以下代码:
rma.mv(yi ~ year, vi,
random = list(~year|study, ~year|district, ~1|school),
struct=c("GEN","GEN"), data = df)
参考链接:https://wviechtb.github.io/metafor/reference/rma.mv.html#specifying-random-effects
英文:
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
通过集体智慧和协作来改善编程学习和解决问题的方式。致力于成为全球开发者共同参与的知识库,让每个人都能够通过互相帮助和分享经验来进步。
评论