I am coding an example and comparing with results obtained by hand however the results are not consistent. The data set is about a study comparing a score for 3 different treatments. The data set can be reproduced as follows (relabelled for convenience).
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)
We can get the group means and the ANOVA as follows:
> 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
Suppose we wanted to conduct the following hypothesis test using a contrast:
# H0 : mu_1 = (mu_2 + mu_3) / 2
# Ha : mu_1 != (mu_2 + mu_3) / 2
where != means "not equal to". We can re-write the hypotheses as:
# H0 : mu_1 - (mu_2)/2 - (mu_3)/2 = 0
# Ha : mu_1 - (mu_2)/2 - (mu_3)/2 != 0
This gives us the contrast coefficients of c1=1, c2=-1/2 and c3=-1/2
If I use the sample means to calculate the contrast gamma-hat by hand, we get
# 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
However I do not get this result using the glht() function in the multcomp library:
> # 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)
From the output obtained, we see that the estimate of gamma-hat is 14.005 and it is not 3.290862 which is the value obtained when I calculated gamma-hat by hand. I can show that the standard error obtained is also different to the standard error calculated by hand if desired.
I used the same technique on a different data set and the results were consistent when calculated by hand and using glht so I'm not sure where my error is.
Can someone please help me figure out where I am going wrong with either my code or calculation please?
