Variance explained for contrasts

Statistics

Learning note of using contrast matrix for sum of squares

Author

Mark Lai

Published

April 12, 2024

One-Way ANOVA

Consider 4 groups, each with 10 observations, and the group means are 5, 7, 9, and 11.

set.seed(2201)
num_obs <- 10
err <- rnorm(num_obs)
err <- (err - mean(err))
mu <- c(11, 7, 5, 9)
y <- rep(mu, each = num_obs) + err
group <- rep(LETTERS[1:4], each = num_obs)
aov(y ~ group)

Call:
   aov(formula = y ~ group)

Terms:
                  group Residuals
Sum of Squares  200.000    46.247
Deg. of Freedom       3        36

Residual standard error: 1.133419
Estimated effects may be unbalanced

Using lm

lm1 <- lm(y ~ group)
anova(lm1)

Analysis of Variance Table

Response: y
          Df  Sum Sq Mean Sq F value    Pr(>F)    
group      3 200.000  66.667  51.895 3.743e-13 ***
Residuals 36  46.247   1.285                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Using contrast matrix

fac <- factor(group)
contrasts(fac) <- contr.sum(4)  # effect coding
lm2 <- lm(y ~ fac)
anova(lm2)

Analysis of Variance Table

Response: y
          Df  Sum Sq Mean Sq F value    Pr(>F)    
fac        3 200.000  66.667  51.895 3.743e-13 ***
Residuals 36  46.247   1.285                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Notice that we get the same sum of squares (SS). Next, we’ll compute SS by hand.

Let \(\mu_j\) be the population mean for group \(j\), and \(\bar Y_{.j}\) be the sample group mean. Let \(\mathbf{L}\) be a matrix such that each column is a contrast vector. For example,

contr.sum(4)

  [,1] [,2] [,3]
1    1    0    0
2    0    1    0
3    0    0    1
4   -1   -1   -1

Note that each column sums to zero. Then, for a balanced design, the SS accounted for by \(\mathbf{L}\) is (e.g., Rencher & Schaalje, 2008)

\[ n [\mathbf{L}^\intercal \bar{\mathbf{Y}}]^\intercal [\mathbf{L}^\intercal \mathbf{L}]^{-1} [\mathbf{L}^\intercal \bar{\mathbf{Y}}], \]

L_mat <- contr.sum(length(mu))
Lmu <- crossprod(L_mat, mu)
num_obs * crossprod(Lmu, solve(crossprod(L_mat), Lmu))

     [,1]
[1,]  200

Interaction

n2 <- num_obs / 2
err <- rnorm(n2)
err <- (err - mean(err))
mu2 <- c(5, 7, 9, 11, 6, 12, 8, 10)
group2 <- rep(rep(1:2, each = n2), length(mu))
y2 <- rep(mu2, each = n2) + err
aov(y2 ~ group * group2)

Call:
   aov(formula = y2 ~ group * group2)

Terms:
                   group   group2 group:group2 Residuals
Sum of Squares  90.00000 90.00000     30.00000  21.47982
Deg. of Freedom        3        1            3        32

Residual standard error: 0.8192951
Estimated effects may be unbalanced

Using contrast coding

fac <- factor(group)
contrasts(fac) <- contr.sum(4)
fac2 <- factor(group2)
contrasts(fac2) <- contr.sum(2)
model.matrix(~ fac * fac2) |> head()

  (Intercept) fac1 fac2 fac3 fac21 fac1:fac21 fac2:fac21 fac3:fac21
1           1    1    0    0     1          1          0          0
2           1    1    0    0     1          1          0          0
3           1    1    0    0     1          1          0          0
4           1    1    0    0     1          1          0          0
5           1    1    0    0     1          1          0          0
6           1    1    0    0    -1         -1          0          0

lm(y2 ~ fac * fac2) |> anova()

Analysis of Variance Table

Response: y2
          Df Sum Sq Mean Sq F value    Pr(>F)    
