如何在嵌套模型的emmeans中使用转换后的参考网格?

Posted

技术标签:

【中文标题】如何在嵌套模型的emmeans中使用转换后的参考网格?【英文标题】:How to use a transformed reference grid in emmeans from a nested model? 【发布时间】:2021-10-20 23:49:51 【问题描述】:

我一直在尝试使用对数转换的参考网格来获得与emmeans 的成对平均比率(遵循先前问题here 的建议解决方案)。

但是,我有一个嵌套模型,我不知道如何让函数 confintpairs 在从嵌套模型创建的对数转换参考网格上工作。这是一个使用来自emmeans messy data vignette 的嵌套示例的示例:

cows <- data.frame (
  route = factor(rep(c("injection", "oral"), c(5, 9))),
  drug = factor(rep(c("Bovineumab", "Charloisazepam", 
                      "Angustatin", "Herefordmycin", "Mollycoddle"), c(3,2,  4,2,3))),
  resp = c(34, 35, 34,   44, 43,      36, 33, 36, 32,   26, 25,   25, 24, 24)
)

cows.lm <- lm(resp ~ route + drug, data = cows)
cows.lrg <- ref_grid(cows.lm, transform="log")
#NOTE: A nesting structure was detected in the fitted model:
#    drug %in% route

confint(cows.lrg, type="response")
#Error in object@linfct[use.elts, , drop = FALSE] : 
#  (subscript) logical subscript too long

pairs(cows.lrg, type = "response", infer = c(TRUE, TRUE), adjust = "none")
#Error in x@linfct[i, , drop = FALSE] : subscript out of bounds

我做错了什么?

【问题讨论】:

什么是rlog?除非我是盲人,否则我看不出你在哪里定义它 哎呀抱歉 - 你不是盲人,但显然我是!我复制了错误的代码...现在编辑。 我同时发布了一个答案;并对其进行了更新,使其对您修改后的问题有意义。 【参考方案1】:

好的,我发现regrid() 函数中有一个错误。如果你甚至做summary(cows.lrg),你会得到一个错误。 问题是嵌套结构可能涉及参考网格中的一些“幻像”行,这些行没有被regrid() 保留。这是一种破解它的方法,但它并不漂亮:

cows.lrg@linfct = matrix(0, nrow = 10, ncol = 5)
cows.lrg@linfct[cows.rg@misc$display, ] = diag(5)
cows.lrg@misc$display = as.logical(cows.lrg@grid$.wgt.)

现在我们有

> summary(cows.lrg)
 route     drug           prediction SE df
 oral      Angustatin           3.53 NA  9
 injection Bovineumab           3.54 NA  9
 injection Charloisazepam       3.77 NA  9
 oral      Herefordmycin        3.24 NA  9
 oral      Mollycoddle          3.19 NA  9

Results are given on the log (not the response) scale.

不幸的是,我们仍然没有有效的协方差值,这使得 SE 全部为 NA

我会努力解决这个问题,并在 emmeans 的 GitHub 存储库上对其进行更新。可能需要一段时间;不幸的是,我今天早上刚刚将包的更新上传到 CRAN,根据墨菲的错误时机定律,它没有解决这个问题。

【讨论】:

谢谢!!期待 GitHub 上的更新,同时有一些小技巧有助于继续进行预测。 @corn_bunting 我创建了一个issue for this on GitHub——所以请参考。我相信我已经修复了它,您可以使用代码页面上的说明安装更新。奇怪的是,这将是对版本 1.6.3 的更新,该版本尚未出现在 CRAN 上。如果 CRAN 出于某种原因拒绝 1.6.3,我也会合并这些更改。 太棒了——我已经安装了它,它可以工作了,感谢你这么快修复它。我仍然在获取比率而不是与嵌套变量的对比差异方面存在问题(即pairs(emmeans(cows.lrg,~drug),type="r") 而不是pairs(emmeans(cows.lrg,~route),type="r"),但会将其发布到我之前更相关的问题中。可能是我的理解问题而不是用函数。 你可能想要~ drug | route 好的,我验证了type = "r" 在第二种方法中不能正常工作。下次我将更新推送到 GitHub 时会修复它。

以上是关于如何在嵌套模型的emmeans中使用转换后的参考网格?的主要内容,如果未能解决你的问题,请参考以下文章

如何在 R 中的 emmeans() 命令中评估字符串变量作为因子?

我可以在 LME 模型中使用 emmeans 吗?

我可以从 glmmTMB 的 emmeans 中获取 p 值吗?

在 emmeans 中定义对比

R中使用emmeans和geepack的每组边际均值和置信水平

使用带有 emmeans 的抽象公式