r - 在 TukeyHSD 中显示家庭课程
问题描述
我习惯在 minitab 中进行 Tukey 事后测试。当我这样做时,我通常会得到因变量/预测变量的家庭分组。
在 R 中,TukeyHSD()
不显示(或计算?)使用族分组。它仅显示每个因/预测变量之间的关系。是否可以像在 minitab 中一样显示家庭分组?
使用diamonds
数据集:
av <- aov(price ~ cut, data = diamonds)
tk <- TukeyHSD(av, ordered = T, which = "cut")
plot(tk)
输出:
Fit: aov(formula = price ~ cut, data = diamonds)
$cut
diff lwr upr p adj
Good-Ideal 471.32248 300.28228 642.3627 0.0000000
Very Good-Ideal 524.21792 401.33117 647.1047 0.0000000
Fair-Ideal 901.21579 621.86019 1180.5714 0.0000000
Premium-Ideal 1126.71573 1008.80880 1244.6227 0.0000000
Very Good-Good 52.89544 -130.15186 235.9427 0.9341158
Fair-Good 429.89331 119.33783 740.4488 0.0014980
Premium-Good 655.39325 475.65120 835.1353 0.0000000
Fair-Very Good 376.99787 90.13360 663.8622 0.0031094
Premium-Very Good 602.49781 467.76249 737.2331 0.0000000
Premium-Fair 225.49994 -59.26664 510.2665 0.1950425
解决方案
这是有关如何为数据集重现 minitab 表的分步示例ggplot2::diamonds
。我已经尽可能地包含了细节/解释。
请注意,据我所知,minitab 表中显示的结果与 Tukey 事后检验的结果无关。它们基于方差分析的结果。Tukey 的诚实显着性差异 (HSD) 检验是一种事后检验,根据 ANOVA 结果确定哪些比较(在所有可能的组均值成对比较中)是(诚实地)具有统计显着性的。
为了重现 minitabs“均值分组”汇总表(参见 minitab Express Support 的“解释结果:步骤 3”的第一个表),我建议(重新)运行线性模型以提取均值和置信区间。请注意,这正是aov
适合每个组的方差分析模型的方式。
拟合线性模型
我们指定一个0
偏移量来获得每个组的绝对估计值(而不是相对于偏移量的变化的估计值)。
fit <- lm(price ~ 0 + cut, data = diamonds)
coef <- summary(fit)$coef;
coef;
# Estimate Std. Error t value Pr(>|t|)
#cutFair 4358.758 98.78795 44.12236 0
#cutGood 3928.864 56.59175 69.42468 0
#cutVery Good 3981.760 36.06181 110.41487 0
#cutPremium 4584.258 33.75352 135.81570 0
#cutIdeal 3457.542 27.00121 128.05137 0
确定家庭分组
为了获得类似于 minitab 的“家庭分组”的东西,我们采用以下方法:
- 计算所有参数的置信区间
- 对所有参数的置信区间数据执行层次聚类分析
- 在与 CI 的标准偏差相对应的高度切割生成的树。这将为我们提供一组基于置信区间的参数估计值。这是一种有点经验的方法,但也是合理的,因为树测量了置信区间之间的成对距离,标准差可以解释为欧几里得距离。
我们首先计算置信区间,然后使用完全链接的层次聚类对得到的距离矩阵进行聚类。
CI <- confint(fit);
hc <- hclust(dist(CI));
我们检查集群树状图
plot(hc);
我们现在在与所有参数估计中所有 CI 的标准差相对应的高度切割树,以获得“家庭分组”
grps <- cutree(hc, h = sd(CI))
总结结果
最后,我们整理所有数量并将结果存储在一个类似于 minitab 的“均值分组”表的表中。
library(tidyverse)
bind_cols(
cut = rownames(coef),
N = as.integer(table(fit$model$cut)),
Mean = coef[, 1],
Groupings = grps) %>%
as.data.frame()
# cut N Mean Groupings
#1 cutFair 1610 4358.758 1
#2 cutGood 4906 3928.864 2
#3 cutVery Good 12082 3981.760 2
#4 cutPremium 13791 4584.258 1
#5 cutIdeal 21551 3457.542 3
请注意,我们的结果与 minitab“平均分组”表中的结果几乎完全一致:单独cut = Ideal
属于组3
(C
minitab 表中的组),而Fair
+Premium
共享组1
(minitab: group A
) 和Good
+Very Good
共享组2
(minitab:组B
)。
推荐阅读
- javascript - 提交表单结果,打印到第二个 url?
- vb.net - VB .NET 检查 FTP 帐户是否已存在
- html - 将 HTML 链接到已删除的 CSS 文件的行为就像文件仍然存在一样
- javascript - DataTables 不能使用 JSON
- c++ - 如何在多个文件之间创建一个全局变量?
- javascript - Highchart 下载图像正在改变纵横比
- angular - Angular2+ 使用/存储 Web 服务凭据的最佳实践是什么?
- twig - 如何在craft cms/twig的搜索页面上显示搜索结果条目?
- html - 如何仅在绝对必要时使浏览器自动连字符?
- c++ - 两个类的前向声明会导致构造函数中的循环依赖吗?