Inverse matrix for contrast coding

Question

I'm trying to understand how "user defined contrast coding" works. My question refers to the example from http://www.ats.ucla.edu/stat/r/library/contrast_coding.htm#User:

#initial contrast matrix including the constant term
mat = matrix(c(1/4, 1/4, 1/4, 1/4, 1, 0, -1, 0, -1/2, 1, 0, -1/2, -1/2, -1/2, 1/2, 1/2), ncol = 4)
mat
     [,1] [,2] [,3] [,4]
[1,] 0.25    1 -0.5 -0.5
[2,] 0.25    0  1.0 -0.5
[3,] 0.25   -1  0.0  0.5
[4,] 0.25    0 -0.5  0.5

mymat = solve(t(mat))
mymat
     [,1] [,2] [,3] [,4]
[1,]    1 -0.5   -1 -1.5
[2,]    1  0.5    1  0.5
[3,]    1 -1.5   -1 -1.5
[4,]    1  1.5    1  2.5

#remove the intercept (constant) term
my.contrasts<-mymat[,2:4]
contrasts(hsb2$race.f) = my.contrasts

Question:

• Why is it necessary to calculate the inverse of the transposed matrix? All other examples on the page doesn't use matrix algebra (e.g "Dummy Coding").

Answer 1

[note: for all those people who got here confused after reading section 9 of the UCLA post on contrast matrices, this will help. after hours of head banging, here's what I figured out]

The above post is correct. You need to actually do an inverse of the transform. In fact, if you look carefully, that is happening for all the other kinds of contrast coding discussed in that post. I will cover the three most common: dummy, simple, and deviation coding.

Please use the following demos to convince yourself of the similarities / difference / relationship between contrasts [as we normally write them for human comprehension] and the thing you see in R [contrast matrix].

First, let look at dummy coding or treatment coding

hsb2 = read.table('http://www.ats.ucla.edu/stat/data/hsb2.csv', header=T, sep=",")

#creating the factor variable race.f
hsb2$race.f = factor(hsb2$race, labels=c("Hispanic", "Asian", "African-Am", "Caucasian"))

#================================================================                       
#                    =: Dummy Coding [normal intercept] := 
#================================================================

# This is just a quick way to create what you would normally consider
# a simple "treatment" coding or "dummy coding" matrix

dummy.coding <- rbind(c(1,-1,-1,-1), 
                      c(0, 1, 0, 0),
                      c(0, 0, 1, 0),
                      c(0, 0, 0, 1))

rownames(dummy.coding)  <- c("control", "L1","L2","L3")
colnames(dummy.coding)  <- c("intercept", "L1 vs Control","L2 vs Control","L3 vs Control")

# Column 2:4 is what you would normally consider the contrast 
# And, as you can see, they add to zero
dummy.coding

# Now, to get the coding matrix that R will use, 
# you need to take the inverse of the transpose of this matrix
# and, drop the first column

coding.matrix <- solve(t(dummy.coding))[,-1]
coding.matrix

# To convince yourself, this is correct, compare it to the 
# inbuilt R function:
contr.treatment(4)

# Now, using the dataset that you cited from the UCLA page
# You can see the what R is doing under the hood:

# First, we assign 
contrasts(hsb2$race.f) = coding.matrix

# Now, we look at the model matrix
View(cbind(as.character(hsb2$race.f), model.matrix(~race.f, data=hsb2)))

# So, essentially, R took that coding matrix and 
# generated dummy variables

# Now, see the linear model:
summary(lm(formula = write ~ race.f, data = temp.hsb2))

OK, having done the basic case, lets extend it. In the above example, the intercept of the model was the "control" variable.

Notice that the final coding matrix did not include a column for the intercept. That is because the inverse of the transpose of the contrast(dummy.coding) baked it into the final outcome: coding.matrix.

Simple Coding

So, now, lets look at a case where we want the intercept to represent the grand mean (i.e. mean of means). This is sometimes called simple coding.

Here, we just set the first column of the contrast matrix to the same value.

#============================================================                        
#=: Dummy Coding [grand mean intercept i.e. simple coding] := 
#============================================================

dummy.coding.GM.intercept <- rbind(c(1,-1,-1,-1), 
c(1, 1, 0, 0),
c(1, 0, 1, 0),
c(1, 0, 0, 1))

rownames(dummy.coding.GM.intercept)  <- c("control", "L1","L2","L3")
colnames(dummy.coding.GM.intercept)  <- c("intercept: GM", "L1 vs Control","L2 vs Control","L3 vs Control")

dummy.coding.GM.intercept

coding.matrix <- solve(t(dummy.coding.GM.intercept))[,-1]
coding.matrix

contrasts(hsb2$race.f) = coding.matrix

# Now, see the linear model and compare the intercepts:
summary(lm(formula = write ~ race.f, data = temp.hsb2))

Deviation Coding or Effects Coding

Finally, lets look at another common case. Here, we want to have contrasts that compare a given level of a factor to the grand mean of the variable. This is sometimes called effects coding, deviation coding, or sum contrasts.

#=======================================================
#     =: Deviation Coding [intercept: grand mean] := 
#=======================================================

deviation.coding <- cbind(c( 0.25,  0.25,  0.25,  0.25), 
      c( 0.75, -0.25, -0.25, -0.25), 
      c(-0.25,  0.75, -0.25, -0.25), 
      c(-0.25, -0.25,  0.75, -0.25))

rownames(deviation.coding)  <- c("L1","L2","L3", "L4")
colnames(deviation.coding)  <- c("intercept: GM", "L1 vs GM","L2 vs GM","L3 vs GM")
deviation.coding

# Now, to get the coding matrix that R will use, 
# you need to take the inverse of the transpose of this matrix
# and, drop the first column

coding.matrix <- solve(t(deviation.coding))[,-1]
coding.matrix

# To convince yourself, this is correct, compare it to the 
# inbuilt R function:
contr.sum(4)

# Now, using the dataset that you cited from the UCLA page
# You can see the what R is doing under the hood:

# First, we assign 
contrasts(hsb2$race.f) = coding.matrix

# Now, we look at the model matrix
View(cbind(as.character(hsb2$race.f), model.matrix(~race.f, data=hsb2)))

# So, essentially, R took that coding matrix and 
# generated variables with EFFECTS CODING 

So, there's that. Hope this helps!!

Answer 2

Ok, i think i have an answer. Please note that I'm not a mathematician, so this isn't a proof.

Let's assume, we have collected a dataset and we want to compare groups.

Our Groups: a, b, c

The comparisons we want to make

• First contrast: a vs b
• Second contrast: b vs c

As an equation this should look like:

\begin{equation*} c_1 = a - b \\ c_2 = b - c \end{equation*}

As a Matrix: \begin{equation*} \begin{pmatrix} c_1 \\ c_2 \end{pmatrix} = \begin{pmatrix} 1 & -1 & 0\\ 0 & 1 & -1 \end{pmatrix} \cdot \begin{pmatrix} a \\ b \\ c \end{pmatrix} \end{equation*}

Or:

\begin{equation*} c = A \cdot v \end{equation*}

But this aren't the droids we were looking for ;-) What we really want is an equation like this:

\begin{equation*} v = X \cdot c \end{equation*}

Or:

\begin{equation*} \begin{pmatrix} a \\ b \\ c \end{pmatrix} = \begin{pmatrix} ? & ? \\ ? & ? \\ ? & ? \end{pmatrix} \cdot \begin{pmatrix} c_1 \\ c_2 \end{pmatrix} \end{equation*} We can achieve this by multiplying by the inverse Matrix (A'):

\begin{equation*} c = A \cdot v \\ A^{-1} \cdot c = A^{-1} \cdot A \cdot v \\ A^{-1} \cdot c = v \\ v = A^{-1} \cdot c \end{equation*}

So, for our example this would be:

\begin{equation*} \begin{pmatrix} a \\ b \\ c \end{pmatrix} = \begin{pmatrix} 2/3 & 1/3 \\ -1/3 & 1/3 \\ -1/3 & -2/3 \end{pmatrix} \cdot \begin{pmatrix} c_1 \\ c_2 \end{pmatrix} \end{equation*}

Please note, that for taking the inverse of a matrix you need a square matrix! If you haven't a square matrix you (or your pc) has to solve the problem numerically. In R this can be done with the function ginv() from the package MASS

Answer 3

This is an old question, and the previous answers were very good, but I will try to answer it, to get a clearer picture. Maybe that can help someone.

Let's get the data and the contrast matrix:

hsb2 = read.table('https://stats.idre.ucla.edu/stat/data/hsb2.csv', header=T, sep=",")
hsb2$race.f = factor(hsb2$race, labels=c("Hispanic", "Asian", "African-Am", "Caucasian"))

mat = matrix(c(1/4, 1/4, 1/4, 1/4, 1, 0, -1, 0, -1/2, 1, 0, -1/2, -1/2, -1/2, 1/2, 1/2), ncol = 4)

mat
     [,1] [,2] [,3] [,4]
[1,] 0.25    1 -0.5 -0.5
[2,] 0.25    0  1.0 -0.5
[3,] 0.25   -1  0.0  0.5
[4,] 0.25    0 -0.5  0.5

The actual contrasts that you want is:

C <- t(mat)
C

      [,1]  [,2]  [,3]  [,4]
[1,]  0.25  0.25  0.25  0.25
[2,]  1.00  0.00 -1.00  0.00
[3,] -0.50  1.00  0.00 -0.50
[4,] -0.50 -0.50  0.50  0.50

And the contrasts are:

\begin{equation} C\mu = \begin{pmatrix} \phantom{..} 1/4 & 1/4 & 1/4 & 1/4 \\ 1 & 0 & -1 & 0\\ -1/2 & 1 & 0 & -1/2\\ -1/2 & -1/2 & 1/2 & 1/2\\ \end{pmatrix}\ \begin{pmatrix}\mu1 \\\mu2 \\\mu3 \\\mu4 \end{pmatrix} \end{equation}

If we calculate them on the data:

means <- with(hsb2, tapply(X = write, INDEX = race.f, FUN = mean))
C %*% means
          [,1]
[1,] 51.678376
[2,] -1.741667
[3,]  7.743247
[4,] -1.101580

Which are the same values given on the website, using their notation:

mymat = solve(t(mat))
summary(lm(write ~ race.f, hsb2, contrasts = list(race.f= mymat[,2:4])))
Call:
lm(formula = write ~ race.f, data = hsb2, contrasts = list(race.f = mymat[, 
    2:4]))

Residuals:
     Min       1Q   Median       3Q      Max 
-23.0552  -5.4583   0.9724   7.0000  18.8000 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  51.6784     0.9821  52.619  < 2e-16 ***
race.f1      -1.7417     2.7325  -0.637  0.52461    
race.f2       7.7432     2.8972   2.673  0.00816 ** 
race.f3      -1.1016     1.9642  -0.561  0.57556 

So we know why we take the transpose, why do we take the inverse?

The model we can use is: \begin{equation} y \ =\ X\mu \ +\ \epsilon \ \end{equation}

with $X$ the design matrix for the cell means model and $\mu$ the vector of means. We can evaluate the means based on this model, using the least squares method:

\begin{equation} \hat{\mu} =(X^{\prime }X)^{-1}\ X^{\prime }y\\ \end{equation}

Now, we constructed the contrast matrix so that C is square and full rank, and we can take its inverse $C^{-1}$, and insert them in the model equation:

\begin{equation} y \ =\ X\mu \ +\ \epsilon \ = \ X I\mu \ +\ \epsilon \ =\ X \ (C^{-1}C)\ \ \mu \ +\ \epsilon = \ (X C^{-1}) \ (C \mu) \ + \epsilon \end{equation}

That means that we can take $XC^{-1}$ as the modified design matrix, to evaluate the contrasts $C\mu$ using the least squares method. In this case, if we name the modified design matrix $X_{1} = XC^{-1}$:

\begin{equation} \hat{C\mu} = (X_{1}^{'}X_{1})^{-1}X_{1}^{'}y \end{equation}

So we use the inverse of the contrast matrix to evaluate the actual contrasts.

To see that it works in R:

X <- model.matrix( ~ hsb2$race.f + 0)  # model matrix for the cell means model
X1 <- X %*% solve(C)                 # this is the modified model matrix, using the inverse.

And we solve by the method of least squares or by lm().

# least squares equations:

solve ( t(X1)  %*% X1 ) %*% t(X1) %*% hsb2$write

          [,1]
[1,] 51.678376
[2,] -1.741667
[3,]  7.743247
[4,] -1.101580

# lm() with modified design matrix:

summary(lm(write ~ X1 + 0, data= hsb2))

Call:
lm(formula = write ~ X1 + 0, data = hsb2)

Residuals:
     Min       1Q   Median       3Q      Max 
-23.0552  -5.4583   0.9724   7.0000  18.8000 

Coefficients:
    Estimate Std. Error t value Pr(>|t|)    
X11  51.6784     0.9821  52.619  < 2e-16 ***
X12  -1.7417     2.7325  -0.637  0.52461    
X13   7.7432     2.8972   2.673  0.00816 ** 
X14  -1.1016     1.9642  -0.561  0.57556 

# lm() with contrast argument:
summary(lm(write ~race.f, data= hsb2, contrasts = list(race.f= MASS::ginv(C[-1,]))))

Call:
lm(formula = write ~ race.f, data = hsb2, contrasts = list(race.f = MASS::ginv(C[-1, 
    ])))

Residuals:
     Min       1Q   Median       3Q      Max 
-23.0552  -5.4583   0.9724   7.0000  18.8000 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  51.6784     0.9821  52.619  < 2e-16 ***
race.f1      -1.7417     2.7325  -0.637  0.52461    
race.f2       7.7432     2.8972   2.673  0.00816 ** 
race.f3      -1.1016     1.9642  -0.561  0.57556    

Which is the same result as above.

Here we generated C so that it is invertible. In the last call, we provided the pseudoinverse of the contrast matrix without the intercept row. lm() adds the intercept column to generate $C^{-1}$. We could have used contrasts = list(race.f = solve(C)[, -1]) .

For pre-defined contrasts, it's the same thing. They are provided without the intercept term in lm(), but it is added internally and the contrasts are evaluated using the least squares method using the modified design matrix. For example, if we use

 contr.treatment(4)
  2 3 4
1 0 0 0
2 1 0 0
3 0 1 0
4 0 0 1

The actual contrasts that are evaluated are:

C <-  solve(cbind(1, contr.treatment(4))) 
    1 0 0 0
   -1 1 0 0
   -1 0 1 0
   -1 0 0 1

And we can evaluate the contrasts as before:

X1 <- X %*% solve(C) 
# least squares equations:

solve ( t(X1)  %*% X1 ) %*% t(X1) %*% hsb2$write

       [,1]
  46.458333
2 11.541667
3  1.741667
4  7.596839

# lm() with modified design matrix:

summary(lm(write ~ X1 +0, data= hsb2))

Call:
lm(formula = write ~ X1 + 0, data = hsb2)

Residuals:
     Min       1Q   Median       3Q      Max 
-23.0552  -5.4583   0.9724   7.0000  18.8000 

Coefficients:
    Estimate Std. Error t value Pr(>|t|)    
X1    46.458      1.842  25.218  < 2e-16 ***
X12   11.542      3.286   3.512 0.000552 ***
X13    1.742      2.732   0.637 0.524613    
X14    7.597      1.989   3.820 0.000179 ***


# lm with pre-defined contrasts:
summary(lm(write ~race.f, data= hsb2, contrasts = list(race.f = contr.treatment(4))))

Call:
lm(formula = write ~ race.f, data = hsb2, contrasts = list(race.f = contr.treatment(4)))

Residuals:
     Min       1Q   Median       3Q      Max 
-23.0552  -5.4583   0.9724   7.0000  18.8000 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   46.458      1.842  25.218  < 2e-16 ***
race.f2       11.542      3.286   3.512 0.000552 ***
race.f3        1.742      2.732   0.637 0.524613    
race.f4        7.597      1.989   3.820 0.000179 ***

And indeed:

C %*% means

       [,1]
  46.458333
2 11.541667
3  1.741667
4  7.596839

$\textbf{Edit in response to @skan}$: When there are more than one categorical variable, the easiest way to proceed is to insert the contrast arguments in the formula, for each categorical variable.

An example is provided here: https://rstudio-pubs-static.s3.amazonaws.com/84177_4604ecc1bae246c9926865db53b6cc29.html in the section "Two factors: one treatment-coded, one deviation-coded".

f <- expand.grid(w=c("left","right"),h=c("low","mid","high")) # 
factor level combinations
d <- data.frame(w=rep(f$w,100),h=rep(f$h,100)) # data frame with 
factor levels only
d$v <-rnorm(nrow(d),100,15)+ifelse(d$w=="left",-3,3)+ifelse(d$h=="low",-5,ifelse(d$h=="high",5,0)) 
str(d)
'data.frame':   600 obs. of  3 variables:
$ w: Factor w/ 2 levels "left","right": 1 2 1 2 1 2 1 2 1 2 ...
$ h: Factor w/ 3 levels "low","mid","high": 1 1 2 2 3 3 1 1 2 2 ...
$ v: num  79.7 121.5 99.5 97.4 117.9 ...

If we use treatment contrast (default) for variable h and sum contrast for variable w, we can enter:

mod <- lm(v ~ w * h, d, contrasts=list(w=contr.sum, h=contr.treatment))
summary(mod)

Call:
lm(formula = v ~ w * h, data = d, contrasts = list(w = contr.sum, 
h = contr.treatment))

Residuals:
Min      1Q  Median      3Q     Max 
-56.194  -9.112   0.287   9.269  47.108 

Coefficients:
        Estimate Std. Error t value Pr(>|t|)    
(Intercept)  93.9555     1.0303  91.191  < 2e-16 ***
w1           -3.9217     1.0303  -3.806 0.000156 ***
h2            5.6757     1.4571   3.895 0.000109 ***
h3           11.8542     1.4571   8.136 2.41e-15 ***
w1:h2         1.0408     1.4571   0.714 0.475312    
w1:h3        -0.6258     1.4571  -0.429 0.667716    

For lme4, this is the same syntax, see for example: How to set custom contrasts with lmer in R

You can use user-defined contrasts, by building each contrast and inverting it, removing the intercept term if needed, and include it in the contrast arg (as described above).

If you want to build your contrasts to run the model without using the contrast arguments, it is more cumbersome (you won't also easily recover the sum of squares for the Anova, using the method below). Here is how to do it for the example above (sum contrast for variable w):

contw <- solve(cbind(1, contr.sum(2)))               # contrast for variable w
conth <- solve(cbind(1, contr.treatment(3)))         # contrast for variable h

contw     
0.5  0.5
0.5 -0.5

conth
 1 0 0
-1 1 0
-1 0 1


f <- expand.grid(w= levels(d$w), h= levels(d$h))     # combination of levels
f
      w    h
1  left  low
2 right  low
3  left  mid
4 right  mid
5  left high
6 right high

# Now we build the contrast matrix C
C <- matrix(nrow= nrow(f), ncol= nrow(f))            
idx <- 0
for (i in 1:nrow(contw)) {
   for (j in 1:nrow(conth)) {
    idx= idx+1
    C[idx,] <- c(contw[i,] %o% conth[j,])  # outer product of contrasts
    }  
}
colnames(C) <- with(f, paste0(w, "_", h))

C
     left_low right_low left_mid right_mid left_high right_high
[1,]      0.5       0.5      0.0       0.0       0.0        0.0
[2,]     -0.5      -0.5      0.5       0.5       0.0        0.0
[3,]     -0.5      -0.5      0.0       0.0       0.5        0.5
[4,]      0.5      -0.5      0.0       0.0       0.0        0.0
[5,]     -0.5       0.5      0.5      -0.5       0.0        0.0
[6,]     -0.5       0.5      0.0       0.0       0.5       -0.5

This is the contrast matrix, which you need to interpret the final results. For example contrast 6 correspond to an interaction term:

1/2 [ (left_high - left_low) - (right_high - right_low ) ]

Now we insert the contrast matrix in the model, using its inverse, as described above:

cond <- with(d, paste0(w, "_", h))
X <- model.matrix(~cond + 0)                     # design matrix
colnames(X) <-  sub("cond", "", colnames(X))
X <- X[, colnames(C)]                            # to make sure the column order is the same as for the contrast matrix
X1 <- X %*% solve(C)                             # modified design matrix

summary(lm(v~X1+0, data= d))                     # solve using the modified design matrix

Call:
lm(formula = v ~ X1 + 0, data = d)

Residuals:
    Min      1Q  Median      3Q     Max 
-56.194  -9.112   0.287   9.269  47.108 

Coefficients:
    Estimate Std. Error t value Pr(>|t|)    
X11  93.9555     1.0303  91.191  < 2e-16 ***
X12   5.6757     1.4571   3.895 0.000109 ***
X13  11.8542     1.4571   8.136 2.41e-15 ***
X14  -3.9217     1.0303  -3.806 0.000156 ***
X15   1.0408     1.4571   0.714 0.475312    
X16  -0.6258     1.4571  -0.429 0.667716    

which corresponds to the results above (using the contrast argument). The order of contrasts is slightly different but values are the same, of course.