Group Comparison

Question: Are male adelie penguins heavier than female adelie penguins ?

par(mfrow = c(1, 2), mar = c(5, 4, 2, 2) + 0.1)
plot(adelie$body_mass_g, col = adelie$sex)      ## scatter plot
boxplot(body_mass_g ~ sex, data = adelie)       ## box plot

Simple \(t\) test

Data:

• \(n_f\) female adelie penguins, \(n_m\) male adelie penguins
• \(y^f_{i}\): mass (g) of a female adelie penguin, \(1 \leq i \leq n_f\)
• \(y^m_{j}\): mass (g) of a male adelie penguin, \(1 \leq j \leq n_m\)

Assumption:

• \(y^f_{i} \sim \mathcal{N}(\mu_f, \sigma^2)\) iid
• \(y^m_{j} \sim \mathcal{N}(\mu_m, \sigma^2)\) iid

Test:

\(\mathcal{H}_0\): \(\mu_f = \mu_m\) vs \(\mathcal{H}_1\): \(\mu_f \neq \mu_m\)

Simple \(t\) test

Test:

\(\mathcal{H}_0\): \(\mu_f = \mu_m\) vs \(\mathcal{H}_1\): \(\mu_f \neq \mu_m\)

Test Statistic:

\[ T = \frac{\bar{\mathbf{y}}^f - \bar{\mathbf{y}}^m}{s_p \sqrt{\frac{1}{n_f} + \frac{1}{n_m}}} \underset{\mathcal{H}_0}{\sim} \mathcal{T}_{n_f + n_m - 2} \] with the “pooled” standard deviation: \[ s_p = \sqrt{\frac{(n_f-1)s^2_{\mathbf{y}^f} + (n_m-1)s^2_{\mathbf{y}^m}}{n_f + n_m - 2}} \]

Model

Assumptions:

• \(y^f_{i} \sim \mathcal{N}(\mu_f, \sigma^2)\) iid
• \(y^m_{j} \sim \mathcal{N}(\mu_m, \sigma^2)\) iid

in other words:

• \(y^f_{i} = \mu_f + \epsilon_i^f\), with \(\epsilon^f_i \sim \mathcal{N}(0, \sigma^2)\) iid
• \(y^m_{j} = \mu_m + \epsilon_j^m\), with \(\epsilon^m_j \sim \mathcal{N}(0, \sigma^2)\) iid

Can we write this as a linear model: “\(\mathbf{y} = \mathbf{X}\boldsymbol{\beta} + \boldsymbol{\epsilon}\)” ?

Model

Write: \[ \mathbf{y} = \begin{pmatrix} \mathbf{y}^f \\ \mathbf{y}^m \\ \end{pmatrix} \quad \text{and} \quad \boldsymbol{\epsilon} = \begin{pmatrix} \boldsymbol{\epsilon}^f \\ \boldsymbol{\epsilon}^m \\ \end{pmatrix} \] \(\mathbf{y}\) and \(\boldsymbol{\epsilon}\) vectors of size \(n = n_f + n_m\).

Then, for \(1 \leq i \leq n\): \[ y_i = \mu_k + \epsilon_i \quad \text{with} \quad \begin{cases} \mu_k = \mu_f & \text{if } 1 \leq i \leq n_f\\ \mu_k = \mu_m & \text{if } n_f+1 \leq i \leq n_f + n_m \end{cases} \]

Regression matrix \(\mathbf{X}\) ?

Linear Model

Hence, the model can be written as: \[ \begin{pmatrix} y^f_1 \\ \vdots\\ y^f_{n_f} \\ y^m_{1} \\ \vdots\\ y^m_{n_m} \\ \end{pmatrix} = \begin{pmatrix} 1 & 0 \\ \vdots & \vdots \\ 1 & 0 \\ 0 & 1 \\ \vdots & \vdots\\ 0 & 1 \\ \end{pmatrix} \begin{pmatrix} \mu_f \\ \mu_m \\ \end{pmatrix} + \begin{pmatrix} \epsilon^f_1 \\ \vdots\\ \epsilon^f_{n_f} \\ \epsilon^m_{1} \\ \vdots\\ \epsilon^m_{n_m} \\ \end{pmatrix} \]

