\(\newcommand{\VS}{\quad \mathrm{VS} \quad}\) \(\newcommand{\and}{\quad \mathrm{and} \quad}\) \(\newcommand{\E}{\mathbb E}\) \(\newcommand{\P}{\mathbb P}\) \(\newcommand{\Var}{\mathbb V}\) \(\newcommand{\Cov}{\mathrm{Cov}}\) \(\newcommand{\1}{\mathbf 1}\)
⬅ Previous: Selection · Next: Introduction to the Generalized Linear Model ➡
Previously, we considered that:
We still assume \(Y\) is quantitative, but explanatory variables can be qualitative and/or quantitative.
We’ll see that these situations reduce to the case of the previous chapter.
\(X\) is qualitative, with \(I\) modalities \(A_1, \dots, A_I\). We use one-hot encoding: \(X_{ki}=\1\{X_k=A_i\}\)
That is, for individual \(k\), \(X_{k\cdot} = (\1\{X_k=A_1\}, \dots, \1\{X_k=A_I\}) \in \{0,1\}^I\)
Or, using previous notations, \(X = (X^{(1)}, \dots, X^{(I)})\) where
\[X^{(i)}= \begin{pmatrix} X^{(i)}_1 \\ \vdots \\ X^{(i)}_n \end{pmatrix}\]
If there are \(3\) categories e.g. blue, orange, green, we replace the column \(X\) by \(3\) columns
Example with \(I=3\) categories and \(n=5\) individuals \[ \begin{pmatrix} \color{blue}{\mathrm{blue}} \\ \mathrm{green} \\ \color{blue}{\mathrm{blue}} \\ \mathrm{orange} \\ \mathrm{orange} \\ \end{pmatrix} \quad \text{becomes} \quad X=\begin{pmatrix} \color{blue}{\mathrm{1}} & \mathrm{0} & \mathrm{0}\\ \mathrm{0} & \mathrm{0} & \mathrm{1}\\ \color{blue}{\mathrm{1}} & \mathrm{0} & \mathrm{0}\\ \mathrm{0} & \mathrm{1} & \mathrm{0}\\ \mathrm{0} & \mathrm{1} & \mathrm{0}\\ \end{pmatrix} \]
We seek to explain \(Y\) using a single qualitative variable \(X\).
We observe \((Y, X)\), and define
\[N_i = \sum_{k=1}^n \mathbf 1\{X_k=i\} \and \overline Y_i = \frac{1}{N_i}\sum_{k=1}^n Y_k \mathbf 1\{X_k=i\}\]
Total mean:
\(\overline Y = \frac{1}{n}\sum_{k=1}^n Y_k = \frac{1}{n}\sum_{i=1}^I N_i \overline Y_i\)
For \(k = 1, \ldots, n\)
\[Y_k = \sum_{i=1}^I \mu_i \mathbf{1}\{X_k=A_i\} + \varepsilon_k\]
where \(\E[\varepsilon_k]=0\), \(\Cov(\varepsilon_k,\varepsilon_l)=\sigma^2\1\{k=l\}\)
We rewrite the model \(Y_k = \sum_{i=1}^I \mu_i \mathbf{1}\{X_k=A_i\} + \varepsilon_k\) as
\[Y = X\mu + \varepsilon\]
If we want to model the constant (intercept) and \(I\) modalities, we assume
\[Y_k = m+\sum_{i=2}^I \alpha_i \mathbf{1}\{X_k=A_i\} + \varepsilon_k\]
Equivalently,
\(Y = m \mathbf{1} + \sum_{i=2}^I \alpha_i\mathbf{1}_{A_i}+\varepsilon\)
\(Y_k = m+\sum_{i=2}^I \alpha_i \mathbf{1}\{X_k=A_i\} + \varepsilon_k\)
\(\sum_{i=1}^I \mathbf{1}_{A_i} = \mathbf{1}\)
We have \(I+1\) parameters instead of \(I\). We have to set one constraint:
Interpretation: gives expectation \(\E[Y_k] = \mu_1\) and coefficients \(\alpha_i= \mu_i - \mu_1\)
Interpretation: gives coefficients \(\alpha_i= \mu_i\)
\(Y\) = weight Loss, \(X\) = Exercise programme: \(I = 4\) modalities, \(N_i = 18\) each, \(n = 72\).
Programme 2 looks best, 1 and 4 almost tied.
Data of the example below: Loss ~ Exercise, \(I=4\) groups of \(N_i=18\).
1 2 3 4
9.6111 16.7778 13.5000 9.0000 that is \(\overline{Y}_1 = 9.6111\), \(\overline{Y}_2 = 16.7778\), \(\overline{Y}_3 = 13.5\), \(\overline{Y}_4 = 9\).
Estimate Std. Error t value Pr(>|t|)
(Intercept) 9.6111 0.7992 12.025 < 2e-16 ***
Exercise2 7.1667 1.1303 6.340 2.14e-08 ***
Exercise3 3.8889 1.1303 3.441 0.000996 ***
Exercise4 -0.6111 1.1303 -0.541 0.5905
Residual standard error: 3.391 on 68 degrees of freedomR prints no Exercise1: the column \(\mathbf{1}_{A_1}\) was dropped, and \(\hat m = \overline{Y}_1 = 9.6111\) is exactly the mean of that removed modality.
The others are differences to it: \(\hat\alpha_i = \overline{Y}_i - \overline{Y}_1\), e.g. \(16.7778 - 9.6111 = 7.1667\).
Estimate Std. Error t value Pr(>|t|)
Exercise1 9.6111 0.7992 12.025 < 2e-16 ***
Exercise2 16.7778 0.7992 20.992 < 2e-16 ***
Exercise3 13.5000 0.7992 16.891 < 2e-16 ***
Exercise4 9.0000 0.7992 11.261 < 2e-16 ***
Residual standard error: 3.391 on 68 degrees of freedomThe four modalities are back, and \(\hat\alpha_i = \overline{Y}_i\) — the same four numbers as tapply.
Both fits are the same model: same \(\widehat{Y}_k = \overline{Y}_i\) for \(X_k = A_i\), same residuals, same \(\hat\sigma = 3.391\) on \(n-I=68\) degrees of freedom. Only the parameterization changes.
So Exercise2 = 7.1667 never means “group 2 has mean \(7.17\)”: it means group 2 sits \(7.17\) above group 1.
Whichever the model we choose (with constant or not), estimation of \(\mu_i=\E[Y_k|X_k=A_i]\) is the same. Same for \(\Var(\varepsilon_k)=\sigma^2\).
Proposition
In category \(i\), OLS estimation of \(\mu_i\) leads to:
An unbiased estimator of \(\sigma^2\) is
We write \(\mathbf 1_i = (\1\{X_1=A_i\}, \dots, \1\{X_n=A_i\})^T\), and \(E = span(\mathbf 1_1, \dots, \mathbf 1_I)\)
Then, the OLS is exactly
\(\hat \mu = \mathrm{argmin}_{\mu' \in E}\|Y - \mu'\|_2^2 = P_E Y\)
For \(\hat \sigma\), \(\mathbb E[((I_n-P_E)\varepsilon)^T(I_n-P_E)\varepsilon]= \sigma^2 (n-I)\)
We want to test \(H_0: \mu_1 = \cdots = \mu_I\).
This is a linear constraints test! (See previous chapters)
Proposition
If \(\varepsilon \sim N(0, \sigma^2 I_n)\), then under \(H_0: \mu_1 = \cdots = \mu_I\):
\[F = \frac{SSB/(I-1)}{SSW/(n-I)} \sim F(I-1, n-I)\]
Critical region at level \(\alpha\): \(CR_\alpha = \{F > f_{I-1,n-I}(1-\alpha)\}\)
One-way ANOVA (a single factor \(A\) with \(I\) groups). Pull the group means apart (separation) or add within-group noise. The partition \(SST = SSB + SSW\) shifts: \(F = \tfrac{SSB/(I-1)}{SSW/(n-I)}\) is large exactly when the between share dominates the within share.
Using the previous projection \(P_E = \sum_{i=1}^I \frac{1}{N_i}\mathbf{1}_i\mathbf{1}_i^T\), we have the decomposition
\[Y - P_0Y = Y-P_E Y + P_E Y - P_0 Y\]
Under \(H_0: \mu_1 = \cdots = \mu_I\), there are \(I-1\) constraints to test.
\[F = \frac{n-I}{I-1} \cdot \frac{SSR_c - SSR}{SSR}\]
We show that
In R: anova(lm(Y~A))
Warning
The previous analysis of variance test tests equality of means between modalities not equality of variances
It is valid under the assumption \(\varepsilon \sim N(0, \sigma^2 I_n)\).
How to test equality of variances in each modality:
In R: leveneTest or bartlett.test from car library
If factor \(A\) is significant, we want to know more:
which modality(ies) differ from the others?
We want to perform all tests:
\[H_{0}^{i,j}: \mu_i = \mu_j \quad \text{vs} \quad H_{1}^{i,j}: \mu_i \neq \mu_j\]
for all \(i \neq j\) in \(\{1, \ldots, I\}\), corresponding to \(I(I-1)/2\) tests.
Perform all Student’s t-tests for mean comparison (1 constraint), each at level \(\alpha\).
Problem: Given the number of tests, this would lead to many false positives.
False positives are a well-known problem in multiple testing.
Solution: Apply a correction to the decision rule, e.g.
For one-way ANOVA: Tukey’s test addresses the problem.
\[Q = \max_{(i,j)} \frac{|\overline{Y}_i - \overline{Y}_j|}{\hat{\sigma}\sqrt{\tfrac{1}{2}\left(\frac{1}{N_i} + \frac{1}{N_j}\right)}}\]
Distribution of Tukey’s Test Statistic
Under \(H_0: \mu_1 = \cdots = \mu_I\) and \(\varepsilon \sim N(0, \sigma^2 I_n)\), \[Q \sim Q_{I,n-I},\] the Tukey distribution with \((I, n-I)\) degrees of freedom — exactly if the \(N_i\) are equal, approximately otherwise.
To test each \(H_0^{i,j}: \mu_i = \mu_j\), we use the critical regions:
\[CR_\alpha^{i,j} = \left\{|\overline{Y}_i - \overline{Y}_j| > \frac{\hat{\sigma}}{\sqrt{2}} \sqrt{\frac{1}{N_i} + \frac{1}{N_j}} \cdot Q_{I,n-I}(1-\alpha)\right\}\]
where \(Q_{I,n-I}(1-\alpha)\) denotes the \((1-\alpha)\) quantile of a tukey distribution \(Q(I,n-I)\) of degrees \((I, n-I)\).
The form of the previous \(CR_\alpha^{i,j}\) ensures that: \[\mathbb{P}_{\mu_1 = \cdots = \mu_I}\left(\bigcup_{(i,j)} \{Y \in CR_\alpha^{i,j}\}\right) = \alpha\]
Interpretation: If all null hypotheses \(H_0^{i,j}\) are true (\(\mu_1 = \cdots = \mu_I\)), then the probability of concluding at least one \(H_1^{i,j}\) equals \(\alpha\).
This is the simultaneous Type I error rate.
With Tukey’s test, two significantly different means are truly different, not just due to false positives.
In R: TukeyHSD
All true means are equal, so every colour is a false positive: left, the FWER by Monte-Carlo (hover for the %); middle, the pairs flagged by the naive test (red) and by Tukey (green); right, this experiment’s means against the two thresholds.
Homoscedasticity Test:
We get
Levene’s Test for Homogeneity of Variance (center = median)
Df F value Pr(>F)
group 3 0.6527 0.584
68We get
Response: Loss
Df Sum Sq Mean Sq F value Pr(>F)
Exercise 3 712.56 237.519 20.657 1.269e-09 ***
Residuals 68 781.89 11.498The means are significantly different. We read in particular:
We finish by analyzing the mean differences more precisely:
Exercise
diff lwr upr p adj
2-1 7.1666667 4.1897551 10.1435782 0.0000001
3-1 3.8888889 0.9119773 6.8658005 0.0053823
4-1 -0.6111111 -3.5880227 2.3658005 0.9487355
3-2 -3.2777778 -6.2546894 -0.3008662 0.0252761
4-2 -7.7777778 -10.7546894 -4.8008662 0.0000000
4-3 -4.5000000 -7.4769116 -1.5230884 0.0009537Conclusion: All differences are significant, except between exercise 4 and exercise 1.
We observe \(Y\) and explanatory variables \(X^{(1)}, X^{(2)}\)
\[\begin{aligned} Y_k &= m + \sum_{i=1}^I \alpha_i\1\{X^{(1)}_k = A_i\} + \sum_{j=1}^J \beta_j\1\{X^{(2)}_k = B_j\} \\ &+ \sum_{i=1}^I\sum_{j=1}^J \gamma_{ij}\1\{X^{(1)}_k = A_i\}\1\{X^{(2)}_k = B_j\} + \varepsilon_k \end{aligned}\]
In other words, in modality \(A_i\) and \(B_j\),
\(Y_k =m+\alpha_i + \beta_j + \gamma_{ij}+ \varepsilon_k\)
This is a model of the type \(Y_k = \mu_{ij} + \varepsilon_k\).
In modalities \((i,j)\), \(\E[Y_k] = \mu_{ij} = m + \alpha_i + \beta_j + \gamma_{ij}\)
\(Y\): employee satisfaction
\(A\): schedule type (flexible or fixed)
\(B\): training level (basic or advanced)
We can imagine:
\(Y\): plant yield
\(A\): fertilizer type (1 or 2)
\(B\): water quantity (low, medium, high)
We can imagine:
Maybe interaction is so strong that the effect due to A seems absent!
In modalities \((i,j)\), \(\E[Y_k] = \mu_{ij} = m + \alpha_i + \beta_j + \gamma_{ij}\)
Initial two-factor ANOVA problem: \(I \times J\) parameters \(\mu_{ij}\).
Now: \(1 + I + J + IJ\) parameters (\(m\), \(\alpha_i\), \(\beta_j\), \(\gamma_{ij}\)).
Therefore, we need \(1 + I + J\) constraints for identifiability:
\(\sum_{i=1}^{I} \alpha_i = 0 \and \sum_{j=1}^{J} \beta_j = 0\)
\(\sum_{i=1}^{I} \gamma_{ij} = 0 \and \sum_{j=1}^{J} \gamma_{ij} = 0\)
These constraints ensure model identifiability by removing the redundant parameters that cause multicollinearity issues.
The complete model (with interaction) is launched with the command:
With these constraints, the parameter interpretation is as follows:
Intercept: \(m = \mu_{11}\)
Main effect A: \(\alpha_i = \mu_{i1} - \mu_{11}\)
Main effect B: \(\beta_j = \mu_{1j} - \mu_{11}\)
Interaction: \(\gamma_{ij} = \mu_{ij} - \mu_{i1} - \mu_{1j} + \mu_{11}\)
Intercept: \(m = \mu_{11}\)
Main effect A: \(\alpha_i = \mu_{i1} - \mu_{11}\)
Main effect B: \(\beta_j = \mu_{1j} - \mu_{11}\)
Interaction: \(\gamma_{ij} = \mu_{ij} - \mu_{i1} - \mu_{1j} + \mu_{11}\)
The choice of constraints does not affect the estimation of the expectation of \(Y\) in each crossed modality \(A_i \cap B_j\).
Proposition
Whatever the linear constraints chosen, the OLS leads to, for all \(i = 1, \ldots, I\), \(j = 1, \ldots, J\) and \(k = 1, \ldots, N_{ij}\), if \(X^{(1)}_{k}=A_i\) and \(X_k^{(2)}=B_j\):
\[\widehat{Y}_{k} = \overline{Y}_{ij}:= \frac{1}{N_{ij}} \sum_{k=1}^{n} Y_{k}\1\{X_k^{(1)}=A_i \and X_k^{(2)}=B_j\}\]
and to the estimation of the residual variance:
\[\hat{\sigma}^2 = \frac{1}{n - IJ} \sum_{i=1}^I\sum_{j=1}^J\sum_{k=1}^n (Y_{k} - \overline{Y}_{ij})^2\1\{X_k^{(1)}=A_i \and X_k^{(2)}=B_j\}\]
Estimate Std. Error t value Pr(>|t|)
(Intercept) 12.000 0.8573 13.997 1.09e-14 ***
fertilizer2 6.000 1.2124 4.949 2.69e-05 ***
watermedium 4.000 1.2124 3.299 0.00250 **
waterhigh 10.000 1.2124 8.248 3.31e-09 ***
fertilizer2:watermedium -3.000 1.7146 -1.750 0.09041 .
fertilizer2:waterhigh -11.000 1.7146 -6.415 4.39e-07 ***
Residual standard error: 2.1 on 30 degrees of freedomEverything is read against the cell \((A_1, B_1)\): \(\hat m = \overline{Y}_{11} = 12\), \(\hat\alpha_2 = \overline{Y}_{21} - \overline{Y}_{11} = 6\), \(\hat\beta_3 = \overline{Y}_{13} - \overline{Y}_{11} = 10\), \(\hat\gamma_{23} = \overline{Y}_{23} - \overline{Y}_{21} - \overline{Y}_{13} + \overline{Y}_{11} = -11\).
Summing them back gives the crossed modality:
\(\hat\mu_{23} = \hat m + \hat\alpha_2 + \hat\beta_3 + \hat\gamma_{23} = 12 + 6 + 10 - 11 = 17 = \overline{Y}_{23}\)
Estimate Std. Error t value Pr(>|t|)
(Intercept) 17.3333 0.3500 49.524 < 2e-16 ***
fertilizer1 -0.6667 0.3500 -1.905 0.06644 .
water1 -2.3333 0.4950 -4.714 5.22e-05 ***
water2 0.1667 0.4950 0.337 0.73870
fertilizer1:water1 -2.3333 0.4950 -4.714 5.22e-05 ***
fertilizer1:water2 -0.8333 0.4950 -1.684 0.10260 \(\hat m = 17.33\) is now the grand mean; the last level of each factor is not printed, being minus the sum of the others: \(\hat\alpha_2 = +0.6667\), \(\hat\beta_3 = +2.1667\), \(\hat\gamma_{23} = -3.1667\). Same fit, same \(\hat\sigma = 2.1\):
\(\hat\mu_{23} = 17.3333 + 0.6667 + 2.1667 - 3.1667 = 17 = \overline{Y}_{23}\)
Yield (q/ha) in hydroponics, \(E\) = fertiliser, \(W\) = water, \(40\) per cell. Truth: \(Y = 50\,EW + \varepsilon\).
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.23 0.791 -0.291 0.77149
fert 0.21 1.118 0.188 0.85125
water 0.45 1.118 0.402 0.68787
fert:water 49.83 1.581 31.515 < 2e-16 ***
Residual standard error: 5 on 156 degrees of freedom| \(\overline{Y}_{ij}\) | no water | water |
|---|---|---|
| no fertiliser | −0.23 | 0.22 |
| fertiliser | −0.02 | 50.26 |
The two factors work better together than the sum of their separate effects.
Loneliness score (0–40), \(C\) = dog, \(H\) = cat, \(60\) per cell. Truth: \(Y = 40 - 20\,C - 20\,H + 20\,CH + \varepsilon\), i.e. \(Y = 40 - 20\max(C,H)\).
Estimate Std. Error t value Pr(>|t|)
(Intercept) 39.70 1.033 38.439 < 2e-16 ***
dog -19.73 1.461 -13.508 < 2e-16 ***
cat -19.41 1.461 -13.289 < 2e-16 ***
dog:cat 19.78 2.066 9.576 < 2e-16 ***
Residual standard error: 8 on 236 degrees of freedom| \(\overline{Y}_{ij}\) | no cat | cat |
|---|---|---|
| no dog | 39.70 | 20.29 |
| dog | 19.97 | 20.34 |
The two factors do the same job, so having one makes the other less useful.
Survival (months), \(T\) = treatment, \(G\) = variant carrier, \(50\) per cell. Truth: \(Y = 24 - 12\,T - 12\,G + 24\,TG + \varepsilon\).
Estimate Std. Error t value Pr(>|t|)
(Intercept) 23.75 0.849 27.990 < 2e-16 ***
treat -11.77 1.200 -9.808 < 2e-16 ***
variant -11.51 1.200 -9.592 < 2e-16 ***
treat:variant 23.81 1.697 14.030 < 2e-16 ***
Residual standard error: 6 on 196 degrees of freedom| \(\overline{Y}_{ij}\) | non-carrier | carrier |
|---|---|---|
| placebo | 23.75 | 12.24 |
| treatment | 11.98 | 24.28 |
Crossover: the effect of one factor changes direction depending on the other.
Drag a cell mean (solid \(\mu_{ij}\), dashed \(\hat Y_{ij}\) with noise); i / j select the active cell, constraints switches the coding. The table holds the whole decomposition, each \(\gamma_{ij}\) green at \(0\). In R: interaction.plot(B, A, Y).
| \(\gamma_{ij}\) | \(j=1\) | \(j=2\) | \(j=3\) | \(\alpha_i\) |
| \(i=1\) | – | – | – | – |
| \(i=2\) | – | – | – | – |
| \(\beta_j\) | – | – | – | \(m\) – |
We test for interaction first. Under \(H_0^{(AB)}\), the model reduces to the additive model \(Y_{ijk} = m + \alpha_i + \beta_j + \varepsilon_{ijk}\)
Do we have \(\gamma_{ij} = 0\) for all \(i, j\)?
We want to test for the presence of an interaction:
\[H_0^{(AB)}: \gamma_{ij} = 0 \text{ for all } i, j\]
If we conclude \(H_0^{(AB)}\) (accept the null hypothesis), we then want to test the marginal effects:
\[H_0^{(A)}: \alpha_i = 0 \text{ for all } i \and H_0^{(B)}: \beta_j = 0 \text{ for all } j\]
\[H_0^{(AB)}: \gamma_{ij} = 0 \text{ for all } i, j\]
Warning
If we reject \(H_0^{(AB)}\), it makes no sense to test whether A or B have an effect: they have one through their interaction.
These tests reduce to constraint tests in the regression model.
We assume that “the design is balanced”: this means that \(N_{ij}:=N\) does not depend on \(i\) or \(j\). (Otherwise, everything becomes complicated).
In this case, we have the analysis of variance formula:
\[S_T^2 = S_A^2 + S_B^2 + S_{AB}^2 + S_R^2\]
\(S_T^2 = \sum_{k=1}^{n} (Y_{k} - \overline{Y})^2\): total sum of squares
\(S_A^2 = \sum_{i=1}^{I} \sum_{j=1}^{J} N_{ij} (\overline{Y}_{i.} - \overline{Y})^2\): \(S^2_{between}\) in the case of one-factor ANOVA where the factor is \(A\)
\(S_B^2 = \sum_{i=1}^{I} \sum_{j=1}^{J} N_{ij} (\overline{Y}_{.j} - \overline{Y})^2\): \(S^2_{between}\) in the case of one-factor ANOVA where the factor is \(B\)
\(S_{AB}^2 = \sum_{i=1}^{I} \sum_{j=1}^{J} N_{ij} (\overline{Y}_{ij} - \overline{Y}_{i.} - \overline{Y}_{.j} + \overline{Y})^2\) quantifies the interaction
\(S_R^2 = \sum_{i=1}^{I} \sum_{j=1}^{J} \sum_{k=1}^{N_{ij}} (Y_{k} - \overline{Y}_{ij})^2\1\{X^{(1)}_k=A_i \and X^{(2)}_k=B_j\}\): \(S_{within}\) in one-factor ANOVA
For testing the interaction, \(H_0^{(AB)}: \gamma_{ij} = 0\) for all \(i, j\)
We use the linear constraint test statistic \(F = \frac{n-p}{q} \frac{(SSR_c - SSR)}{SSR}\)
where \(SSR_c\) corresponds to the sum of squares of the residuals in the space \(\gamma_{ij}=0\) for all \((i,j)\).
\[F^{(AB)} = \frac{S_{AB}^2/\big((I-1)(J-1)\big)}{S_R^2/(n-IJ)}\]
When \(\varepsilon \sim \mathcal N(0, \sigma^2 I_n)\), \(F^{(AB)} \sim \mathcal F((I-1)(J-1), n-IJ)\)
For testing the main effect of \(A\), \(H_0^{(A)}: \alpha_i = 0\) for all \(i\):
We use \(F = \frac{n-p}{q} \frac{(SSR_c - SSR)}{SSR}\)
where \(SSR_c\) corresponds to the sum of squares of the residuals in the space \(\alpha_i=0\) for all \(i\).
\[F^{(A)} = \frac{S_A^2/(I-1)}{S_R^2/(n-IJ)}\]
When \(\varepsilon \sim \mathcal{N}(0, \sigma^2 I_n)\), \(F^{(A)} \sim \mathcal{F}(I-1, n-IJ)\)
For testing the main effect of \(B\), \(H_0^{(B)}: \beta_j = 0\) for all \(j\):
We use \(F = \frac{n-p}{q} \frac{(SSR_c - SSR)}{SSR}\)
where \(SSR_c\) corresponds to the sum of squares of the residuals in the space \(\beta_j=0\) for all \(j\).
\[F^{(B)} = \frac{S_B^2/(J-1)}{S_R^2/(n-IJ)}\]
When \(\varepsilon \sim \mathcal{N}(0, \sigma^2 I_n)\), \(F^{(B)} \sim \mathcal{F}(J-1, n-IJ)\)
In software, these tests are summarized in a table as shown below.
In R: anova(lm(Y ~ A*B))
| Source | df | Sum Sq | Mean Sq | F value | Pr(>F) |
|---|---|---|---|---|---|
| \(A\) | \(I-1\) | \(S_A^2\) | \(S_A^2/(I-1)\) | \(F^{(A)}\) | … |
| \(B\) | \(J-1\) | \(S_B^2\) | \(S_B^2/(J-1)\) | \(F^{(B)}\) | … |
| \(A:B\) | \((I-1)(J-1)\) | \(S_{AB}^2\) | \(S_{AB}^2/((I-1)(J-1))\) | \(F^{(AB)}\) | … |
| Residuals | \(n-IJ\) | \(S_R^2\) | \(S_R^2/(n-IJ)\) |
Fisher tests are based on the assumption \(\varepsilon \sim N(0, \sigma^2 I_n)\)
Test equality of variances
Form the ANOVA table (independent of chosen constraints)
If the \(AB\) interaction is significant: don’t change anything
If the interaction is not significant: analyze the marginal effects of A and B
lm(Y ~ A + B)Once effects are identified: perform post-hoc analysis by examining differences between (crossed) modalities more closely, using Tukey’s test as in one-factor ANOVA
We seek to explain \(Y\) using \(k\) qualitative variables \(A=X^{(1)}\), \(B=X^{(2)}\), \(C=X^{(3)}\), …
We can assume that \(\E[Y_k]\) depends on each factor and their interactions:
Warning
This results in \(2^k - 1\) possible effects for \(k\) factors.
Each effect can be tested using linear constraint tests. However, this approach presents several challenges:
In practice, choices are made to include only a limited number of effects and interactions in the analysis.
Run time (min) of a job: algo (new vs old), ssd (SSD vs HDD), ram (32 vs 8 GB), \(25\) per cell.
Estimate Std. Error t value Pr(>|t|)
(Intercept) 57.76 0.776 74.473 < 2e-16 ***
algo -5.79 0.776 -7.461 2.72e-12 ***
ssd -8.06 0.776 -10.397 < 2e-16 ***
ram -5.94 0.776 -7.665 8.09e-13 ***
Residual standard error: 5.484 on 196 degrees of freedomThree switches, three savings, read anywhere: \(-5.8\) min for the new algorithm, \(-8.1\) for the SSD, \(-5.9\) for the RAM.
Estimate Std. Error t value Pr(>|t|)
(Intercept) 59.72 0.810 73.715 < 2e-16 ***
algo -9.72 1.025 -9.486 < 2e-16 ***
ssd -12.00 1.025 -11.708 < 2e-16 ***
ram -5.94 0.725 -8.203 3.11e-14 ***
algo:ssd 7.87 1.449 5.430 1.66e-07 ***
Residual standard error: 5.124 on 195 degrees of freedomThe new algorithm saves \(9.7\) min on a HDD but only \(9.72 - 7.87 \approx 1.9\) on an SSD: it was mostly avoiding disk reads. The RAM keeps its own \(-5.9\) whatever the rest, so it stays additive.
Most general situation: we seek to explain Y using both quantitative and qualitative variables.
Quantitative variable effects: Through each regression coefficient \(\beta_j\) associated with each variable
Factor effects and interactions: As in ANOVA analysis
Mixed interactions: Effects of interactions between factors and quantitative variables
\[Y = \beta_0 + \beta_1 \1\{A = \text{Lab}\} + \beta_2 X^{(1)} + \beta_3 \1\{A = \text{Lab}\}X^{(1)} + \varepsilon\]
\(\beta_3\): one extra km costs a Labrador more grams than a Chihuahua — the two slopes differ.
Each effect can be tested using a Fisher test.
Obviously, choices must be made regarding which effects and interactions to include in the model.
\(Y\): quantitative response variable
\(X\): quantitative variable,
\(Z\): factor with \(I\) modalities \(\{A_1, \dots, A_I\}\)
lm(Y~X+Z) estimates the model without interaction, where, for each individual \(k = 1, \ldots, n\):
\[Y_k = m + \beta X_k + \sum_{i=2}^{I} \alpha_i \mathbf{1}\{Z_k=A_i\} + \varepsilon_k\]
Constraint \(\alpha_1 = 0\) is adopted to make the model identifiable.
lm(Y~X+A+X:A) or lm(Y~X*A) estimates the model with interaction:
\[ \begin{aligned} Y_k &= m + \beta X_k + \sum_{i=2}^{I} \beta_i X_k \mathbf{1}\{Z_k=A_i\} \\ &+ \sum_{i=2}^{I} \alpha_i \mathbf{1}\{Z_k=A_i\} + \varepsilon_k \end{aligned} \]
Key Insight: In the interaction model, the coefficient \(\beta\) associated with \(X\) varies according to the modality \(A_i\).
Daily food (g) of \(n = 200\) dogs: breed (lab / chihuahua), ster (sterilised yes / no), km walked per day, \(0\) to \(10\) for both breeds.
Estimate Std. Error t value Pr(>|t|)
(Intercept) 88.31 8.364 10.559 < 2e-16 ***
breedlab 486.40 6.703 72.571 < 2e-16 ***
steryes -59.29 6.766 -8.763 9.02e-16 ***
km 24.99 1.170 21.364 < 2e-16 ***
Residual standard error: 47.35 on 196 degrees of freedomOne slope for every dog: \(+25\) g per km, whatever the breed.
Estimate Std. Error t value Pr(>|t|)
(Intercept) 148.88 4.840 30.763 < 2e-16 ***
breedlab 348.67 6.544 53.280 < 2e-16 ***
km 12.16 0.783 15.528 < 2e-16 ***
steryes -58.92 3.364 -17.514 < 2e-16 ***
breedlab:km 28.26 1.156 24.455 < 2e-16 ***
Residual standard error: 23.54 on 195 degrees of freedomTwo slopes now: \(12.2\) g per km for a chihuahua, \(12.16 + 28.26 = 40.4\) for a labrador. Sterilisation stays a flat \(-59\) g, and \(\hat\sigma\) falls from \(47\) to \(24\).
Introduction to the Generalized Linear Model
⬅ Previous: Selection · Next: Introduction to the Generalized Linear Model ➡