与手动计算相比,R 中的 glht 函数给出不同的结果
glht function in R giving different results compared to calculations done by hand
我正在编写示例并与手动获得的结果进行比较,但结果并不一致。该数据集是关于一项比较 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)
我们可以得到组均值和方差分析如下:
> 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
其中 !=
表示 "not equal to"。我们可以将假设重写为:
# 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
但是我使用 multcomp
库中的 glht()
函数没有得到这个结果:
> # 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,而不是我手算gamma-hat
时得到的值3.290862。如果需要,我可以证明获得的标准误差也不同于手工计算的标准误差。
我在不同的数据集上使用了相同的技术,并且使用手算和使用 glht
结果是一致的,所以我不确定我的错误在哪里。
有人可以帮我弄清楚我的代码或计算哪里出了问题吗?
14.005 是正确的。如果您查看 anova 模型的拟合度,您没有包括截距,因此将 T1 作为参考,系数反映了不同组偏离参考均值 (T1)
的程度
coefficients(anova.model)
它returns
(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)
我正在编写示例并与手动获得的结果进行比较,但结果并不一致。该数据集是关于一项比较 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)
我们可以得到组均值和方差分析如下:
> 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
其中 !=
表示 "not equal to"。我们可以将假设重写为:
# 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
但是我使用 multcomp
库中的 glht()
函数没有得到这个结果:
> # 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,而不是我手算gamma-hat
时得到的值3.290862。如果需要,我可以证明获得的标准误差也不同于手工计算的标准误差。
我在不同的数据集上使用了相同的技术,并且使用手算和使用 glht
结果是一致的,所以我不确定我的错误在哪里。
有人可以帮我弄清楚我的代码或计算哪里出了问题吗?
14.005 是正确的。如果您查看 anova 模型的拟合度,您没有包括截距,因此将 T1 作为参考,系数反映了不同组偏离参考均值 (T1)
的程度coefficients(anova.model)
它returns
(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)