metafor
包(即
rma.mv()
对象)中的随机斜率:
https://besjournals.onlinelibrary.wiley.com/doi/10.1111/2041-210X.12225
深入研究此引文的代码,该代码使用lme4
/
nlme
包来演示如何执行计算,它看起来像是
metafor
相当于
VarCorr()
/
lme4
的
nlme
输出需要拟合模型。我正在努力弄清楚如何获得
rma.mv()
对象随机效应每个级别的方差(或 SD)摘要。任何见解都值得赞赏!这是引文中用于计算伪 R 平方的代码:
# 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)
除了使用 VarCorr() 的步骤之外,所有这些步骤都可以轻松地为 rma.mv()
对象重现,因为该函数不接受
rma.mv()
对象:
# 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
这是一个简单的模型,它说明了这一点并突出显示了 rma.mv() 模型的两个组件,我认为它们可能具有必要的信息,但我不确定如何将其转换为最终的 VarCorr() 格式.
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]
~year|study
对应的随机效应的方差-协方差矩阵可以用
m1$G
提取。如果您只需要方差,那么
m1$tau2
就可以了。
~year|district
的 var-cov 矩阵是
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://wvichtb.github.io/metafor/reference/rma.mv.html#specifying-random-effects