从 R 中的 emmeans 中删除一个对比
Posted
技术标签:
【中文标题】从 R 中的 emmeans 中删除一个对比【英文标题】:Remove one contrast from emmeans in R 【发布时间】:2020-06-09 07:21:12 【问题描述】:在阅读了 emmeans 的小插曲后,我仍在为可能有一个非常简单的解决方案而苦苦挣扎。
我模拟了 2 组 6 名受试者的一些数据。测量时间长达 360 分钟 (ExpDelta)。我的 lme 模型如下:
library(lme4)
library(emmeans)
lme.model = lmer(Value ~ Treatment*ExpDelta + Baseline + (1 | SubjectNr), data = df)
现在,我可以计算每个时间点的 emmeans 对比度:
emm.s <- emmeans(lme.model, pairwise ~ Treatment | ExpDelta) # emmeans for every time point
或仅通过治疗:
emm.s <- emmeans(lme.model, 'Treatment') # emmeans over the whole investigation period
pairwise_emm<-pairs(emm.s)
两个结果都符合预期。但现在我只想比较 2 个处理组,而排除 ExpDelta 240 和 360 组,我不知道如何。
所以我的问题是:排除 ExpDelta 240 和 360 的数据时,安慰剂与 1 mg 药物 Y 的 p 值是多少?
参考数据集如下:
df<-structure(list(SubjectNr = c(1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L,
2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 3L, 4L, 4L, 4L, 4L, 4L, 4L,
5L, 5L, 5L, 5L, 5L, 5L, 6L, 6L, 6L, 6L, 6L, 6L, 7L, 7L, 7L, 7L,
7L, 7L, 8L, 8L, 8L, 8L, 8L, 8L, 9L, 9L, 9L, 9L, 9L, 9L, 10L,
10L, 10L, 10L, 10L, 10L, 11L, 11L, 11L, 11L, 11L, 11L, 12L, 12L,
12L, 12L, 12L, 12L), Treatment = structure(c(2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L), .Label = c("Placebo", "1 mg drug Y"), class = "factor"),
ExpDelta = c("30", "60", "90", "120", "240", "360", "30",
"60", "90", "120", "240", "360", "30", "60", "90", "120",
"240", "360", "30", "60", "90", "120", "240", "360", "30",
"60", "90", "120", "240", "360", "30", "60", "90", "120",
"240", "360", "30", "60", "90", "120", "240", "360", "30",
"60", "90", "120", "240", "360", "30", "60", "90", "120",
"240", "360", "30", "60", "90", "120", "240", "360", "30",
"60", "90", "120", "240", "360", "30", "60", "90", "120",
"240", "360"), Baseline = c(9.64, 9.64, 9.64, 9.64, 9.64,
9.64, 7.92, 7.92, 7.92, 7.92, 7.92, 7.92, 5.88, 5.88, 5.88,
5.88, 5.88, 5.88, 11.79, 11.79, 11.79, 11.79, 11.79, 11.79,
11.07, 11.07, 11.07, 11.07, 11.07, 11.07, 9.38, 9.38, 9.38,
9.38, 9.38, 9.38, 12.37, 12.37, 12.37, 12.37, 12.37, 12.37,
8.51, 8.51, 8.51, 8.51, 8.51, 8.51, 10.86, 10.86, 10.86,
10.86, 10.86, 10.86, 8.13, 8.13, 8.13, 8.13, 8.13, 8.13,
11.79, 11.79, 11.79, 11.79, 11.79, 11.79, 9.3, 9.3, 9.3,
9.3, 9.3, 9.3), Value = c(10.72, 11.58, 11.3, 11.28, 10.39,
10.09, 8.78, 10.71, 11.01, 9.98, 8.15, 7.85, 6.6, 8.65, 7.86,
7.7, 6.61, 6.88, 12.91, 13.3, 14.13, 14.57, 12.31, 11.02,
10.78, 12.93, 13.07, 12.07, 11.92, 11.8, 10.62, 10.62, 12.26,
11.7, 10.86, 8.97, 13.03, 12.86, 13.5, 11.45, 12.78, 12.7,
9.14, 9.08, 7.81, 8.56, 8.51, 7.73, 10.86, 11.25, 11.5, 11.21,
10.6, 11.59, 8.57, 7.54, 7.87, 8.07, 7.56, 8.7, 11.46, 11.33,
12.1, 12.18, 11.69, 11.53, 9.73, 10.01, 8.85, 9.91, 10.02,
9.01)), row.names = c(NA, -72L), class = "data.frame")
【问题讨论】:
你的模型适合整个数据集;因此,如果您真的想排除这两个处理级别的 数据,如问题中所述,那么您需要拟合排除这些处理级别的不同模型(例如,使用subset
在模型拟合阶段)。
但是,如果您只是不想考虑某些级别,您可以在调用 contrast()‘ or
pairs() 时使用 exclude‘ or
include 参数。有关详细信息,请参阅pairwise.emmc
的文档。或者在 emmeans() 调用中使用at
(参见“ref_grid”的文档)
【参考方案1】:
我认为您的问题最有可能的解释是
emm <- emmeans(lme.model, "Treatment",
at = list(ExpDelta = c("30", "60", "90", "120")))
pairs(emm)
有关at
参数的详细信息,请参阅? ref_grid
。
【讨论】:
感谢您的回答。这给了我以下错误。Error in ref_grid(object, ...) : Non-conformable elements in reference grid. Probably due to rank deficiency not handled as expected.
。我需要将 ExpDelta 更改为一个因子才能正常工作!
嗯。好吧,我试过了,得到了同样的错误。也许它与拟合模型时的消息有关:boundary (singular) fit: see ?isSingular。我以前从未见过这种情况;看看我能不能弄明白,但不要屏住呼吸,因为我不知道这需要多长时间
糟糕,我错了。我越过了所有的 t,发现它确实可以使用 ExpDelta
作为一个因素。这将帮助我解决这个问题。
我发现并修复了这个错误。请参阅 githib 站点上的问题 #175 -- github.com/rvlenth/emmeans/issues/175。固定版本现在可以从 githuib 获得,并且很快会在下一次更新到 CRAN。以上是关于从 R 中的 emmeans 中删除一个对比的主要内容,如果未能解决你的问题,请参考以下文章
如何在 R 中的 emmeans() 命令中评估字符串变量作为因子?
r包emmeans中的weights="cell"和weights="proportional"有啥区别