英文:
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和visdat <- 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)res4library(emmeans) #加载emmeanssav <- 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 & visdat <- escalc(measure="RR", ai=tpos, bi=tneg, ci=cpos, di=cneg, data=dat,slab=paste0(author, ", ", year))#run a nested modelres4 <- rma.mv(yi, vi, mods = ~ alloc/ablat, data=dat)res4library(emmeans) #load emmeanssav <- 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 5mod <- 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 5mod <- 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 essentialRG## '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>
通过集体智慧和协作来改善编程学习和解决问题的方式。致力于成为全球开发者共同参与的知识库，让每个人都能够通过互相帮助和分享经验来进步。
评论