fac        3  90.00  30.000  44.693 1.508e-11 ***
fac2       1  90.00  90.000 134.079 5.582e-13 ***
fac:fac2   3  30.00  10.000  14.898 3.037e-06 ***
Residuals 32  21.48   0.671                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Using contrast matrix

# Extract the contrast matrix
df2 <- unique(data.frame(fac, fac2))
(X <- model.matrix(~ fac * fac2, data = df2))

   (Intercept) fac1 fac2 fac3 fac21 fac1:fac21 fac2:fac21 fac3:fac21
1            1    1    0    0     1          1          0          0
6            1    1    0    0    -1         -1          0          0
11           1    0    1    0     1          0          1          0
16           1    0    1    0    -1          0         -1          0
21           1    0    0    1     1          0          0          1
26           1    0    0    1    -1          0          0         -1
31           1   -1   -1   -1     1         -1         -1         -1
36           1   -1   -1   -1    -1          1          1          1
attr(,"assign")
[1] 0 1 1 1 2 3 3 3
attr(,"contrasts")
attr(,"contrasts")$fac
  [,1] [,2] [,3]
A    1    0    0
B    0    1    0
C    0    0    1
D   -1   -1   -1

attr(,"contrasts")$fac2
  [,1]
1    1
2   -1

# For main effect of group
L_mat <- X[, 2:4]
Lmu <- crossprod(L_mat, mu2)
n2 * crossprod(Lmu, solve(crossprod(L_mat), Lmu))

     [,1]
[1,]   90

# For main effect of group2
L_mat2 <- X[, 5, drop = FALSE]
Lmu2 <- crossprod(L_mat2, mu2)
n2 * crossprod(Lmu2, solve(crossprod(L_mat2), Lmu2))

     [,1]
[1,]   90

# For interaction
L_mat3 <- X[, 6:8]
Lmu3 <- crossprod(L_mat3, mu2)
n2 * crossprod(Lmu3, solve(crossprod(L_mat3), Lmu3))

     [,1]
[1,]   30

Unbalanced data

\[ [\mathbf{L}^\intercal \bar{\mathbf{Y}}]^\intercal [\mathbf{L}^\intercal \mathbf{W}^{-1} \mathbf{L}]^{-1} [\mathbf{L}^\intercal \bar{\mathbf{Y}}], \]

where \(\mathbf{W}\) is the diagonal matrix of group sizes.

set.seed(2203)
num_obs <- c(10, 5, 16, 9)
err <- rnorm(sum(num_obs))
err <- err - ave(err, rep.int(seq_along(num_obs), times = num_obs))
mu <- c(9, 5, 7, 11)
y <- rep(mu, num_obs) + err
group <- rep(LETTERS[1:4], num_obs)
aov(y ~ group)

Call:
   aov(formula = y ~ group)

Terms:
                    group Residuals
Sum of Squares  151.10000  37.66745
Deg. of Freedom         3        36

Residual standard error: 1.022897
Estimated effects may be unbalanced

Using lm

lm1 <- lm(y ~ group)
anova(lm1)

Analysis of Variance Table

Response: y
          Df  Sum Sq Mean Sq F value    Pr(>F)    
group      3 151.100  50.367  48.137 1.107e-12 ***
Residuals 36  37.667   1.046                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Using contrast matrix

fac <- factor(group)
contrasts(fac) <- contr.helmert(4)
lm2 <- lm(y ~ fac)
anova(lm2)

Analysis of Variance Table

Response: y
          Df  Sum Sq Mean Sq F value    Pr(>F)    
fac        3 151.100  50.367  48.137 1.107e-12 ***
Residuals 36  37.667   1.046                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

L_mat <- contr.sum(length(mu))
Lmu <- crossprod(L_mat, mu)
crossprod(Lmu, solve(crossprod(L_mat, L_mat / num_obs), Lmu))

      [,1]
[1,] 151.1