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

sxpgvts3  于 2023-07-31  发布在  其他
关注(0)|答案(1)|浏览(116)

我很乐意在metafor中使用rma.mv()的复杂模型,为此有emmprep()函数。然而,我的模型是不平衡和嵌套的,导致冗余的预测因子由于秩不足而从模型中删除。在这种情况下,emmprep()不工作。根据帮助文件:在这种情况下,在使用此功能之前,应删除原始模型中的所有冗余。
我该怎么做?我是否必须使用虚拟变量来避免秩不足,或者有更简单的解决方案?
下面是一个可重复的例子:

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.我以前问过这个问题,为了避免混淆,我已经删除了那个版本)

juzqafwq

juzqafwq1#

可能的解决方法

看起来你并不需要emmprep()那么多。emmeans中的qdrg()函数几乎是开箱即用的。唯一的问题是rma模型将拦截命名为intrcpt而不是(Intercept),我们必须解决这个问题。下面是我为feature request创建的一个数据集的示例,该数据集是我在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)   # 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

字符串
创建于2023-07-22带有reprex v2.0.2

相关问题