Testing contrasts of factorial variables is notoriously difficult in R. Most of the things that go on beneath the surface in friendlier programs like SPSS must be spelled out specifically in R. This is both a bad and a good thing. Bad because it takes you a long time to learn how, but good because your understanding is vastly improved as a result.
I don't know where the Mon-Fri factor comes from. But, using the treat 1-5 levels from the original dataset supplied, I want to take issue with something @StatNoodle wrote in their (mostly excellent) answer. To do this we need to make the treatment
variable a factor and specify the contr.sum
set of contrasts.
df$treatment <- factor(df$treatment)
contrasts(df$treatment) <- contr.sum(5)
contrasts(df$treatment)
[,1] [,2] [,3] [,4]
treat1 1 0 0 0
treat2 0 1 0 0
treat3 0 0 1 0
treat4 0 0 0 1
treat5 -1 -1 -1 -1
I want to clarify that this set of contrasts, does not compare treat1 to treat5, treat2 to treat5, treat3 to treat5 and treat4 to treat5. It is deviation coding (see here) and compares the mean for each of treatments 1 to 4 to the grand mean, that is the mean of all five groups.
This is easy to show if we just get the means for each level
(sumDF <- aggregate(score ~ treatment, data = df, mean))
treatment score
1 treat1 6.130
2 treat2 6.528
3 treat3 7.966
4 treat4 8.744
5 treat5 10.566
And calculate the grand mean by taking the average of the group means
(grandMean <- mean(sumDF$score))
[1] 7.9868
Now run an lm()
to get the coefficients to compare to the descriptive means
contLM <- lm(score ~ treatment, df)
contLM
Coefficients:
(Intercept) treatment1 treatment2 treatment3 treatment4
7.9868 -1.8568 -1.4588 -0.0208 0.7572
The intercept is our grand mean, 7.9868
. The treatment1
coefficient is the average of scores for the treat1
group minus the grand mean 6.130-7.9868 = -1.8568
. The treatment2
coefficient is the average of scores for the treat2
group minus the grand mean 6.528-7.9868 = -1.4588
. The treatment3
coefficient is the average of scores for the treat3
group minus the grand mean 7.966-7.9868 = -0.0208
. The treatment4
coefficient is the average of scores for the treat4
group minus the grand mean 8.744-7.9868 = 0.7572
. There is no coefficient comparing the treat5
group to the grand mean (it is taken out of the matrix by the row of -1
s).
If you want to compare the average of each group to a reference group, use the default treatment coding...
(contLM <- lm(score ~ treatment, df, contrasts = list(treatment = contr.treatment(5))))
Coefficients:
(Intercept) treatment2 treatment3 treatment4 treatment5
6.130 0.398 1.836 2.614 4.436
...which if you do the arithmetic on the descriptive means above will show that intercept is the average of scores in treat1
, the treatment2
coefficient the difference between the mean score in treat2
and the mean score in treat1
, the treatment3
coefficient the difference between the mean score of treat3
and the mean score of treat1
etc.
Answers to Questions
If you want to perform custom (i.e. non-default) contrasts in R I would strongly recommend using a package designed specifically for such things. The lsmeans
package is excellent.
library(lsmeans)
First get a reference grid using the lsmeans()
function
(lsm <- lsmeans(contLM, ~treatment))
treatment lsmean SE df lower.CL upper.CL
treat1 6.130 0.3086882 20 5.486088 6.773912
treat2 6.528 0.3086882 20 5.884088 7.171912
treat3 7.966 0.3086882 20 7.322088 8.609912
treat4 8.744 0.3086882 20 8.100088 9.387912
treat5 10.566 0.3086882 20 9.922088 11.209912
The values in the lsmean
column of the output are covariate adjusted means for each group. The descending order of these coeffcients - treat1
first to treat5
last` - is important because we use these relative positions to specify what contrasts we wish to test.
Then specify the contrast you want in a list. We want to test the first four treatments (positions 1 to 4) to treat5
(fifth position).
contr <- list("treats1-4 vs treat5" = c(-1, -1, -1, -1, 4))
We then pass this list into the lsmeans::contrast()
function
contrast(lsm, contr)
contrast estimate SE df t.ratio p.value
treats1-4 vs treat5 12.896 1.380496 20 9.342 <.0001
This will yield the correct t ratio for the contrasts (square this value if you want the F ratio) and the correct p value. But notice that the estimate is not what we would expect if we perform the procedure manually on the raw group means, subtracting the mean for treat5
from the average mean scores across groups 1 to 4.
mean(c(6.130, 6.528, 7.966, 8.744)) - 10.566
[1] -3.224
In fact it the coefficient is out by a factor of exactly 4
4*-3.224
[1] -12.896
To get accurate estimates we need to transform the contrast coefficients so each of the two sides of the contrast sum to |1| and the two sides together sum to 0.
contr2 <- list("treats1-4 vs treat5" = c(-.25, -.25, -.25, -.25, 1))
contrast(lsm, contr2)
contrast estimate SE df t.ratio p.value
Question1: trt1-4 vs 5 3.224 0.3451239 20 9.342 <.0001
A great feature of lsmeans
is that you can pass multiple contrasts into the list
contr3 <- list("Question1: trt1-4 vs 5" = c(-.25, -.25, -.25, -.25, 1),
"Question2: trt1 vs trt5" = c(1, 0, 0, 0, -1),
"Question3: linear trend" = c(-0.2, -0.1, 0, 0.1, 0.2))
contrast(lsm, contr3)
contrast estimate SE df t.ratio p.value
Question1: trt1-4 vs 5 3.2240 0.34512389 20 9.342 <.0001
Question2: trt1 vs trt5 -4.4360 0.43655103 20 -10.161 <.0001
Question3: linear trend 1.1088 0.09761578 20 11.359 <.0001
These custom contrasts should answer all of your questions, with accurate estimates. Note the transformed coefficients for the polynomial linear trend contrasts. These sorts of transformations require all kinds of weird voodoo which mostly go on under the hood in SPSS (and which I delved into briefly a few years ago, see here) but are necessary to get accurate estimates. In this case we divide the trend coefficients by their squared sum, then multiply by the incremental change across coefficients (they increase by 1 each time), thus
c(-2, -1, 0, 1, 2)/(sum(c(-2, -1, 0, 1, 2)^2))*1
[1] -0.2 -0.1 0.0 0.1 0.2
We can verify that this estimate of rate of change across groups, 1.109
, is accurate because if we transform the treatment
factor to a numerical variable
df$treatNum <- rep(1:5, each = 5)
And then run the regression
lm(score ~ treatNum, df)
Coefficients:
(Intercept) treatNum
4.660 1.109
The coefficient is exactly the same as the estimate returned by the polynomial trend contrast in the contr
list, and also matches the rate of increase across factors we can see if we graph score against the numerical treatment variable
ggplot(df, aes(treatNum, score)) +
geom_line(stat = "summary", fun.y = "mean") +
geom_smooth(method = "lm")
UPDATE: the lsmeans
package has been superseded by the emmeans
package but they were created by the same person and do more or less the same thing.