i.e: \[ \mathbf{y} = \mathbf{X}\boldsymbol{\beta} + \boldsymbol{\epsilon} \quad \text{with} \quad \boldsymbol{\epsilon} \sim \mathcal{N}(\mathbf{0}, \sigma^2 \mathbf{I}_n) \]

\(\to\) classical linear model with Gaussian assumption !

Linear Model - \(t\) test

\[ \mathbf{y} = \mathbf{X}\boldsymbol{\beta} + \boldsymbol{\epsilon} \quad \text{with} \quad \boldsymbol{\epsilon} \sim \mathcal{N}(\mathbf{0}, \sigma^2 \mathbf{I}_n) \]

As \[ \boldsymbol{\beta} = \begin{pmatrix} \mu_f \\ \mu_m - \mu_f \\ \end{pmatrix}, \]

the test \[ \mathcal{H}_0: \mu_f = \mu_m \quad \text{vs} \quad \mathcal{H}_1: \mu_f \neq \mu_m \]

becomes: \[ \mathcal{H}_0: \beta_2 = 0 \quad \text{vs} \quad \mathcal{H}_1: \beta_2 \neq 0 \]

\(\to\) \(t\)-test on the \(\beta_2\) coefficient !

Linear Model

fit <- lm(body_mass_g ~ sex, data = adelie)

model.matrix(fit)[1:10, ]

   (Intercept) sexmale
1            1       1
2            1       0
3            1       0
4            1       0
5            1       1
6            1       0
7            1       1
8            1       0
9            1       1
10           1       1

General Model Over-Parametrization

\[ y_{i,k} = \mu + \beta_k + \epsilon_i \]

The model has \(K+1\) parameters.

The associated regression matrix would be: \[ \mathbf{X} = (\mathbb{1}, \mathbb{1}(S_1), \dotsc, \mathbb{1}(S_K)) \]

\(\mathbf{X}\) is not of full rank: \(\mathbb{1} = \sum_{k=1}^K \mathbb{1}(S_k)\).

Exemple with two groups

\[ \mathbf{X} = \begin{pmatrix} 1 & 1 & 0 \\ \vdots & \vdots & \vdots \\ 1 & 1 & 0 \\ 1 & 0 & 1 \\ \vdots & \vdots & \vdots\\ 1 & 0 & 1 \\ \end{pmatrix} \]

We need some constraint on the parameters.

Constraint : Group 1 as reference

\[ y_{i,k} = \mu + \beta_k + \epsilon_i \]

In the previous studies, we assumed \[\beta_1 = 0\]

• group 1 has the “reference” mean \(\mu\)
• group \(k>1\) has mean \(\mu + \beta_k\)
• \(\beta_k = \mu_k - \mu_1\) is the difference between group 1 and \(k\).
• \(\mathbf{X} = (\mathbb{1}, \mathbb{1}(S_2), \dotsc, \mathbb{1}(S_K))\) is full rank.

\(\to\) This is the constraint used by default in R.

Constraint : No intercept

\[ y_{i,k} = \mu + \beta_k + \epsilon_i \]

We can assume \[\mu = 0\]

• \(\beta_k = \mu_k\) models the mean of group \(k\) directly.
• \(\mathbf{X} = (\mathbb{1}(S_1), \dotsc, \mathbb{1}(S_K))\) is full rank.
• estimators of \(\beta_k\) are changed, but not the estimator of \(\sigma^2\), nor the fitted values.

Constraint : Null Sum

\[ y_{i,k} = \mu + \beta_k + \epsilon_i \]

We can assume \[ \sum_{k = 1}^K \beta_k= 0 \qquad \text{i.e.} \qquad \beta_K = - \sum_{k = 1}^{K-1} \beta_k \]

