Solved – Interpreting non-orthogonal contrasts R

anovacategorical-encodingcontrastsrregression

I'm having trouble finding information on helping me understand how to interpret non-orthogonal contrasts in a regression model.

I have simulated some data to show what I am working with:

set.seed(5552)
Subject <- 1:30
Group <- c(rep("Control", times = 10), rep("Experimental 1", times = 10), rep("Experimental 2", times = 10))
Pre <- round(rnorm(n = Subject, mean = 25, sd = 5), 1)
Post <- c(round(rnorm(n = 10, mean = 25, sd = 5), 1), 
            round(rnorm(n = 10, mean = 28, sd = 3),1),
            round(rnorm(n = 10, mean = 33, sd = 4),1))

dat <- data.frame(Group, Subject, Pre, Post)
dat$Diff <- with(dat, Post - Pre)
dat

Model 1: Here is the regression model using the default contrasts in R:

model1.lm <- lm(Diff ~ Group, data = dat)
summary(model1.lm)

Call:
lm(formula = Diff ~ Group, data = dat)

Residuals:
    Min      1Q  Median      3Q     Max 
-17.250  -2.958  -0.015   3.743  13.950 

Coefficients:
                    Estimate Std. Error t value Pr(>|t|)   
(Intercept)            0.750      1.960   0.383  0.70495   
GroupExperimental 1    3.530      2.772   1.274  0.21365   
GroupExperimental 2    8.430      2.772   3.042  0.00519 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 6.197 on 27 degrees of freedom
Multiple R-squared:  0.2569,    Adjusted R-squared:  0.2018 
F-statistic: 4.666 on 2 and 27 DF,  p-value: 0.01817

The interpretation here is easy as each coefficient is the difference between that group and the reference (control) group, reflected as the intercept.

Model 2: Now, I change the contrasts to orthogonal:

contrast1 <- c(-2,1,1)
contrast2 <- c(0,1,-1)
contrasts(dat$Group) <- cbind(contrast1, contrast2)

Now re-run the regression model:

model1.lm2 <- lm(Diff ~ Group, data = dat)
summary(model1.lm2)

Call:
lm(formula = Diff ~ Group, data = dat)

Residuals:
    Min      1Q  Median      3Q     Max 
-17.250  -2.958  -0.015   3.743  13.950 

Coefficients:
               Estimate Std. Error t value Pr(>|t|)    
(Intercept)      4.7367     1.1315   4.186  0.00027 ***
Groupcontrast1   1.9933     0.8001   2.491  0.01917 *  
Groupcontrast2  -2.4500     1.3858  -1.768  0.08837 .  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 6.197 on 27 degrees of freedom
Multiple R-squared:  0.2569,    Adjusted R-squared:  0.2018 
F-statistic: 4.666 on 2 and 27 DF,  p-value: 0.01817

Okay, this output makes sense. The intercept now reflects the grand mean. The coefficient for contrast 1 (1.99) is the difference between the pooled mean of the 2 experimental groups and the control group, divided by3. The coefficient for contrast 2 now represents the difference between experimental group 2 and experimental group 3, divided by 2.

Model 3: Now, what if I make the contrasts non-orthogonal. I'll keep the first contrast the same (Control Group vs Both Experimental Groups) and change the second contrast to the Control Group vs Experimental Group 2. Because the Control group is entered back into another contrast, this analysis is non-orthogonal. Here are the contrasts:

contrast1 <- c(-2,1,1)
contrast2 <- c(-1,0,1)
contrasts(dat$Group) <- cbind(contrast1, contrast2)

Now re-run the regression:

model1.lm3 <- lm(Diff ~ Group, data = dat)
summary(model1.lm3)

Call:
lm(formula = Diff ~ Group, data = dat)

Residuals:
    Min      1Q  Median      3Q     Max 
-17.250  -2.958  -0.015   3.743  13.950 

Coefficients:
               Estimate Std. Error t value Pr(>|t|)    
(Intercept)      4.7367     1.1315   4.186  0.00027 ***
Groupcontrast1  -0.4567     1.6002  -0.285  0.77753    
Groupcontrast2   4.9000     2.7716   1.768  0.08837 .  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 6.197 on 27 degrees of freedom
Multiple R-squared:  0.2569,    Adjusted R-squared:  0.2018 
F-statistic: 4.666 on 2 and 27 DF,  p-value: 0.01817

How do I best interpret these non-orthogonal contrasts? The intercept is the same as the orthogonal model, the grand mean. Coefficient 1 has changed yet the contrast is the same as the orthogonal model, it doesn't really make sense to me as to why? I can't make sense of the 2nd coefficient either.

