Handout on Contrasts Stat 705, 4/5/08
--------------------
## We start with a factor consisting of
## levels A,..,E
> facE = factor(rep(c("A","B","C","D","E"),3))
> contrasts(facE)
B C D E ### These are the default
A 0 0 0 0 ### "treatment" contrasts
B 1 0 0 0
C 0 1 0 0
D 0 0 1 0
E 0 0 0 1
### One R function to set contrasts by name is "C"
> facE2 = C(facE, contr=helmert)
> cmat = contrasts(facE2)
dimnames(cmat)[[2]] = c("B","C","D","E")
cmat
B C D E
A -1 -1 -1 -1
B 1 -1 -1 -1
C 0 2 -1 -1
D 0 0 3 -1
E 0 0 0 4
### Now cmat is a 5 by 4 matrix, which could have
### been defined directly, and then the contrasts
### of the factor could directly have been defined
### equal to that matrix
> facE3 = facE
contrasts(facE3) = contr.helmert(5)
> facE3
[1] A B C D E A B C D E A B C D E
attr(,"contrasts")
B C D E
A -1 -1 -1 -1
B 1 -1 -1 -1
C 0 2 -1 -1
D 0 0 3 -1
E 0 0 0 4
Levels: A B C D E
### Let us have a precise understanding of what
### these "contrasts" matrices are. Suppose we have
### a response variable Y in a vector (in this case
### of length 15, same as facE etc.)
> Yv = runif(15)
### If we do a linear model fit of Y versus one of
### these factors (with intercept), the first five
### rows of the model matrix will be defined equal
### to the corresponding contrast matrix !
> model.matrix(lm(Yv ~ facE))[1:5,]
(Intercept) facEB facEC facED facEE
1 1 0 0 0 0
2 1 1 0 0 0
3 1 0 1 0 0
4 1 0 0 1 0
5 1 0 0 0 1
### If there were no intercept, then the first
### column is replace by the coded dummy indicator
### for level A:
> model.matrix(lm(Yv ~ facE - 1))[1:5,]
facEA facEB facEC facED facEE
1 1 0 0 0 0
2 0 1 0 0 0
3 0 0 1 0 0
4 0 0 0 1 0
5 0 0 0 0 1
### Similarly for the factor with recoded "Helmert" contrasts
> model.matrix(lm(Yv ~ facE3))[1:5,]
(Intercept) facE3B facE3C facE3D facE3E
1 1 -1 -1 -1 -1
2 1 1 -1 -1 -1
3 1 0 2 -1 -1
4 1 0 0 3 -1
5 1 0 0 0 4
### Finally, let's explain exactly how we could have arrived
### at the Helmert contrasts by specifying the interpretation
### we want for the resulting coefficients, as linear
### combinations of group means: "Helmert contrasts" connote
### regression coefficients: beta_1 = mu_A = intercept, and
### beta_2 = mu_B-mu_A, beta_3 = mu_C-mu_B-mu_A,
### beta_4 = mu_D-mu_C-mu_B-mu_A, and
### beta_5 = mu_E-mu_D-mu_C-mu_B-mu_A.
### More generally, let the row-index correspond to the
### indices of a K-component vector Xv with successive
### entries as.character(1:K). Let mu[1:K] be the
### corresponding means for the K successive entries of Yv.
### Now suppose we want to re-code the contrasts, ie make a
### KxK design matrix M for a regression of Yv on Xv, with
### first column of M all 1's, so that : E{ Yv } = mu
### but so that the regression is Yv = M b + mean-0 errors.
### The idea is that we choose b so that b[1]=mu[1] and the
### regression coefficients for its 2:K components have specific
### interpretations as linear combinations of the group means.
### Then M is determined by inverting a matrix.
### In the case of "Helmert contrasts", we want to specify
b[j] = (1/j)*(mu[j]-mean(mu[1:(j-1)]) for j > 1
b[1] = mean(mu)
## so that M has first column 1's and all other columns sum to 0
### after which we obtain M by: M b = mu ,for example when K=5,
> Minv = array(0,c(5,5))
Minv[1,] = rep(.2,5)
for (i in 2:5) Minv[i,] = c(rep(-1,i-1),i-1,rep(0,5-i))/(i*(i-1))
> round(solve(Minv),10)
[,1] [,2] [,3] [,4] [,5]
[1,] 1 -1 -1 -1 -1
[2,] 1 1 -1 -1 -1
[3,] 1 0 2 -1 -1
[4,] 1 0 0 3 -1
[5,] 1 0 0 0 4 ### Exactly the Helmert contrasts.
# ------------ Another example: -----------------------------
### In the case where we want b_j = mu_j-mu_{j-1} for j=2,...,K
### it is easy to check that sum(b[1:j]) = mu[j] so that
### M is the matrix with all 1's on and below the main diagonal
### and all 0's above.