R中的glht函数与手工计算相比得出不同的结果

Question

我正在编码示例，并与手工获得的结果进行比较，但是结果不一致。数据集是关于一项研究，比较了3种不同治疗方法的得分。数据集可以如下复制（为方便起见重新标记）。

Score = c(12, 11, 15, 11, 10, 13, 10, 4, 15, 16, 9, 14, 10, 6, 10, 8, 11, 12,
          12, 8,12, 9, 11, 15, 10, 15, 9, 13, 8, 12, 10, 8, 9, 5, 12,4, 6, 9, 4, 7, 7, 7,
          9, 12, 10, 11, 6, 3, 4, 9, 12, 7, 6, 8, 12, 12, 4, 12,13, 7, 10, 13, 9, 4, 4,
          10, 15, 9,5, 8, 6, 1, 0, 8, 12, 8, 7, 7, 1, 6, 7, 7, 12, 7, 9, 7, 9, 5, 11, 9, 5,
          6, 8,8, 6, 7, 10, 9, 4, 8, 7, 3, 1, 4, 3)

Treatment = c(rep("T1",35), rep("T2",33), rep("T3",37))

medicine = data.frame(Score, Treatment)

我们可以通过以下方法获得组均值和ANOVA：

> aggregate(medicine$Score ~ medicine$Treatment, FUN = mean)
  medicine$Treatment medicine$Score
1                 T1      10.714286
2                 T2       8.333333
3                 T3       6.513514

> anova.model = aov(Score ~ Treatment, dat = medicine)
> anova(anova.model)
Analysis of Variance Table

Response: Score
           Df Sum Sq Mean Sq F value    Pr(>F)    
Treatment   2 318.51 159.255   17.51 2.902e-07 ***
Residuals 102 927.72   9.095                      
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

假设我们想使用对比进行以下假设检验：

# H0 : mu_1 = (mu_2 + mu_3) / 2
# Ha : mu_1 != (mu_2 + mu_3) / 2

其中!=表示“不等于”。我们可以将假设改写为：

# H0 : mu_1 - (mu_2)/2 - (mu_3)/2 = 0
# Ha : mu_1 - (mu_2)/2 - (mu_3)/2 != 0

这为我们提供了c1 = 1，c2 = -1 / 2和c3 = -1 / 2的对比度系数

[如果我使用样本方法手动计算对比度伽玛帽，我们得到

# gamma-hat = (c1)(x-bar_1) + (c2)(x-bar_2) + (c3)(x-bar_3)
# gamma-hat = (1)(10.714286) - (1/2)(8.333333) - (1/2)(6.513514) = 3.290862

但是我没有使用glht()库中的multcomp函数得到此结果：

> # run code above
>
> library(multcomp)
>
> # declare contrast with coefficients corresponding to those in hypothesis test
> contrast = matrix(c(1, -1/2, -1/2), nrow = 1)
>
> # anova model declared earlier
> contrast.model = glht(anova.model , linfct = contrast)
> summary(contrast.model)

     Simultaneous Tests for General Linear Hypotheses

Fit: aov(formula = Score ~ Treatment, data = medicine)

Linear Hypotheses:
       Estimate Std. Error t value Pr(>|t|)    
1 == 0   14.005      1.082   12.95   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Adjusted p values reported -- single-step method)

[从获得的输出中，我们看到gamma-hat的估计值为14.005，而不是3.290862，这是我手工计算gamma-hat时获得的值。我可以证明，如果需要，获得的标准误差也与手工计算的标准误差不同。

我在不同的数据集上使用了相同的技术，并且通过手工计算并使用glht时结果是一致的，所以我不确定我的错误在哪里。

有人可以帮助我确定我的代码或计算出什么问题吗？

Answer 1

14.005是正确的。如果您查看方差分析模型的拟合，则不包括截距，因此将T1用作参考，系数反映出不同组有多少偏离参考平均值（T1）]

coefficients(anova.model)

返回

(Intercept) TreatmentT2 TreatmentT3 
  10.714286   -2.380952   -4.200772

例如，T2的系数是-2.38，因为它的平均值是-2.38 + 10.712 = 8.3如果使用指定的对比度来计算差异：

coefficients(anova.model)%*%t(contrast)

您得到与系数（contrast.model）相同的估计。

要获得上面想要的东西，您必须做：

anova.model = aov(Score ~ 0+Treatment, dat = medicine)
contrast.model <- glht(anova.model , linfct = contrast)