Learning note of using contrast matrix for sum of squares
Mark Lai
April 12, 2024
Consider 4 groups, each with 10 observations, and the group means are 5, 7, 9, and 11.
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 unbalancedUsing lm
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 ' ' 1Using contrast matrix
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 ' ' 1Notice 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,
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}}], \]
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 unbalancedUsing contrast coding
(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 0Analysis 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 ' ' 1Using contrast matrix
(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 [,1]
[1,] 90 [,1]
[1,] 90 [,1]
[1,] 30\[ [\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.
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 unbalancedUsing lm
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 ' ' 1Using contrast matrix
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---
title: "Variance explained for contrasts"
description: |
Learning note of using contrast matrix for sum of squares
author: Mark Lai
date: "2024-04-12"
categories:
- Statistics
---
## One-Way ANOVA
Consider 4 groups, each with 10 observations, and the group means are 5, 7, 9, and 11.
```{r}
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)
```
Using `lm`
```{r}
lm1 <- lm(y ~ group)
anova(lm1)
```
Using contrast matrix
```{r}
fac <- factor(group)
contrasts(fac) <- contr.sum(4) # effect coding
lm2 <- lm(y ~ fac)
anova(lm2)
```
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,
```{r}
contr.sum(4)
```
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](https://onlinelibrary.wiley.com/doi/book/10.1002/9780470192610))
$$
n [\mathbf{L}^\intercal \bar{\mathbf{Y}}]^\intercal [\mathbf{L}^\intercal \mathbf{L}]^{-1} [\mathbf{L}^\intercal \bar{\mathbf{Y}}],
$$
```{r}
L_mat <- contr.sum(length(mu))
Lmu <- crossprod(L_mat, mu)
num_obs * crossprod(Lmu, solve(crossprod(L_mat), Lmu))
```
## Interaction
```{r}
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)
```
Using contrast coding
```{r}
fac <- factor(group)
contrasts(fac) <- contr.sum(4)
fac2 <- factor(group2)
contrasts(fac2) <- contr.sum(2)
model.matrix(~ fac * fac2) |> head()
lm(y2 ~ fac * fac2) |> anova()
```
Using contrast matrix
```{r}
# Extract the contrast matrix
df2 <- unique(data.frame(fac, fac2))
(X <- model.matrix(~ fac * fac2, data = df2))
# For main effect of group
L_mat <- X[, 2:4]
Lmu <- crossprod(L_mat, mu2)
n2 * crossprod(Lmu, solve(crossprod(L_mat), Lmu))
# For main effect of group2
L_mat2 <- X[, 5, drop = FALSE]
Lmu2 <- crossprod(L_mat2, mu2)
n2 * crossprod(Lmu2, solve(crossprod(L_mat2), Lmu2))
# For interaction
L_mat3 <- X[, 6:8]
Lmu3 <- crossprod(L_mat3, mu2)
n2 * crossprod(Lmu3, solve(crossprod(L_mat3), Lmu3))
```
## 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.
```{r}
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)
```
Using `lm`
```{r}
lm1 <- lm(y ~ group)
anova(lm1)
```
Using contrast matrix
```{r}
fac <- factor(group)
contrasts(fac) <- contr.helmert(4)
lm2 <- lm(y ~ fac)
anova(lm2)
```
```{r}
L_mat <- contr.sum(length(mu))
Lmu <- crossprod(L_mat, mu)
crossprod(Lmu, solve(crossprod(L_mat, L_mat / num_obs), Lmu))
```