如何在R中对不平衡的嵌套rma.mv元分析模型使用emmprep?

huangapple go评论62阅读模式
英文:

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 &lt;- dat.bcg #take some example data to modify

#modify to make it work for this example (so now we have two factors)
dat$ablat&lt;-c(&quot;a&quot;,&quot;b&quot;,&quot;c&quot;,&quot;a&quot;,&quot;b&quot;,&quot;b&quot;,&quot;a&quot;,&quot;b&quot;,&quot;c&quot;,&quot;a&quot;,&quot;b&quot;,&quot;c&quot;,&quot;a&quot;)

#get the yis &amp; vis
dat &lt;- escalc(measure=&quot;RR&quot;, ai=tpos, bi=tneg, ci=cpos, di=cneg, data=dat,
              slab=paste0(author, &quot;, &quot;, year)) 

#run a nested model
res4 &lt;- rma.mv(yi, vi, mods = ~ alloc/ablat, data=dat) 
res4
 
library(emmeans) #load emmeans
sav &lt;- emmprep(res4) #this returns the redundant predictors error

#&quot;Error in emmprep(res4) : 
#    Cannot use function when some redundant predictors were dropped from the model.&quot;

(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 &#39;metafor&#39; package (version 4.2-0). For an
## introduction to the package please type: help(metafor)
library(emmeans)

dat &lt;- escalc(measure=&quot;RR&quot;, ai=tpos, bi=tneg, ci=cpos, di=cneg, data=dat.bcg)
dat &lt;- transform(dat,
                 acat = factor(ablat &lt; 40, labels = c(&quot;high&quot;, &quot;low&quot;)),
                 era = factor(1+as.integer((year-1940)/15), 
                              labels = c(&quot;early&quot;,&quot;mid&quot;,&quot;late&quot;)))

with(dat, table(acat, era))   # there&#39;s an empty cell
##       era
## acat   early mid late
##   high     3   2    1
##   low      0   2    5

mod &lt;- 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] = &quot;(Intercept)&quot;
RG = qdrg(object = mod)   # 命名参数 &#39;object&#39; 是必需的
RG
## &#39;emmGrid&#39; 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 &#39;metafor&#39; package (version 4.2-0). For an
## introduction to the package please type: help(metafor)
library(emmeans)

dat &lt;- escalc(measure=&quot;RR&quot;, ai=tpos, bi=tneg, ci=cpos, di=cneg, data=dat.bcg)
dat &lt;- transform(dat,
                 acat = factor(ablat &lt; 40, labels = c(&quot;high&quot;, &quot;low&quot;)),
                 era = factor(1+as.integer((year-1940)/15), 
                              labels = c(&quot;early&quot;,&quot;mid&quot;,&quot;late&quot;)))

with(dat, table(acat, era))   # there&#39;s an empty cell
##       era
## acat   early mid late
##   high     3   2    1
##   low      0   2    5

mod &lt;- 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] = &quot;(Intercept)&quot;
RG = qdrg(object = mod)   # named argument &#39;object&#39; is essential
RG
## &#39;emmGrid&#39; 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>

huangapple
  • 本文由 发表于 2023年7月17日 15:37:54
  • 转载请务必保留本文链接:https://go.coder-hub.com/76702365.html
匿名

发表评论

匿名网友

:?: :razz: :sad: :evil: :!: :smile: :oops: :grin: :eek: :shock: :???: :cool: :lol: :mad: :twisted: :roll: :wink: :idea: :arrow: :neutral: :cry: :mrgreen:

确定