使用 multcomp 在双向 ANOVA 中测试对比度

机器算法验证 r 方差分析 对比
2022-03-26 13:00:35

我正在尝试使用 multcomp 包比较一个因子的两个级别与另一个级别的聚合。

我使用 afex 包中包含的数据集 obk.long 的“主题间版本”，我计算如下：

library(dplyr)
library(afex)
data(obk.long)

obk <- obk.long%>%
group_by(id)%>%
mutate(treatment,value = mean(value))%>%
distinct(id)

obk <- data.frame(obk)

现在可以通过执行以下操作来测试几个对比。

fit <- lm(value~treatment*gender,data=obk)
summary(fit)

# test for difference in A and B when gender = F
K <- matrix(c(0, 1,-1,0,0,0),1)
t <- glht(fit, linfct = K)
summary(t)

到目前为止一切都很好，但是如果我想比较 A 和 B 在性别水平上的聚合呢？我正在考虑这种方法：

K <- matrix(c(0, 0.5,-0.5,0,0.5,-0.5),1)
t <- glht(fit, linfct = K)
summary(t)

或者，我可以安装另一个没有性别的模型

fit2 <- lm(value~treatment,data=obk)
summary(fit2)

K <- matrix(c(0, 1,-1),1)
t <- glht(fit2, linfct = K)
summary(t)

但结果不一样。这可能是由于交互？

使用 multcomp 测试这种对比度的正确（最佳）方法是什么？

1个回答

这可以通过使用 with 的巧妙组合来解决afex（如果有人愿意，lsmeans也可以multcomp这样做，但这通常不是必需的）。此外，由于afex具有自动聚合的功能，dplyr因此不需要。

library(afex)
require(lsmeans)
require(multcomp)
data(obk.long)

# Step 1: set up the model using afex
# but use return = "aov" to obtain an object lsmeans can handle.
fit <- aov.car(value~treatment*gender + Error(id),data=obk.long, return = "aov")


# Step 2: set up reference grid on which we 
# we can test any type of tests using lsmeans functionality:
(ref1 <- lsmeans(fit, c("treatment", "gender")))
##  treatment gender   lsmean        SE df lower.CL upper.CL
##  control   F      4.333333 0.8718860 10 2.390650 6.276016
##  A         F      4.500000 0.8718860 10 2.557317 6.442683
##  B         F      5.833333 0.6165165 10 4.459649 7.207018
##  control   M      4.111111 0.7118919 10 2.524917 5.697305
##  A         M      8.000000 0.8718860 10 6.057317 9.942683
##  B         M      6.222222 0.7118919 10 4.636028 7.808416
## 
## Confidence level used: 0.95 

# we simply define the contrasts as a list on the reference grid:
c_list <- list(c1 = c(0, -1, 1, 0, 0, 0),
               c2 = c(0, -0.5, 0.5, 0, -0.5, 0.5))

# because we want to control for Type I errors we test this using
# the Bonferroni-Holm correction
summary(contrast(ref1, c_list), adjust = "holm")
##  contrast   estimate        SE df t.ratio p.value
##  c1        1.3333333 1.0678379 10   1.249  0.4805
##  c2       -0.2222222 0.7757662 10  -0.286  0.7804
## 
## P value adjustment: holm method for 2 tests 

# alternatively, we can pass it to multcomp for even cooler corrections:
summary(as.glht(contrast(ref1, c_list)), test = adjusted("free"))
## Note: df set to 10
## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Linear Hypotheses:
##         Estimate Std. Error t value Pr(>|t|)
## c1 == 0   1.3333     1.0678   1.249    0.359
## c2 == 0  -0.2222     0.7758  -0.286    0.780
## (Adjusted p values reported -- free method)

请注意，第一个对比（c1当性别 = F 时：A 和 B）对应于您的第一个结果，而后一个则不对应。你一定犯了一个错误。

您希望的第二个对比也可以更容易获得（只需忽略其他两行）：

pairs(lsmeans(fit, "treatment"), adjust = "none")
## NOTE: Results may be misleading due to involvement in interactions
##  contrast      estimate        SE df t.ratio p.value
##  control - A -2.0277778 0.8347673 10  -2.429  0.0355
##  control - B -1.8055556 0.7338014 10  -2.461  0.0336
##  A - B        0.2222222 0.7757662 10   0.286  0.7804
## 
## Results are averaged over the levels of: gender
其它你可能感兴趣的问题
上一篇二项分布的 Berry-Esseen 界 下一篇可视化 SVM 结果