英文:
How to use emmprep on an unbalanced nested rma.mv meta-analysis model in R?
问题
我想使用emmeans
来处理metafor
中从rma.mv()
得到的复杂模型,并且为此存在emmprep()
函数。然而,我拥有的模型是不平衡和嵌套的,导致由于秩亏而删除了冗余的预测变量。在这种情况下,emmprep()
不起作用。根据帮助文件的描述:“在这种情况下,应该在使用此函数之前删除原始拟合模型中的任何冗余项。”
我该如何做?我是否需要使用虚拟变量来避免秩亏,还是有更简单的解决方案?
这里有一个可重现的示例:
library(metafor)
dat <- dat.bcg #获取一些示例数据进行修改
# 修改以使其适用于此示例(现在我们有两个因子)
dat$ablat <- c("a", "b", "c", "a", "b", "b", "a", "b", "c", "a", "b", "c", "a")
# 获取yis和vis
dat <- escalc(measure = "RR", ai = tpos, bi = tneg, ci = cpos, di = cneg, data = dat,
slab = paste0(author, ", ", year))
# 运行一个嵌套模型
res4 <- rma.mv(yi, vi, mods = ~ alloc/ablat, data = dat)
res4
library(emmeans) #加载emmeans
sav <- emmprep(res4) #这会返回冗余预测变量错误
# "Error in emmprep(res4) :
# Cannot use function when some redundant predictors were dropped from the model."
英文:
I would love to use emmeans
with a complicated model from rma.mv()
in metafor
, and for this there is the function emmprep()
. However, the model I have is unbalanced and nested, leading to redundant predictors being dropped from the model due to rank deficiencies. In this case, emmprep()
does not work. As per the helpfile: "In this case, one should remove any redundancies in the original model fitted before using this function."
How do I do that? Do I have to use dummy variables to avoid rank deficiencies, or is there a simpler solution?
Here's a reproducible example:
library(metafor)
dat <- dat.bcg #take some example data to modify
#modify to make it work for this example (so now we have two factors)
dat$ablat<-c("a","b","c","a","b","b","a","b","c","a","b","c","a")
#get the yis & vis
dat <- escalc(measure="RR", ai=tpos, bi=tneg, ci=cpos, di=cneg, data=dat,
slab=paste0(author, ", ", year))
#run a nested model
res4 <- rma.mv(yi, vi, mods = ~ alloc/ablat, data=dat)
res4
library(emmeans) #load emmeans
sav <- emmprep(res4) #this returns the redundant predictors error
#"Error in emmprep(res4) :
# Cannot use function when some redundant predictors were dropped from the model."
(PS. I asked this question badly before and have deleted that version to avoid confusion)
答案1
得分: 1
可能的解决方法
似乎你不太需要 emmprep()
。emmeans 中的 qdrg()
函数几乎可以直接使用。唯一的问题是 rma
模型将截距命名为 intrcpt
而不是 (Intercept)
,我们需要修复这个问题。这里有一个我为在 metafor 的 GitHub 网站上发布的 功能请求 创建的数据集的示例
library(metafor)
## Loading required package: Matrix
## Loading required package: metadat
## Loading required package: numDeriv
##
## Loading the 'metafor' package (version 4.2-0). For an
## introduction to the package please type: help(metafor)
library(emmeans)
dat <- escalc(measure="RR", ai=tpos, bi=tneg, ci=cpos, di=cneg, data=dat.bcg)
dat <- transform(dat,
acat = factor(ablat < 40, labels = c("high", "low")),
era = factor(1+as.integer((year-1940)/15),
labels = c("early","mid","late")))
with(dat, table(acat, era)) # there's an empty cell
## era
## acat early mid late
## high 3 2 1
## low 0 2 5
mod <- rma(yi, vi, mods = ~ acat * era, data = dat)
## Warning: Redundant predictors dropped from the model.
emmprep(mod) # 失败
## Error in emmprep(mod): Cannot use function when some redundant predictors were dropped from the model.
# 不过...
rownames(mod$beta)[1] = "(Intercept)"
RG = qdrg(object = mod) # 命名参数 'object' 是必需的
RG
## 'emmGrid' object with variables:
## acat = high, low
## era = early, mid, late
## Some things are non-estimable (null space dim = 1)
emmeans(RG, ~era | acat)
## acat = high:
## era emmean SE df lower.CL upper.CL
## early -0.971 0.244 8 -1.533 -0.409
## mid -1.366 0.347 8 -2.166 -0.566
## late -1.442 0.325 8 -2.190 -0.693
##
## acat = low:
## era emmean SE df lower.CL upper.CL
## early nonEst NA NA NA NA
## mid -0.299 0.340 8 -1.082 0.485
## late -0.268 0.161 8 -0.641 0.104
##
## Confidence level used: 0.95
<sup>使用 reprex v2.0.2 在2023年07月22日创建</sup>
英文:
Possible workaround
It seems you don't really need emmprep()
all that much. The qdrg()
function in emmeans almost works out-of-the-box. The only hitch is that rma
models name the intercept intrcpt
instead of (Intercept)
, and we have to fix that. Here's an example with a dataset I created for the feature request I posted on metafor's GitHub site
library(metafor)
## Loading required package: Matrix
## Loading required package: metadat
## Loading required package: numDeriv
##
## Loading the 'metafor' package (version 4.2-0). For an
## introduction to the package please type: help(metafor)
library(emmeans)
dat <- escalc(measure="RR", ai=tpos, bi=tneg, ci=cpos, di=cneg, data=dat.bcg)
dat <- transform(dat,
acat = factor(ablat < 40, labels = c("high", "low")),
era = factor(1+as.integer((year-1940)/15),
labels = c("early","mid","late")))
with(dat, table(acat, era)) # there's an empty cell
## era
## acat early mid late
## high 3 2 1
## low 0 2 5
mod <- rma(yi, vi, mods = ~ acat * era, data = dat)
## Warning: Redundant predictors dropped from the model.
emmprep(mod) # fails
## Error in emmprep(mod): Cannot use function when some redundant predictors were dropped from the model.
# However ...
rownames(mod$beta)[1] = "(Intercept)"
RG = qdrg(object = mod) # named argument 'object' is essential
RG
## 'emmGrid' object with variables:
## acat = high, low
## era = early, mid, late
## Some things are non-estimable (null space dim = 1)
emmeans(RG, ~era | acat)
## acat = high:
## era emmean SE df lower.CL upper.CL
## early -0.971 0.244 8 -1.533 -0.409
## mid -1.366 0.347 8 -2.166 -0.566
## late -1.442 0.325 8 -2.190 -0.693
##
## acat = low:
## era emmean SE df lower.CL upper.CL
## early nonEst NA NA NA NA
## mid -0.299 0.340 8 -1.082 0.485
## late -0.268 0.161 8 -0.641 0.104
##
## Confidence level used: 0.95
<sup>Created on 2023-07-22 with reprex v2.0.2</sup>
通过集体智慧和协作来改善编程学习和解决问题的方式。致力于成为全球开发者共同参与的知识库,让每个人都能够通过互相帮助和分享经验来进步。
评论