• \(\mathbf{X} = (\mathbb{1}, \mathbb{1}(S_1) - \mathbb{1}(S_K), \dotsc, \mathbb{1}(S_{K-1}) - \mathbb{1}(S_K))\) is full rank.
• estimators of \(\mu\), \(\beta_k\) are changed, but not the estimator of \(\sigma^2\), nor the fitted values.

Constraint : Weighted Null Sum

\[ y_{i,k} = \mu + \beta_k + \epsilon_i \]

We can assume \[ \sum_{k = 1}^K n_k \beta_k= 0 \qquad \text{i.e.} \qquad \beta_K = -\frac{1}{n_K} \sum_{k = 1}^{K-1} n_k \beta_k \]

• \(\mathbf{X} = (\mathbb{1}, \mathbb{1}(S_1) - \frac{n_1}{n_K}\mathbb{1}(S_K), \dotsc, \mathbb{1}(S_{K-1}) - \frac{n_{K-1}}{n_K}\mathbb{1}(S_K))\) is full rank.
• estimators of \(\mu\), \(\beta_k\) are changed, but not the estimator of \(\sigma^2\), nor the fitted values.

One-way ANOVA - \(F\) test

\[ y_{i,k} = \mu + \beta_k + \epsilon_i \]

ANOVA Test: Has the grouping any effect at all ? I.e. \[ \begin{aligned} \mathcal{H}_0&: \beta_{1} = \cdots = \beta_K = 0 &\text{vs}\\ \mathcal{H}_1&: \exists~k\in \{1, \dotsc, K\} ~|~ \beta_k \neq 0 \end{aligned} \]

No matter which set of constraint we choose, the \(F\) statistic associated with the one factor ANOVA is:

\[ F = \frac{ \|\hat{\mathbf{y}} - \bar{y}\mathbb{1}\|^2 / (K - 1) }{ \|\mathbf{y} - \hat{\mathbf{y}}\|^2 / (n - K) } \underset{\mathcal{H}_0}{\sim} \mathcal{F}^{K - 1}_{n - K}. \]

\(\to\) Global \(F\) test on the coefficients of the linear regression.

ANOVA = “ANalyse Of VAriance”

Model

Observation \(i\) of group \(k\) in factor 1 and group \(l\) in factor 2 is written as: \[ y_{i,k,l} = \mu + \alpha_k + \beta_l + \gamma_{kl} + \epsilon_i \]

• \(\mu\) the “grand mean” (intercept)
• \(\alpha_k\) the specific effect of group \(k\) in factor 1 (\(1 \leq k \leq K\)).
• \(\beta_l\) the specific effect of group \(l\) in factor 2 (\(1 \leq l \leq L\)).
• \(\gamma_{kl}\) the interaction effect of groups \(k\) and \(l\).

The full model has regression matrix: \[ {\small \mathbf{X} = (\mathbb{1}, \mathbb{1}_{S_1}, \dotsc, \mathbb{1}_{S_K}, \mathbb{1}_{T_1}, \dotsc, \mathbb{1}_{T_L}, \mathbb{1}_{S_1, T_1}, \dotsc, \mathbb{1}_{S_K, T_L} ) } \]

\(\to\) \(1 + K + L + KL\) predictors. \(\to\) not identifiable !

Constraint : Groups 1 as reference

\[ y_{i,k,l} = \mu + \alpha_k + \beta_l + \gamma_{kl} + \epsilon_i \]

Assume, for all \(1\leq k \leq K\) and \(1\leq l \leq L\): \[ \alpha_1 = 0, \qquad \beta_1 = 0, \qquad \gamma_{1,l} = 0, \qquad \gamma_{k,1} = 0. \]

