PSYCH 710 Scheffe Test
PJB
23/09/2023
1 ANOVA
In this note I want to highlight an issue that occurs when you are
using the emmeans package to perform a Scheffe test.
I’m going to illustrate the issue with an analysis on the following data:
options(contrasts=c("contr.sum","contr.poly")) # IMPORTANT
options(digits=6)
a <- 5
n <- 10
N <- a*n
y.bar <- 0.9*c(-1/2,-1/2,0,0,1)
set.seed(143)
score <- 100+5*(rep(x=y.bar,each=n) + rnorm(N,0,1))
group <- factor(x=rep(x=seq(1,a),each=n),labels="g",ordered=F)
df0 <- data.frame(group,score)
summary(df0)
## group score
## g1:10 Min. : 88.1
## g2:10 1st Qu.: 97.0
## g3:10 Median :101.1
## g4:10 Mean :101.2
## g5:10 3rd Qu.:105.7
## Max. :112.8boxplot(score~group,data=df0)
First I will do a plain-vanilla ANOVA.
aov.01 <- aov(score~group,data=df0)
summary(aov.01)
## Df Sum Sq Mean Sq F value Pr(>F)
## group 4 200 50.1 1.63 0.18
## Residuals 45 1382 30.72 Four planned contrasts
Next, we will perform four orthogonal contrasts using the following weights.
c1 <- c(-1/4,-1/4,-1/4,-1/4,1)
c2 <- c(-1/2,-1/2,1/2,1/2,0)
c3 <- c(-1,1,0,0,0)
c4 <- c(0,0,-1,1,0)
library(MASS)
cMat <- ginv(rbind(c1,c2,c3,c4)) # our coding matrixFirst, I will use the aov and split
commands.
contrasts(df0$group) <- cMat
contrasts(df0$group) # our coding matrix
## [,1] [,2] [,3] [,4]
## g1 -0.2 -5.00000e-01 -5.00000e-01 1.11041e-17
## g2 -0.2 -5.00000e-01 5.00000e-01 3.63225e-17
## g3 -0.2 5.00000e-01 -7.10914e-17 -5.00000e-01
## g4 -0.2 5.00000e-01 -8.11958e-17 5.00000e-01
## g5 0.8 1.35912e-17 -5.34703e-17 -2.80000e-17
aov.02 <- aov(score~group,data=df0) # re-do anova
summary(aov.02,split=list(group=list(c1=1,c2=2,c3=3,c4=4)))
## Df Sum Sq Mean Sq F value Pr(>F)
## group 4 200 50.1 1.63 0.183
## group: c1 1 152 151.8 4.94 0.031 *
## group: c2 1 44 43.6 1.42 0.240
## group: c3 1 5 4.9 0.16 0.691
## group: c4 1 0 0.0 0.00 0.982
## Residuals 45 1382 30.7
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
coef(aov.02)
## (Intercept) group1 group2 group3 group4
## 101.1590132 4.3553954 2.0878661 -0.9926173 0.0577542Next, I will perform the contrasts using emmeans. Notice
that the p values are the same as those in the ANOVA. If the contrasts
were planned, then only contrast c1 would be significant.
Also note that the values of the contrasts (in the estimate
column) are the same as the coefficients/alpha’s from our linear
model.
library(emmeans)
aov.em <- emmeans(aov.02,specs="group")
con <- contrast(aov.em,method=list(c1,c2,c3,c4))
summary(con,adjust="none")3 Four post-hoc contrasts
Now let’s consider the case where the contrasts were performed after looking at the data. The Scheffe test can be performed by recalculating the critical value of F that is used to reject the null hypothesis. Here is the critical value of F used in the omnibus test:
a <- 5
n <- 10
N <- a*n
alpha <- .05
( F.omnibus <- qf(1-alpha,df1=(a-1),df2=N-a) )
## [1] 2.57874The critical value of F for the Scheffe test is obtained by
multiplying F.omnibus by the degrees of freedom for
group:
( F.Scheffe <- F.omnibus * (a-1))
## [1] 10.315If c1 was post-hoc, and we used the Scheffe test to
control Type I error at \(\alpha=0.05\), then contrast
c1 would not be significant because the observed F for that
contrast, 4.94, is less than the Scheffe critical F of 10.31.
Next we repeat the Scheffe test using the emmeans
package. The commands in emmeans perform the Scheffe test
by adjusting the p values, rather than adjusting the critical value of t
or F used to reject the null hypothesis. Nevertheless, the conclusions
about the post-hoc contrasts are the same: none of the contrasts are
significant when we use the Scheffe adjustment.
library(emmeans)
aov.em <- emmeans(aov.02,specs="group")
contrast(aov.em,method=list(c1,c2,c3,c4),adjust="scheffe")
## contrast estimate SE df t.ratio p.value
## c(-0.25, -0.25, -0.25, -0.25, 1) 4.3554 1.96 45 2.223 0.3096
## c(-0.5, -0.5, 0.5, 0.5, 0) 2.0879 1.75 45 1.191 0.8394
## c(-1, 1, 0, 0, 0) -0.9926 2.48 45 -0.400 0.9968
## c(0, 0, -1, 1, 0) 0.0578 2.48 45 0.023 1.0000
##
## P value adjustment: scheffe method with rank 4I want to call your attention to the last line in the output, which
says the p values have been adjusted with the
"scheffe method with rank 4". We will return to this point
in a moment.
4 One post-hoc contrast
In the previous example we used the Scheffe test to evaluate three
post-hoc contrasts. However, let’s suppose we had conducted only the
first contrast (c1) and never considered the other two. Our
analysis using the aov split command would
have been something like the following:
contrasts(df0$group) <- c1 # define our 1 contrast
aov.03 <- aov(score~group,data=df0) # re-do anova
summary(aov.03,split=list(group=list(c1=1)))
## Df Sum Sq Mean Sq F value Pr(>F)
## group 4 200 50.1 1.63 0.183
## group: c1 1 152 151.8 4.94 0.031 *
## Residuals 45 1382 30.7
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1The F and p values for this contrast are exactly the same as before, as they should be, and the critical Scheffe F also remains unchanged. Not surprisingly, the results of our Scheffe test are the same: the contrast is not significant. This result makes sense because the Scheffe adjustment does not depend on the number of contrasts that we perform. Instead, we assume that our inspection of the data allowed us to do many contrasts implicitly, and we selected the one that we thought had a chance of being significant.
Now let’s try it with the emmeans commands:
library(emmeans)
aov.em <- emmeans(aov.03,specs="group")
contrast(aov.em,method=list(c1),adjust="scheffe")
## contrast estimate SE df t.ratio p.value
## c(-0.25, -0.25, -0.25, -0.25, 1) 4.36 1.96 45 2.223 0.0313
##
## P value adjustment: scheffe method with rank 1Surprise! The contrast is significant using the Scheffe
adjustment. What happened? Notice the last line in the output: the p
values have been adjusted with the
"scheffe method with rank 1", not rank 4 as before. It
seems that the Scheffe adjustment used by contrast is
adjusting for multiple comparisons, rather than performing a traditional
Scheffe test. In the previous case we had four p values and here we have
only one. Hence, in this case, the p value has not been adjusted at all.
This point can be illustrated by comparing the output to the results
obtained with no adjustment:
contrast(aov.em,method=list(c1),adjust="none")
## contrast estimate SE df t.ratio p.value
## c(-0.25, -0.25, -0.25, -0.25, 1) 4.36 1.96 45 2.223 0.0313The p value is the same as the one obtained with the
scheffe adjustment. Clearly, this is not what we want to
happen because the Scheffe test should control Type I error in the case
where we looked at the data and implicitly performed many
contrasts. This point was discussed by the author of the
emmeans package here.
Fortunately, there is a simple way of addressing this issue. In the
following code, I save the results of the contrast in the variable
con and list the results with the summary
command which includes the parameter scheffe.rank. By
setting scheffe.rank to 4, which corresponds to the degrees
of freedom of group, I am telling emmeans to
provide the full Scheffe adjustment.
con <- contrast(aov.em,method=list(c1),adjust="scheffe")
summary(con,adjust="scheffe",scheffe.rank=4) # set rank to df for group!Note that the p value is the same as the one obtained when we
performed the Scheffe test with four orthogonal contrasts. So this is
the correct way to perform the post-hoc Scheffe adjustment.
Unfortunately, the scheffe.rank parameter works in the
summary command but not in the contrast
command, at least with emmeans version 1.8.8.
In conclusion, the correct way to the Scheffe test with
emmeans is to set scheffe.rank to equal the
degrees of freedom of the grouping variable, even if you perform only
one post-hoc contrast.