View the latest recording of this lecture
We have mostly focused on (i) models containing just numerical covariates (linear regression models) and (ii) models containing just factors (‘ANOVA type’ models).
Naturally there will be situations in which we wish to model a response in terms of both numerical covariates and factors. We saw this in the case of ANCOVA, e.g. the English Exam data and the Retail Space Data.
Generalizing factorial and regression models to include models with potentially both factors and numerical covariates produces the class of models (historically) referred to as General Linear Models.
The abbreviation GLM was used for these models for many years, and still is in some software. R uses this acronym to refer to generalised linear models which are covered in another course. As a consequence, we often find some confusion, but all linear models fit into the wider group of generalised linear models
We will look at general linear models in this lecture.
As we saw when we looked at factorial models, any factor can be coded using dummy variables so as to produce a linear regression model.
It follows that any general linear model can be fitted as a linear regression.
All the ideas about least squares estimation, fitted values and residuals remain unchanged.
In order to understand how to use and interpret models with a combination of factors and numerical covariates, consider the simple case where we have a response variable, y, a covariate x and a factor A.
There are a number of different models that we could fit.
\[Y_{ij} = \mu + \alpha_i + \beta_i x_{ij} + \varepsilon_{ij}\]
Note: intercept and slope depends on factor level i.
Treatment constraint fixes \(\alpha_1 =
0\).
R formula is Y ~ A + A:x
Here the notation of an interaction between the covariate and the factor is telling R to allow the regression slope to depend on the factor level.
\[Y_{ij} = \mu + \alpha_i + \beta x_{ij} + \varepsilon_{ij}\]
Note: constant slope, but intercept depends on factor level
i.
Treatment constraint sets \(\alpha_1 =
0\).
R formula is Y ~ A + x
\[Y_{ij} = \mu + \beta x_{ij} + \varepsilon_{ij}\]
Note: factor A plays no role here.
R formula is the familiar Y ~ x
Technical aside:
Another possible model would be different regression slopes but a common intercept, but this is not commonly used.
The ideas just presented can be extended in a natural fashion to models with two or more factors and covariates.
Suppose that the R formula is Y ~ A + A:B + A:x + z
,
where A
and B
are factors, and x
and z
are covariates. Then the model is defined by:
\[Y_{ijk} = \mu +\alpha_i + \beta_{ij} + \gamma_i x_{ijk} + \delta z_{ijk} + \varepsilon_{ijk}.\]
Suppose that the R formula is Y ~ A*B + A:x
This model is defined by \(Y_{ijk} = \mu +\alpha_i + \beta_j + (\alpha\beta)_{ij} + \gamma_i x_{ijk} + \varepsilon_{ijk}.\)
Suppose that the R formula is Y ~ A + A:x + B + B:z
This model is defined by \(Y_{ijk} = \mu +\alpha_i + \beta_j + \gamma_{i} x_{ijk} + \delta_j z_{ijk} + \varepsilon_{ijk}.\)
Samara are small winged fruit on maple trees. In Autumn these fruit fall to the ground, spinning as they go. Research on the aerodynamics of the fruit has applications for helicopter design.
In one study the following variables were measured on individual samara:
Velocity
: speed of fallTree
: data collected from 3 treesLoad
: ‘disk loading’ (an aerodynamical quantity based
on each fruit’s size and weight).The aim of this Case Study is to model Velocity
in terms
of Load
(a numerical covariate) and Tree
(a
factor).
'data.frame': 35 obs. of 3 variables:
$ Tree : int 1 1 1 1 1 1 1 1 1 1 ...
$ Load : num 0.239 0.208 0.223 0.224 0.246 0.213 0.198 0.219 0.241 0.21 ...
$ Velocity: num 1.34 1.06 1.14 1.13 1.35 1.23 1.23 1.15 1.25 1.24 ...
N.B. It might prove useful to have the values of Tree
as
both a numeric and a factor variable.
Samara |>
ggplot(mapping = aes(y = Velocity, x = Load)) + geom_point(mapping = aes(col = TreeF,
pch = TreeF))
We will fit a model with separate regression lines at different levels of TreeF. The following shows the effect parameterising the model two ways:
The first way gives estimates for the intercept and slope for each line separately.
Call:
lm(formula = Velocity ~ TreeF + TreeF:Load, data = Samara)
Residuals:
Min 1Q Median 3Q Max
-0.120023 -0.049465 -0.001298 0.049938 0.145571
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.5414 0.2632 2.057 0.0488 *
TreeF2 -0.8408 0.3356 -2.505 0.0181 *
TreeF3 -0.2987 0.4454 -0.671 0.5078
TreeF1:Load 3.0629 1.1599 2.641 0.0132 *
TreeF2:Load 6.7971 0.9511 7.147 7.26e-08 ***
TreeF3:Load 3.8834 1.9672 1.974 0.0580 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.07554 on 29 degrees of freedom
Multiple R-squared: 0.8436, Adjusted R-squared: 0.8167
F-statistic: 31.29 on 5 and 29 DF, p-value: 7.656e-11
Note the first model has TreeF1:Load, and the coefficients for TreeF2:Load and TreeF3.Load are tested against the hypotheses that the slope is zero for each of those trees. This may not be a sensible hypothesis.
The second way specifies a baseline (for TreeF1), and then adjustments to the intercept and slope for the second and third line. The second way gives us insight into where there are significant differences to the baseline, (if any).
Call:
lm(formula = Velocity ~ TreeF * Load, data = Samara)
Residuals:
Min 1Q Median 3Q Max
-0.120023 -0.049465 -0.001298 0.049938 0.145571
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.5414 0.2632 2.057 0.0488 *
TreeF2 -0.8408 0.3356 -2.505 0.0181 *
TreeF3 -0.2987 0.4454 -0.671 0.5078
Load 3.0629 1.1599 2.641 0.0132 *
TreeF2:Load 3.7343 1.5000 2.490 0.0188 *
TreeF3:Load 0.8205 2.2837 0.359 0.7220
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.07554 on 29 degrees of freedom
Multiple R-squared: 0.8436, Adjusted R-squared: 0.8167
F-statistic: 31.29 on 5 and 29 DF, p-value: 7.656e-11
As usual, the treatment constraint is used by default for
TreeF
(the factor version of Tree
). Hence the
Intercept
listed under the coefficients is the intercept
for TreeF1
, with adjustments for treeF2 and 3 given in the
summary table. The relevant hypothesis tests are for whether the
intercepts and slopes are the same or different to that of the
baseline.
Either way we get the fitted model is as follows:
\[\begin{aligned} &\mbox{Tree 1}&~~~E[ \mbox{Velocity}] = 0.5414 + 3.0629 \mbox{Load}\\ &\mbox{Tree 2}&~~~E[ \mbox{Velocity}] = -0.2994 + 6.7971 \mbox{Load}\\ &\mbox{Tree 3}&~~~E[ \mbox{Velocity}] = 0.2427 + 3.8834 \mbox{Load}\end{aligned}\]
Check that you can get the same numbers with either parameterisation.
Write down the design matrix for the model
samara.lm
.
Verify your answer using the model.matrix()
command.
(Intercept) TreeF2 TreeF3 TreeF1:Load TreeF2:Load TreeF3:Load
1 1 0 0 0.239 0.000 0.000
2 1 0 0 0.208 0.000 0.000
3 1 0 0 0.223 0.000 0.000
4 1 0 0 0.224 0.000 0.000
5 1 0 0 0.246 0.000 0.000
6 1 0 0 0.213 0.000 0.000
7 1 0 0 0.198 0.000 0.000
8 1 0 0 0.219 0.000 0.000
9 1 0 0 0.241 0.000 0.000
10 1 0 0 0.210 0.000 0.000
11 1 0 0 0.224 0.000 0.000
12 1 0 0 0.269 0.000 0.000
13 1 1 0 0.000 0.238 0.000
14 1 1 0 0.000 0.206 0.000
15 1 1 0 0.000 0.172 0.000
16 1 1 0 0.000 0.235 0.000
17 1 1 0 0.000 0.247 0.000
18 1 1 0 0.000 0.239 0.000
19 1 1 0 0.000 0.233 0.000
20 1 1 0 0.000 0.234 0.000
21 1 1 0 0.000 0.189 0.000
22 1 1 0 0.000 0.192 0.000
23 1 1 0 0.000 0.209 0.000
24 1 0 1 0.000 0.000 0.192
25 1 0 1 0.000 0.000 0.200
26 1 0 1 0.000 0.000 0.175
27 1 0 1 0.000 0.000 0.187
28 1 0 1 0.000 0.000 0.181
29 1 0 1 0.000 0.000 0.195
30 1 0 1 0.000 0.000 0.155
31 1 0 1 0.000 0.000 0.179
32 1 0 1 0.000 0.000 0.184
33 1 0 1 0.000 0.000 0.177
34 1 0 1 0.000 0.000 0.177
35 1 0 1 0.000 0.000 0.186
attr(,"assign")
[1] 0 1 1 2 2 2
attr(,"contrasts")
attr(,"contrasts")$TreeF
[1] "contr.treatment"
(Intercept) TreeF2 TreeF3 Load TreeF2:Load TreeF3:Load
1 1 0 0 0.239 0.000 0.000
2 1 0 0 0.208 0.000 0.000
3 1 0 0 0.223 0.000 0.000
4 1 0 0 0.224 0.000 0.000
5 1 0 0 0.246 0.000 0.000
6 1 0 0 0.213 0.000 0.000
7 1 0 0 0.198 0.000 0.000
8 1 0 0 0.219 0.000 0.000
9 1 0 0 0.241 0.000 0.000
10 1 0 0 0.210 0.000 0.000
11 1 0 0 0.224 0.000 0.000
12 1 0 0 0.269 0.000 0.000
13 1 1 0 0.238 0.238 0.000
14 1 1 0 0.206 0.206 0.000
15 1 1 0 0.172 0.172 0.000
16 1 1 0 0.235 0.235 0.000
17 1 1 0 0.247 0.247 0.000
18 1 1 0 0.239 0.239 0.000
19 1 1 0 0.233 0.233 0.000
20 1 1 0 0.234 0.234 0.000
21 1 1 0 0.189 0.189 0.000
22 1 1 0 0.192 0.192 0.000
23 1 1 0 0.209 0.209 0.000
24 1 0 1 0.192 0.000 0.192
25 1 0 1 0.200 0.000 0.200
26 1 0 1 0.175 0.000 0.175
27 1 0 1 0.187 0.000 0.187
28 1 0 1 0.181 0.000 0.181
29 1 0 1 0.195 0.000 0.195
30 1 0 1 0.155 0.000 0.155
31 1 0 1 0.179 0.000 0.179
32 1 0 1 0.184 0.000 0.184
33 1 0 1 0.177 0.000 0.177
34 1 0 1 0.177 0.000 0.177
35 1 0 1 0.186 0.000 0.186
attr(,"assign")
[1] 0 1 1 2 3 3
attr(,"contrasts")
attr(,"contrasts")$TreeF
[1] "contr.treatment"
Note the difference in the column for TreeF1:Load