• \(y_{i,1,1}\) has mean \(\mu\) (reference)
• \(y_{i,1,l}\) has mean \(\mu + \beta_l\)
• \(y_{i,k,1}\) has mean \(\mu + \alpha_k\)
• \(y_{i,k,l}\) has mean \(\mu + \alpha_k + \beta_l + \gamma_{k,l}\) (\(k,l>1\))
• \(1 + K + L\) constraints \(\to\) \(KL\) free parameters
• \(\mathbf{X} = (\mathbb{1}, \mathbb{1}_{S_2}, \dotsc, \mathbb{1}_{S_K}, \mathbb{1}_{T_2}, \dotsc, \mathbb{1}_{T_L}, \mathbb{1}_{S_2, T_2}, \dotsc, \mathbb{1}_{S_K, T_L})\) full

\(\to\) This is the constraint used by default in R.

Other constraint

\[ y_{i,k,l} = \mu + \alpha_k + \beta_l + \gamma_{kl} + \epsilon_i \]

Interactions only: for all \(1\leq k \leq K\) and \(1\leq l \leq L\): \[ \mu = 0,\qquad \alpha_k = 0, \qquad \beta_l = 0. \]

Null sum: \[ \sum_{k=1}^K\alpha_k = 0, \qquad \sum_{l=1}^L\beta_l = 0, \qquad \sum_{k=1}^K\gamma_{kl} = 0, \qquad \sum_{l=1}^L\gamma_{kl} = 0. \]

\(\to\) changes the estimation of the coefficients, but not of the variance, nor the predicted values.

Nested F Test - Reminder

Recall that (from CM5), in general, to test: \[ \begin{aligned} \mathcal{H}_0&: \beta_{p_0 + 1} = \cdots = \beta_p = 0 &\text{vs}\\ \mathcal{H}_1&: \exists~k\in \{p_0+1, \dotsc, p\} ~|~ \beta_k \neq 0 \end{aligned} \]

We can use the statistic: \[ F = \frac{ \|\hat{\mathbf{y}} - \hat{\mathbf{y}}_0\|^2 / (p - p_0) }{ \|\mathbf{y} - \hat{\mathbf{y}}\|^2 / (n - p) } = \frac{n - p}{p - p_0}\frac{RSS_0 - RSS}{RSS} \underset{\mathcal{H}_0}{\sim} \mathcal{F}^{p - p_0}_{n - p}. \]

Writing \(\mathbf{X}_{p_0} = (\mathbf{X}_1, \dotsc, \mathbf{X}_{p_0})\), this is the same as testing: \[ \begin{aligned} \mathcal{H}_0&: \mathbf{P}^{\mathbf{X}_{p_0}}\mathbf{y} = \mathbf{P}^{\mathbf{X}}\mathbf{y} &\text{vs}\\ \mathcal{H}_1&: \mathbf{P}^{\mathbf{X}_{p_0}}\mathbf{y} \neq \mathbf{P}^{\mathbf{X}}\mathbf{y} \end{aligned} \]

Anova model

Model: \[ y_{i,k,l} = \mu + \alpha_k + \beta_l + \gamma_{kl} + \epsilon_i \] with, for all \(1\leq k \leq K\) and \(1\leq l \leq L\): \[ \alpha_1 = 0, \qquad \beta_1 = 0, \qquad \gamma_{1,l} = 0, \qquad \gamma_{k,1} = 0. \] (\(\to KL\) free parameters in total.)

Tests:

• Is factor 1 useful ?
• Is factor 2 useful ?
• Is the interaction between factor 1 and 2 useful ?

Nested F Test - Anova - Factor 1

To test: \[ \begin{aligned} \mathcal{H}_0&: \mathbf{P}^{\mathbf{X}_{\mu}}\mathbf{y} = \mathbf{P}^{\mathbf{X}_{\mu,\alpha}}\mathbf{y} &\text{v.s.}&& \mathcal{H}_1&: \mathbf{P}^{\mathbf{X}_{\mu}}\mathbf{y} \neq \mathbf{P}^{\mathbf{X}_{\mu,\alpha}}\mathbf{y} \end{aligned} \] i.e. \(\mathcal{H}_0\): “factor 1 does not bring any information compared to intercept alone”.

