I wanted to add some more basic information to the previous (great) responses, and clarify a little (also for myself) how contrast coding works in R, and why we need to calculate the inverse of the contrast coding matrix to understand which comparisons are performed.
I'll start with the description of the linear model and contrasts in terms of matrix algebra, and then go through an example in R.
The cell means model for ANOVA is:
\begin{equation} y = X\mu + \epsilon = X\begin{pmatrix} \mu1 \\\mu2 \\\mu3 \\\mu4 \end{pmatrix} + \epsilon \end{equation}
With X as the design matrix and u as the vector of means. An example is this, where we have 4 groups coded in each column:
\begin{equation} X=\begin{pmatrix} 1 & 0 & 0 & 0 \\ 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \\ 0 & 0 & 0 & 1 \\ \end{pmatrix} \end{equation}
In this case, we can estimate the means by the least square method, using the equations:
\begin{equation} \hat{\mu} =(X^{\prime }X)^{-1}\ X^{\prime }y\\ \end{equation}
This is all good, but let's imagine we want specific comparisons rather than the means, like differences in means compared to a reference group. In the case of 4 groups, we could express this as a matrix C of comparisons, multiplied by the vector of means:
\begin{equation} C\mu = \begin{pmatrix} \phantom{..} 1 & 0 & 0 &0 \\ -1 & 1 & 0 & 0\\ -1 & 0 & 1 & 0\\ -1 & 0 & 0 & 1\\ \end{pmatrix}\ \begin{pmatrix}\mu1 \\\mu2 \\\mu3 \\\mu4 \end{pmatrix} = \begin{pmatrix} \mu1 \\\mu2-\mu1 \\\mu3-\mu1 \\\mu4-\mu1 \end{pmatrix} \end{equation}
The first group serves as reference, and we calculate the deviations from it. The matrix C serves to describe the comparisons, it acts as the original matrix of contrasts, the contrast matrix. Technically here these are not contrasts, because the sum in each row should be zero by definition, but that will serve our purpose, and this is the matrix referred to in the contr.treatment() function in R (its inverse, see below).
The matrix C defines the contrasts. We want to evaluate contrasts from the data, in the context of the same model. We note that:
\begin{equation} y \ =\ X\mu \ +\ \epsilon \ =\ X \ (C^{-1}C)\ \ \mu \ +\ \epsilon = \ (X C^{-1}) \ (C \mu) \ + \epsilon \end{equation}
Therefore we can use the first term in parentheses to evaluate the second term (our comparisons), using the least squares method, just as we did for the original equation above. This is why we use the inverse of the contrast matrix C (it is invertible in this case).
We use the least square method to evaluate the contrasts, with the same equation as above, using the modified design matrix:
\begin{equation} (X C^{-1}) \end{equation}
And we evaluate: \begin{equation} C\mu \end{equation} using the method of least squares. The coefficients for this model can be evaluated as before using least squares, replacing the original design matrix by the new one. Or working it out in details, by replacing X and u by the new matrices: \begin{equation} [(X C^{-1})'(X C^{-1})]^{-1}(X C^{-1})'y=[(C^{-1})' X'(XC^{-1})]^{-1}(X C^{-1})'y= (C^{-1})^{-1}(X'X)^{-1}[(C^{-1})']^{-1}(X C^{-1})'y= C(X'X)^{-1}[(C^{-1})^{-1}]'(XC^{-1})'y= C(X'X)^{-1}C'(C^{-1})'X'y=C(X'X)^{-1}C'(C')^{-1}X'y= C(X^{\prime }X)^{-1}\ X^{\prime }y=\\C\hat{\mu}= \begin{pmatrix} \hat{\mu1} \\\hat{\mu2}-\hat{\mu1} \\\hat{\mu3}-\hat{\mu1} \\\hat{\mu4}-\hat{\mu1} \end{pmatrix} \end{equation}
(cf the least square equation above for the last step)
As expected from the model equation, using the modified design matrix (with the inverse of the contrast matrix) and the least squares method, we evaluate the desired constrasts. Of course, to get the original contrast matrix, we need to invert the contrast coding matrix used in R.
Let's try and make it work on an example in R:
x <- rnorm(20,7,2) + 7
y <- rnorm(20,7,2)
z <- rnorm(20,7,2) + 15
t <- rnorm(20,7,2) + 10
df <- data.frame(Score=c(x,y,z,t), Group = c(rep("A",20),rep("B",20),rep("C",20),rep("D",20)))
df$Group <- as.factor(df$Group)
head(df)
Score Group
1 12.83886 A
2 11.49714 A
3 16.27147 A
4 11.84989 A
5 16.00455 A
6 13.78611 A
We have four teams A, B, C, D and the scores of each individual. Let's make the design matrix X for the cell means model:
D <- diag(length(levels(df$Group)))
idx <- as.numeric(df$Group)
X <- matrix(nrow= length(idx), ncol= length(unique(idx)))
for (i in 1:nrow(X)) {
j <- idx[i]
X[i, ] <- D[j,]
}
colnames(X) <- c("A", "B", "C", "D")
head(X)
A B C D
[1,] 1 0 0 0
[2,] 1 0 0 0
[3,] 1 0 0 0
[4,] 1 0 0 0
[5,] 1 0 0 0
[6,] 1 0 0 0
We can find the means of each group by the least squares equation \begin{equation} \hat{\mu} =(X^{\prime }X)^{-1}\ X^{\prime }y\\ \end{equation} in R:
solve( t(X) %*% X) %*% t(X) %*% df$Score
[,1]
A 14.189628
B 7.021692
C 21.668745
D 17.595326
with(df, tapply(X= Score, FUN = mean, INDEX = Group))
A B C D
14.189628 7.021692 21.668745 17.595326
But we want comparisons of means to the first group (treatment contrasts). We use the matrix C of contrasts defined earlier.
Based on what was said before, what we really want is the inverse of C, to evaluate the contrasts.
R has a built-in function for this, called contr.treament(), where we specificy the number of factors.
We build the inverse of C, the contrast coding matrix, this way:
cbind(1, contr.treatment(4) )
2 3 4
1 1 0 0 0
2 1 1 0 0
3 1 0 1 0
4 1 0 0 1
if we invert this matrix, we get C, the comparisons we want:
solve(cbind(1, contr.treatment(4)))
1 0 0 0
-1 1 0 0
-1 0 1 0
-1 0 0 1
Now we construct the modified design matrix for the model:
X1 <- X %*% cbind(1, contr.treatment(4) )
colnames(X1) <- unique(levels(df$Group))
And we solve for the contrasts, either by plugging the modified design matrix into the least squares equation, or using the lm() function:
solve(t(X1) %*% X1) %*% t(X1) %*% df$Score
[,1]
A 14.189628
B -7.167936
C 7.479117
D 3.405698
summary( lm(formula = Score ~ 0 + X1 , data = df) )
Call:
lm(formula = Score ~ 0 + X1, data = df)
Residuals:
Min 1Q Median 3Q Max
-3.5834 -1.2433 -0.1077 1.3763 4.5317
Coefficients:
Estimate Std. Error t value Pr(>|t|)
X1A 14.1896 0.3851 36.845 < 2e-16 ***
X1B -7.1679 0.5446 -13.161 < 2e-16 ***
X1C 7.4791 0.5446 13.732 < 2e-16 ***
X1D 3.4057 0.5446 6.253 2.16e-08 ***
summary( lm(formula = Score ~ Group , data = df, contrasts = list(Group = "contr.treatment")) )
Call:
lm(formula = Score ~ Group, data = df, contrasts = list(Group = "contr.treatment"))
Residuals:
Min 1Q Median 3Q Max
-3.5834 -1.2433 -0.1077 1.3763 4.5317
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 14.1896 0.3851 36.845 < 2e-16 ***
GroupB -7.1679 0.5446 -13.161 < 2e-16 ***
GroupC 7.4791 0.5446 13.732 < 2e-16 ***
GroupD 3.4057 0.5446 6.253 2.16e-08 ***
We get the mean of the first group and the deviations for the others, as defined in the contrast matrix C.
We can define any type of contrast in this way, either using the built-in functions contr.treatment(), contr.sum() etc or by specifying which comparisons we want. There are many refinements on this scheme (orthogonal contrasts, more complex contrasts etc), but this is the gist of it.