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.