We can use: \[ \begin{aligned} F &= \frac{n - KL}{K - 1}\frac{RSS_{\mu} - RSS_{\mu, \alpha}}{RSS}\\ &= \frac{R(\alpha | \mu) / (K - 1)}{\hat{\sigma}^2}\\ &\underset{\mathcal{H}_0}{\sim} \mathcal{F}^{K - 1}_{n - KL}. \end{aligned} \]

Nested F Test - Anova - Factor 2

To test: \[ \begin{aligned} \mathcal{H}_0&: \mathbf{P}^{\mathbf{X}_{\mu,\alpha}}\mathbf{y} = \mathbf{P}^{\mathbf{X}_{\mu,\alpha, \beta}}\mathbf{y} &\text{v.s.}&& \mathcal{H}_1&: \mathbf{P}^{\mathbf{X}_{\mu,\alpha}}\mathbf{y} \neq \mathbf{P}^{\mathbf{X}_{\mu,\alpha,\beta}}\mathbf{y} \end{aligned} \] i.e. \(\mathcal{H}_0\): “factor 2 does not bring any information compared to factor 1”.

We can use: \[ \begin{aligned} F &= \frac{n - KL}{L - 1}\frac{RSS_{\mu, \alpha} - RSS_{\mu, \alpha, \beta}}{RSS}\\ &= \frac{R(\beta | \mu, \alpha) / (L - 1)}{\hat{\sigma}^2}\\ &\underset{\mathcal{H}_0}{\sim} \mathcal{F}^{L - 1}_{n - KL}. \end{aligned} \]

Nested F Test - Anova - Interaction

To test: \[ \begin{aligned} \mathcal{H}_0&: \mathbf{P}^{\mathbf{X}_{\mu,\alpha,\beta}}\mathbf{y} = \mathbf{P}^{\mathbf{X}_{\mu,\alpha, \beta,\gamma}}\mathbf{y} &\text{v.s.}&& \mathcal{H}_1&: \mathbf{P}^{\mathbf{X}_{\mu,\alpha,\gamma}}\mathbf{y} \neq \mathbf{P}^{\mathbf{X}_{\mu,\alpha,\beta,\gamma}}\mathbf{y} \end{aligned} \] i.e. \(\mathcal{H}_0\): “Interactions do not bring any information compared to factors 1 and 2”.

We can use: \[ \begin{aligned} F &= \frac{n - KL}{(L - 1)(K - 1)}\frac{RSS_{\mu, \alpha, \beta} - RSS_{\mu, \alpha, \beta, \gamma}}{RSS}\\ &= \frac{R(\gamma | \mu, \alpha, \beta) / ((K - 1)(L - 1))}{\hat{\sigma}^2}\\ &\underset{\mathcal{H}_0}{\sim} \mathcal{F}^{(K - 1)(L - 1)}_{n - KL}. \end{aligned} \]

Ancova - One Regressor - Model

For any observation \(1 \leq i \leq n\) in group \(1 \leq k \leq K\), we write: \[ y_{i,k} = \mu + \alpha_k + (\beta + \gamma_k) \times x_{i} + \epsilon_i \qquad \epsilon_i\sim\mathcal{N}(0, \sigma^2)~iid \]

Identifiability: \(\alpha_1 = \gamma_1 = 0\).

\[ \text{With $K=2$:}\qquad \mathbf{X} = \begin{pmatrix} 1 & 0 & x_1 & 0\\ \vdots & \vdots & \vdots & \vdots\\ 1 & 0 & x_{n_1} & 0 &\\ 1 & 1 & x_{n_1 + 1} & x_{n_1 + 1} \\ \vdots & \vdots & \vdots & \vdots\\ 1 & 1 & x_{n_1 + n_2} & x_{n_1 + n_2}\\ \end{pmatrix} \]