I’m trying to understand some of the basics of contrasts in analysis of variance. In particular, I want to figure out some alternative parameterizations of the coefficients that are estimated as part of the ANOVA model.
library(MASS)
library(matlib)## Warning: package 'matlib' was built under R version 4.1.3library(tidyverse)## -- Attaching packages --------------------------------------- tidyverse 1.3.1 --## v ggplot2 3.3.5 v purrr 0.3.4
## v tibble 3.1.7 v dplyr 1.0.9
## v tidyr 1.1.3 v stringr 1.4.0
## v readr 2.0.1 v forcats 0.5.1## Warning: package 'tibble' was built under R version 4.1.3## Warning: package 'dplyr' was built under R version 4.1.3## -- Conflicts ------------------------------------------ tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()
## x dplyr::select() masks MASS::select()knitr::opts_chunk$set(echo = TRUE)First, let’s set up an artificial dataset that is nice and easy to work with.
means <- c(52, 54, 62, 72)
mean(means)## [1] 60y <-
rep(means, each=3) +
rep(c(0, 1, -1), 4)
g <- rep(1:4, each=3)
df <- cbind(g, y)Notice the pattern in the group means.
data.frame(g, y) %>%
group_by(g) %>%
summarize(mean_y=mean(y))## # A tibble: 4 x 2
## g mean_y
## <int> <dbl>
## 1 1 52
## 2 2 54
## 3 3 62
## 4 4 72I want to use the lm function to fit the data, because it is easier to peek at the underlying coefficients.
The simplest parameterization
The simplest parameterization is to let each beta coefficient represent a different mean. We have to drop the intercept term for this model to work properly.
x1 <- as.numeric(g==1)
x2 <- as.numeric(g==2)
x3 <- as.numeric(g==3)
x4 <- as.numeric(g==4)This is what the X matrix looks like.
A <- cbind(x1, x2, x3, x4)
A## x1 x2 x3 x4
## [1,] 1 0 0 0
## [2,] 1 0 0 0
## [3,] 1 0 0 0
## [4,] 0 1 0 0
## [5,] 0 1 0 0
## [6,] 0 1 0 0
## [7,] 0 0 1 0
## [8,] 0 0 1 0
## [9,] 0 0 1 0
## [10,] 0 0 0 1
## [11,] 0 0 0 1
## [12,] 0 0 0 1This is what the regression results look like.
beta_A <- coef(lm(y~x1+x2+x3+x4-1))
beta_A## x1 x2 x3 x4
## 52 54 62 72It may help to view the underlying matrix calculations. Recall that
\(\hat{\beta}=(X'X)^{-1}X'Y\)
So the matrix
\((X'X)^{-1}X'\)
will tell you how each parameter in the regression model is a linear combination of the Y’s.
inv(t(A)%*%A)%*%t(A) %>% fractions## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12]
## [1,] 1/3 1/3 1/3 0 0 0 0 0 0 0 0 0
## [2,] 0 0 0 1/3 1/3 1/3 0 0 0 0 0 0
## [3,] 0 0 0 0 0 0 1/3 1/3 1/3 0 0 0
## [4,] 0 0 0 0 0 0 0 0 0 1/3 1/3 1/3So the first coefficient is an average of the first three y values, the second is an average of the next three y values, etc.
Indicator variables
A common way to parameterize the ANOVA model is to use an intercept term and indicator variables for all but one of the groups.
x0 <- rep(1, 12)
x1 <- as.numeric(g==2)
x2 <- as.numeric(g==3)
x3 <- as.numeric(g==4)This is what the X matrix looks like.
B <- cbind(x0, x1, x2, x3)
B## x0 x1 x2 x3
## [1,] 1 0 0 0
## [2,] 1 0 0 0
## [3,] 1 0 0 0
## [4,] 1 1 0 0
## [5,] 1 1 0 0
## [6,] 1 1 0 0
## [7,] 1 0 1 0
## [8,] 1 0 1 0
## [9,] 1 0 1 0
## [10,] 1 0 0 1
## [11,] 1 0 0 1
## [12,] 1 0 0 1This is what the regression results look like.
beta_B <- coef(lm(y~x0+x1+x2+x3-1))
beta_B## x0 x1 x2 x3
## 52 2 10 20inv(t(B)%*%B)%*%t(B) %>% fractions## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12]
## [1,] 1/3 1/3 1/3 0 0 0 0 0 0 0 0 0
## [2,] -1/3 -1/3 -1/3 1/3 1/3 1/3 0 0 0 0 0 0
## [3,] -1/3 -1/3 -1/3 0 0 0 1/3 1/3 1/3 0 0 0
## [4,] -1/3 -1/3 -1/3 0 0 0 0 0 0 1/3 1/3 1/3So the first coefficient is an average of the first three y values (Y1, Y2, Y3), the second is the difference between the average of the second three y values (Y4, Y5, Y6) and the average of the first three y values (Y1, Y2, Y3), etc.
In other words,
\(\hat{\beta}_0=\bar{Y}_1\)
\(\hat{\beta}_1=\bar{Y}_2-\bar{Y}_1\)
\(\hat{\beta}_2=\bar{Y}_3-\bar{Y}_1\)
\(\hat{\beta}_3=\bar{Y}_4-\bar{Y}_1\)
This seems counter intuitive. The slope coefficient for \(\beta_2\) involves a column with zeros everywhere except for some ones for Y4, Y5, and Y6. How could that coefficient end up involving Y1, Y2, and Y3?
The point to remember is that the slope coefficient in a linear regression with multiple independent variables has a slightly more complex interpretation than the slope coefficient in a linear regression with a single independent variable. In the later case, the interpretation would be
- the estimated average change in \(Y\) when \(X\) increases by one unit.
In the former case, the interpretation would be
- the estimated average change in \(Y\) when \(X_i\) increases by one unit while holding all the other \(X\)’s constant.
So the slope associated with X1 requires you to hold X2 and X3 constant. You can hold them constant at 0 or hold them constant at 1, but the former is simpler to understand. Think of it as moving from one of the first three rows of your X matrix.
B[1:3, ]## x0 x1 x2 x3
## [1,] 1 0 0 0
## [2,] 1 0 0 0
## [3,] 1 0 0 0to one of the next three rows of your X matrix
B[4:6, ]## x0 x1 x2 x3
## [1,] 1 1 0 0
## [2,] 1 1 0 0
## [3,] 1 1 0 0This is effectively a change from the mean of the first three observations (Y1, Y2, Y3) subtracted from the mean of the next three observations (Y4, Y5, Y6).
Deviation coding
If you tried to code the difference in means directly, you get a third result.A common way to parameterize the ANOVA model is to use an intercept term and indicator variables for all but one of the groups.
x0 <- rep(1, 12)
x1 <- as.numeric(g==2) - as.numeric(g==1)
x2 <- as.numeric(g==3) - as.numeric(g==1)
x3 <- as.numeric(g==4) - as.numeric(g==1)This is what the X matrix looks like.
D <- cbind(x0, x1, x2, x3)
D## x0 x1 x2 x3
## [1,] 1 -1 -1 -1
## [2,] 1 -1 -1 -1
## [3,] 1 -1 -1 -1
## [4,] 1 1 0 0
## [5,] 1 1 0 0
## [6,] 1 1 0 0
## [7,] 1 0 1 0
## [8,] 1 0 1 0
## [9,] 1 0 1 0
## [10,] 1 0 0 1
## [11,] 1 0 0 1
## [12,] 1 0 0 1This is what the regression results look like.
beta_D <- coef(lm(y~x0+x1+x2+x3-1))
beta_D## x0 x1 x2 x3
## 60 -6 2 12inv(t(D)%*%D)%*%t(D) %>% fractions## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12]
## [1,] 1/12 1/12 1/12 1/12 1/12 1/12 1/12 1/12 1/12 1/12 1/12 1/12
## [2,] -1/12 -1/12 -1/12 1/4 1/4 1/4 -1/12 -1/12 -1/12 -1/12 -1/12 -1/12
## [3,] -1/12 -1/12 -1/12 -1/12 -1/12 -1/12 1/4 1/4 1/4 -1/12 -1/12 -1/12
## [4,] -1/12 -1/12 -1/12 -1/12 -1/12 -1/12 -1/12 -1/12 -1/12 1/4 1/4 1/4Okay, so what is this 1/4 doing here? There are only three data values in each group. You need to rewrite
1/4 = 1/3 - 1/12
Then the -1/12 represents subtracting the overall mean of all 12 observations. So the estimates involve an average of a group minus the overall average. In other words,
\(\hat{\beta}_0=\bar{Y}_.\)
\(\hat{\beta}_1=\bar{Y}_2-\bar{Y}_.\)
\(\hat{\beta}_2=\bar{Y}_3-\bar{Y}_.\)
\(\hat{\beta}_3=\bar{Y}_4-\bar{Y}_.\)
where \(\bar{Y}_.\) is the overall mean.
Again, this seems counter-intuitive, but remember that the interpretation of the regression coefficient in a multivariate setting is the effect of a one unit change in a variable holding all the other variables constant.
Recode
Suppose you want to recombine the columns of your X matrix to transform it using some linear combinations. This is the same as post-multiplying by a matrix.
\(B=AV\)
or equivalently
\(A=BV^{-1}\)
So to convert from
A## x1 x2 x3 x4
## [1,] 1 0 0 0
## [2,] 1 0 0 0
## [3,] 1 0 0 0
## [4,] 0 1 0 0
## [5,] 0 1 0 0
## [6,] 0 1 0 0
## [7,] 0 0 1 0
## [8,] 0 0 1 0
## [9,] 0 0 1 0
## [10,] 0 0 0 1
## [11,] 0 0 0 1
## [12,] 0 0 0 1into
B## x0 x1 x2 x3
## [1,] 1 0 0 0
## [2,] 1 0 0 0
## [3,] 1 0 0 0
## [4,] 1 1 0 0
## [5,] 1 1 0 0
## [6,] 1 1 0 0
## [7,] 1 0 1 0
## [8,] 1 0 1 0
## [9,] 1 0 1 0
## [10,] 1 0 0 1
## [11,] 1 0 0 1
## [12,] 1 0 0 1you would multiply by this V
## [,1] [,2] [,3] [,4]
## [1,] 1 0 0 0
## [2,] 1 1 0 0
## [3,] 1 0 1 0
## [4,] 1 0 0 1to get
A%*%V## [,1] [,2] [,3] [,4]
## [1,] 1 0 0 0
## [2,] 1 0 0 0
## [3,] 1 0 0 0
## [4,] 1 1 0 0
## [5,] 1 1 0 0
## [6,] 1 1 0 0
## [7,] 1 0 1 0
## [8,] 1 0 1 0
## [9,] 1 0 1 0
## [10,] 1 0 0 1
## [11,] 1 0 0 1
## [12,] 1 0 0 1Note the V is not a symmetric matrix. Symmetry is not requirement here.
You can go in the other direction using
U=inv(V)
U## [,1] [,2] [,3] [,4]
## [1,] 1 0 0 0
## [2,] -1 1 0 0
## [3,] -1 0 1 0
## [4,] -1 0 0 1which yields
B%*%U## [,1] [,2] [,3] [,4]
## [1,] 1 0 0 0
## [2,] 1 0 0 0
## [3,] 1 0 0 0
## [4,] 0 1 0 0
## [5,] 0 1 0 0
## [6,] 0 1 0 0
## [7,] 0 0 1 0
## [8,] 0 0 1 0
## [9,] 0 0 1 0
## [10,] 0 0 0 1
## [11,] 0 0 0 1
## [12,] 0 0 0 1Now the formula for the estimates for a given parameterization is
\(\hat{\beta}_B=(B'B)^{-1}B'Y\)
Replace B with AV to get
\(\hat{\beta}_B=(V'A'AV)^{-1}V'A'Y\)
which simplifies, after a bit of work, to
\(\hat{\beta}_B=V^{-1}(A'A)^{-1}A'Y\)
or
\(\hat{\beta}_B=V^{-1}\hat{\beta}_A\)
Let’s check this out with the coefficients we calculated earlier.
beta_A## x1 x2 x3 x4
## 52 54 62 72beta_B## x0 x1 x2 x3
## 52 2 10 20inv(V)%*%beta_A## [,1]
## [1,] 52
## [2,] 2
## [3,] 10
## [4,] 20You can also go in the reverse direction.
V%*%beta_B## [,1]
## [1,] 52
## [2,] 54
## [3,] 62
## [4,] 72Check the deviation coding
The same works with any coding. The V matrix that gets you from
A## x1 x2 x3 x4
## [1,] 1 0 0 0
## [2,] 1 0 0 0
## [3,] 1 0 0 0
## [4,] 0 1 0 0
## [5,] 0 1 0 0
## [6,] 0 1 0 0
## [7,] 0 0 1 0
## [8,] 0 0 1 0
## [9,] 0 0 1 0
## [10,] 0 0 0 1
## [11,] 0 0 0 1
## [12,] 0 0 0 1to
D## x0 x1 x2 x3
## [1,] 1 -1 -1 -1
## [2,] 1 -1 -1 -1
## [3,] 1 -1 -1 -1
## [4,] 1 1 0 0
## [5,] 1 1 0 0
## [6,] 1 1 0 0
## [7,] 1 0 1 0
## [8,] 1 0 1 0
## [9,] 1 0 1 0
## [10,] 1 0 0 1
## [11,] 1 0 0 1
## [12,] 1 0 0 1is
## [,1] [,2] [,3] [,4]
## [1,] 1 -1 -1 -1
## [2,] 1 1 0 0
## [3,] 1 0 1 0
## [4,] 1 0 0 1Let’s check this in one direction
A%*%V## [,1] [,2] [,3] [,4]
## [1,] 1 -1 -1 -1
## [2,] 1 -1 -1 -1
## [3,] 1 -1 -1 -1
## [4,] 1 1 0 0
## [5,] 1 1 0 0
## [6,] 1 1 0 0
## [7,] 1 0 1 0
## [8,] 1 0 1 0
## [9,] 1 0 1 0
## [10,] 1 0 0 1
## [11,] 1 0 0 1
## [12,] 1 0 0 1and the opposite direction
U=inv(V)
U## [,1] [,2] [,3] [,4]
## [1,] 0.25 0.25 0.25 0.25
## [2,] -0.25 0.75 -0.25 -0.25
## [3,] -0.25 -0.25 0.75 -0.25
## [4,] -0.25 -0.25 -0.25 0.75D%*%U## [,1] [,2] [,3] [,4]
## [1,] 1 0 0 0
## [2,] 1 0 0 0
## [3,] 1 0 0 0
## [4,] 0 1 0 0
## [5,] 0 1 0 0
## [6,] 0 1 0 0
## [7,] 0 0 1 0
## [8,] 0 0 1 0
## [9,] 0 0 1 0
## [10,] 0 0 0 1
## [11,] 0 0 0 1
## [12,] 0 0 0 1We can establish that
beta_A## x1 x2 x3 x4
## 52 54 62 72can be converted to
beta_D## x0 x1 x2 x3
## 60 -6 2 12using
inv(V)%*%beta_A## [,1]
## [1,] 60
## [2,] -6
## [3,] 2
## [4,] 12and the reverse direction works as well.
V%*%beta_D## [,1]
## [1,] 52
## [2,] 54
## [3,] 62
## [4,] 72Make your own contrast
Suppose you wanted to develop a parameterization that is meaningful for your particular setting and none of the ones discussed here or elsewhere seem to fit the bill. With a bit of work, you can calculate your own.
Suppose for example that you wanted
\(\hat{\beta}_0=\bar{Y}_1\)
\(\hat{\beta}_1=\bar{Y}_2-\bar{Y}_1\)
\(\hat{\beta}_2=\bar{Y}_3-\bar{Y}_2\)
\(\hat{\beta}_3=\bar{Y}_4-\bar{Y}_3\)
Instead of comparing each mean to the first mean, you want to compare each mean to the previous mean.
The U matrix that converts
\(\hat{\beta}_0=\bar{Y}_1\)
\(\hat{\beta}_1=\bar{Y}_2\)
\(\hat{\beta}_2=\bar{Y}_3\)
\(\hat{\beta}_3=\bar{Y}_4\)
to the new parameterization is
## [,1] [,2] [,3] [,4]
## [1,] 1 0 0 0
## [2,] -1 1 0 0
## [3,] 0 -1 1 0
## [4,] 0 0 -1 1Now, you want to travel from the parameterization represented by
A## x1 x2 x3 x4
## [1,] 1 0 0 0
## [2,] 1 0 0 0
## [3,] 1 0 0 0
## [4,] 0 1 0 0
## [5,] 0 1 0 0
## [6,] 0 1 0 0
## [7,] 0 0 1 0
## [8,] 0 0 1 0
## [9,] 0 0 1 0
## [10,] 0 0 0 1
## [11,] 0 0 0 1
## [12,] 0 0 0 1to a new parameterization represented by E. This would need V, the inverse of U, which is
V <- inv(U)
V## [,1] [,2] [,3] [,4]
## [1,] 1 0 0 0
## [2,] 1 1 0 0
## [3,] 1 1 1 0
## [4,] 1 1 1 1which provides a matrix
E <- A%*%inv(U)
E## [,1] [,2] [,3] [,4]
## [1,] 1 0 0 0
## [2,] 1 0 0 0
## [3,] 1 0 0 0
## [4,] 1 1 0 0
## [5,] 1 1 0 0
## [6,] 1 1 0 0
## [7,] 1 1 1 0
## [8,] 1 1 1 0
## [9,] 1 1 1 0
## [10,] 1 1 1 1
## [11,] 1 1 1 1
## [12,] 1 1 1 1which gives the correct results.
x0 <- E[ , 1]
x1 <- E[ , 2]
x2 <- E[ , 3]
x3 <- E[ , 4]
beta_E <- coef(lm(y~x0+x1+x2+x3-1))
beta_E## x0 x1 x2 x3
## 52 2 8 10Acknowledgements
In writing this webpage, I relied on a very helpful tutorial paper:
Daniel J. Schad, Shravan Vasishth, Sven Hohenstein, Reinhold Kliegl. How to capitalize on a priori contrasts in linear (mixed) models: A tutorial. Journal of Memory and Language, 2020, 110. Available in html format or pdf format.
as well as a page on the UCLA Statistics site
R Library Contrast Coding Systems for Categorical Variables. Available in html format