Best Answer

The $\beta$s can be interpreted as weighted group means. However, the weights are not given by the coding matrix $\mathbf{B}$ that you specified, but by its inverse, namely, the contrast matrix $\mathbf{C} = \mathbf{B}^{-1}$: $$\beta = \mathbf{C} \mu = \mathbf{B}^{-1} \mu.$$ You can find details in the vignette of the R package codingMatrices, and I have a blog post about it here.

I'll go through your three models to illustrate the interpretation of the $\beta$s. (Instead of calculating the inverse by hand, I'll use codingMatrices::mean_contrasts(), which is really just a wrapper around solve(cbind(1, B)), but has nicer output.)

Dummy Coding

Looking at $\mathbf{C}$ instead of $\mathbf{B}$ we directly see that $\beta_0=1\mu_1$ (m1), $\beta_1=-1\mu_1+1\mu_2$ (-m1+m2), and $\beta_2=-1\mu_1+1\mu_3$.

(B1 <- contr.treatment(3))
#>   2 3
#> 1 0 0
#> 2 1 0
#> 3 0 1
codingMatrices::mean_contrasts(B1)
#>     m1 m2 m3
#> Ave  1  .  .
#> 2   -1  1  .
#> 3   -1  .  1

Orthogonal Contrasts

Here, $\mathbf{C}$ and $\mathbf{B}$ are structurally so similar that you knew the meaning of the $\beta$s without explicitly knowing $\mathbf{C}$.

(B2 <- cbind(c(-2, 1, 1), c(0, 1, -1)))
#>      [,1] [,2]
#> [1,]   -2    0
#> [2,]    1    1
#> [3,]    1   -1
codingMatrices::mean_contrasts(B2)
#>     m1   m2   m3  
#> Ave  1/3  1/3  1/3
#>     -1/3  1/6  1/6
#>        .  1/2 -1/2

Nonorthogonal Contrasts I

Here, $\mathbf{C}$ and $\mathbf{B}$ are structurally completely different. You wanted to compare $\mu_3$ with $\mu_1$ but actually compared it with $\mu_2$ ($9.18 - 4.28 = 4.90$). As you can see, multiplying the contrasts of $\mathbf{C}$ with the group means gives you the $\beta$s.

(B3 <- cbind(c(-2, 1, 1), c(-1, 0, 1)))
#>      [,1] [,2]
#> [1,]   -2   -1
#> [2,]    1    0
#> [3,]    1    1
codingMatrices::mean_contrasts(B3)
#>     m1   m2   m3  
#> Ave  1/3  1/3  1/3
#>     -1/3  2/3 -1/3
#>        .   -1    1

(group_means <- tapply(dat$Diff, dat$Group, mean))
#>        Control Experimental 1 Experimental 2 
#>           0.75           4.28           9.18
codingMatrices::mean_contrasts(B3) %*% array(group_means, dim = c(3, 1))
#>           [,1]
#> Ave  4.7366667
#>     -0.4566667
#>      4.9000000

model1.lm3 <- lm(Diff ~ Group, data = dat, contrasts = list(Group = B3))
coef(model1.lm3)
#> (Intercept)      Group1      Group2 
#>   4.7366667  -0.4566667   4.9000000

Nonorthogonal Contrasts II

We can use this the other way round and specify a matrix $\mathbf{C}$ to to compare $\mu_3$ with $\mu_1$. Then, we use its inverse as input to lm(), and this will give you the desired result, namely the difference between group 3 and group 1 which is $9.18 - 0.75 = 8.43$.

(C4 <- rbind(1/3, c(-1, 0.5, 0.5), c(-1, 0, 1)))
#>            [,1]      [,2]      [,3]
#> [1,]  0.3333333 0.3333333 0.3333333
#> [2,] -1.0000000 0.5000000 0.5000000
#> [3,] -1.0000000 0.0000000 1.0000000
solve(C4)
#>      [,1]       [,2] [,3]
#> [1,]    1 -0.6666667    0
#> [2,]    1  1.3333333   -1
#> [3,]    1 -0.6666667    1

model1.lm4 <- lm(Diff ~ Group, data = dat,
                 contrasts = list(Group = solve(C4)[, -1]))
coef(model1.lm4)
#> (Intercept)      Group1      Group2 
#>    4.736667    5.980000    8.430000

In summary, if you want to know how to interpret the estimated $\beta$s, you have to look at the inverse of your design/coding matrix. It may be obvious for standard schemes (e.g., dummy) or for orthogonal contrasts, but definitely not for non-standard non-orthongonal coding schemes.

Related Question