r - 如何从拟合模型中获取成对 p 值表
问题描述
我在下面有一个数据dt
,我已经为这个数据拟合了模型。我想获得模型中所有可能对(所有级别)的成对 p 值表。
数据:
dt <- structure(list(treatment = structure(c(1L, 1L, 1L, 1L, 2L, 2L,
2L, 2L, 3L, 3L, 3L, 3L, 4L, 4L, 4L, 4L, 1L, 1L, 1L, 1L, 2L, 2L,
2L, 2L, 3L, 3L, 3L, 3L, 4L, 4L, 4L, 4L), .Label = c("CC", "CCS",
"CS", "SCS"), class = "factor"), block = structure(c(1L, 2L,
3L, 4L, 1L, 2L, 3L, 4L, 1L, 2L, 3L, 4L, 1L, 2L, 3L, 4L, 1L, 2L,
3L, 4L, 1L, 2L, 3L, 4L, 1L, 2L, 3L, 4L, 1L, 2L, 3L, 4L), .Label = c("1",
"2", "3", "4"), class = "factor"), yield = c(5156L, 5157L, 5551L,
5156L, 4804L, 4720L, 4757L, 5021L, 4826L, 4807L, 4475L, 4596L,
4669L, 4588L, 4542L, 4592L, 5583L, 5442L, 5693L, 5739L, 5045L,
4902L, 5006L, 5086L, 4639L, 4781L, 4934L, 4857L, 4537L, 4890L,
4842L, 4608L)), row.names = c(NA, -32L), class = "data.frame")
fit <- lmerTest::lmer(yield ~ treatment + (1|block), data = dt)
summary(fit)$coeff
预期的结果将类似于(全部填充):
treatmentCC treatmentCCS treatmentCS treatmentSCS
treatmentCC 2.44E-06 1.25E-08 1.34E-09
treatmentCCS 2.44E-06
treatmentCS 1.25E-08
treatmentSCS 1.34E-09
我可以重新调整因子并多次拟合模型,并手动填写 p 值表,但我的原始数据中有很多级别。有没有简单的方法来做到这一点?
解决方案
为什么您认为多重比较的 p 值与简单比较相同?(如果您想了解更多信息,也许 Cross Validated 是更好的询问网站)。我同意 Stéphane Laurent 的评论。
fit <- lmerTest::lmer(yield ~ treatment + (1|block), data = dt)
emmeans_mod <- emmeans::emmeans(fit, pairwise ~ treatment)
## the results of lmerTest
summary(fit) %>%
pluck("coefficients") %>%
as.data.frame()
# Estimate Std. Error df t value Pr(>|t|)
# treatmentCCS -517.000 87.73629 28 -5.892658 2.443788e-06
## the result of summary method of emmeans
summary(emmeans_mod, adjust = "none") %>% # default adjust is tukey
purrr::pluck("contrasts") %>%
as.data.frame()
# contrast estimate SE df t.ratio p.value
# 1 CC - CCS 517.000 87.73629 25 5.8926583 3.783457e-06
pt(-5.892658, df = 28, lower.tail = TRUE) * 2 # 2.44379e-06
pt(-5.892658, df = 28 - 3, lower.tail = TRUE) * 2 # 3.78346e-06
推荐阅读
- javascript - 在 React.js 中使用相同的索引来映射两个数组
- javascript - 为什么 jquery 在 bootstrap4 中不起作用?
- kotlin - Kotlin 数据类:用作改造响应时成员未初始化
- reactjs - React 将 dom 节点作为输入传递给另一个组件
- performance - 如何使用 jmeter+influxDB+Grafana Stack 在 grafana 中生成“连接时间随时间变化”图?
- phpstorm - Xdebug 已安装,但似乎无法与 PhpStorm 一起使用
- sockets - 为什么 UDP writer 拨号?
- python - JupyterLab vs. JupyterNotebook 以及如何快速编辑代码
- c++ - 如何在 C++ 中有条件地包含标头?
- laravel - 如何构建后端以最小化所需的 axios 调用数量