1 Introduction

Regression is one of the fundamental topics in statistical modelling, providing a mechanism for expressing the potential dependence of some response of interest on a number of covariates. Its uses are widespread, permeating almost every scientific field, providing the bedrock for analysis of data arising from medicine, population health, demography, agriculture, sports and many more. In this module we will focus on the particular modelling case of when the response variable is continuous.

We will begin with a recap of the multiple linear regression model in chapter 1, extending some of the ideas you have seen in MAS2902. Subsequently, we will consider checking our core assumptions and using regression diagnostics to identify any unusual observations. The method of analysis of variance (anova) will then be introduced in chapter 3, allowing us to remove predictors from a model using the extra sum of squares technique. The fourth chapter will introduce indicator variables to handle binary explanatory variables and consider model selection techniques. Focus then shifts in chapter 5 to looking at models to handle factors in the context of designed experiments, giving rise to classic one-way and two-way anova models.

1.1 Multiple linear regression

As you have seen in MAS2902, there are often several possible explanatory variables to consider in a regression model. For example, in modelling someone’s weight we may want to include height and age (as well as a few other things) as covariates. We can extend the simple linear regression model to the multiple linear regression model \[ Y_i = \beta_0 + \beta_1 x_{i1} + \beta_2 x_{i2} + \ldots \beta_p x_{ip} + \epsilon_i \tag{1.1} \] for \(i = 1, \ldots n\). As in simple linear regression, we make assumptions about the (unobserved) error terms, namely that they have zero mean and are independently normally distributed, i.e. \(\epsilon_i \sim N(0, \sigma_{\epsilon}^2)\) and \(\textrm{Cov}(\epsilon_i, \epsilon_j) = 0\).

Note our model in (1.1) is called a multiple linear regression model because it is linear in the parameters, \(\underline{\beta} = \left(\beta_0, \ldots, \beta_p\right)^T\) - it does not have to be linear in the explanatory variables. For example, the following are multiple linear regression models \[\begin{eqnarray*} Y_i &=& \beta_0 + \beta_1 x_i + \beta_2 x_i^2 + \epsilon_i \\ Y_i &=& \beta_0 + \beta_1 x_{i1} + \beta_2 x_{i2}^2 + \beta_3 \ln x_{i3} + \epsilon_i \end{eqnarray*}\] for \(i = 1, \ldots n\). To see this, we can substitute \(w_i\), say, for \(x_i^2\) in the first model. However, the below example is not a multiple linear regression model: \[ Y_i = \beta_0 + e^{\beta_1 x_{i1} + \beta_2 x_{i2}} + \epsilon_i \;\; \textrm{for}\;\; i = 1, \ldots n. \] Sometimes seemingly nonlinear models can be transformed to linear models, but this will have a knock-on effect on the assumptions.

1.2 Matrix form of the model

Suppose we have \(n\) observations then equation (1.1) gives us the following system of equations: \[\begin{eqnarray*} Y_1 &=& \beta_0 + \beta_1 x_{11} + \beta_2 x_{12} + \ldots \beta_p x_{1p} + \epsilon_1 \\ Y_2 &=& \beta_0 + \beta_1 x_{21} + \beta_2 x_{22} + \ldots \beta_p x_{2p} + \epsilon_2 \\ &\vdots& \\ Y_n &=& \beta_0 + \beta_1 x_{n1} + \beta_2 x_{n2} + \ldots \beta_p x_{np} + \epsilon_n. \end{eqnarray*}\]

Note that the sample size is usually far greater than the number of covariates (\(n>>p\)), i.e. the problem is regular. This set-up is more conveniently represented in matrix format by placing

  1. \(Y_i\) values in a vector, \(\underline{Y}\);
  2. \(x_{ij}\) values in a (design) matrix \(\mathrm{X}\);
  3. \(\epsilon_i\) values in a vector \(\underline{\epsilon}\);
  4. \(\beta_j\) values in a vector \(\underline{\beta}\).

This leads to

\[\begin{align*} \begin{pmatrix} Y_1 \\ Y_2 \\ \vdots \\ Y_n \\ \end{pmatrix} = \begin{pmatrix} 1 & x_{11} & x_{12}& \ldots& x_{1p} \\ 1 & x_{21} & x_{22}& \ldots& x_{2p} \\ \vdots& \vdots& \vdots& \vdots& \vdots \\ 1 & x_{n1} & x_{n2}& \ldots& x_{np} \\ \end{pmatrix} \begin{pmatrix} \beta_0 \\ \beta_1 \\ \vdots \\ \beta_p \end{pmatrix} + \begin{pmatrix} \epsilon_1 \\ \epsilon_2 \\ \vdots \\ \epsilon_n \end{pmatrix}. \end{align*}\] This can be compactly written as \[\begin{equation} \underline{Y} = \mathrm{X}\underline{\beta} + \underline{\epsilon}. \end{equation}\] The associated dimensions are.
(i) \([n \times 1]\) for the response vector, \(\underline{Y}\),
(ii) \([n \times (p+1)]\) for the design matrix, \(\mathrm{X}\),
(iii) \([(p+1) \times 1]\) for the parameter vector, \(\underline{\beta}\),
(iv) \([n \times 1]\) for the error vector, \(\underline{\epsilon}\).
The column of ‘1’s in the design matrix is necessary to give the constant term for each observation, i.e. the intercept. We will discuss the role of the intercept in greater detail later. We can collect the covariate information for each item or individual into a vector alongside the ’1’ for the intercept: \[ \underline{x}_i = \begin{pmatrix} 1 \\ x_{i1} \\ x_{i2} \\ \vdots \\ x_{ip} \end{pmatrix}, \] for \(i = 1, \ldots, n\). This allows us to express the design matrix as \[ \mathrm{X} = \begin{pmatrix} \underline{x}_1^T \\ \underline{x}_2^T \\ \vdots \\ \underline{x}_n^T \end{pmatrix}. \]

It is sometimes convenient to set \(\underline{x}_i^* = (x_{i1}, x_{i2}, \ldots x_{ip})\) for \(i = 1, \ldots, n\). This is the vector of explanatory variables alone (i.e. no intercept term) for the \(i^{th}\) observation, which we will make use of later. We will now illustrate the multivariate regression model with a few examples.

Example - Matrix form for pre-diabetes data

We now introduce a specific example, which we shall return to throughout the chapter. In a study on body weight in 24 patients over 50, who have been diagnosed with pre-diabetes. Data were collected on their weight (in kilograms), food consumption (in calories), and a(n) (ordinal) measure of how much exercise each patient takes, on average per week: 0 = no exercise, 1 = some exercise, 2 = moderate exercise, 3 = heavy exercise. The food consumption was calculated by averaging over a week’s consumption and we will treat the ordinal covariate (exercise) as if it was continuous. The full data are given in Table 1.1

Table 1.1: Weight, average food consumption and exercise score for twenty-four pre-diabetes patients in body weight study.
Weight (kg) Consumption (cal) Exercise Weight (kg) Consumption (cal) Exercise
60.4 2680 3 84.4 3160 0
81.1 3280 1 93.0 3330 0
94.9 3890 2 61.3 2360 2
86.4 3170 0 74.9 3030 0
90.3 3390 1 94.3 3390 0
60.4 2670 2 61.8 2700 2
77.8 2770 0 78.1 3090 0
85.0 3330 1 74.8 3020 1
71.6 2710 0 59.0 2410 0
64.6 2600 1 69.2 2830 2
75.1 2880 0 67.1 2620 0
89.6 3430 0 82.4 2820 1

We can produce scatterplots of weight against consumption and exercise:

Scatterplot of weight against (average) food consumption (left) and exercise (right) for the pre-diabetes dataScatterplot of weight against (average) food consumption (left) and exercise (right) for the pre-diabetes data

Figure 1.1: Scatterplot of weight against (average) food consumption (left) and exercise (right) for the pre-diabetes data

Comments:

  1. We see that body weight is approximately linearly related to food consumption with a positive slope. There may be one unusual observation (top-right)? Could this be influential? Could others? We will consider these questions in more detail in the next chapter.

  2. There is some evidence that body weight declines with increasing exercise but the effect is possibly being masked by the variability in food consumption. Note that we observe vertical strips when dealing with ordinal variables so trends are perhaps harder to see.

We wish to fit the model: \[ \color{red}{Y_i = \beta_0 + \beta_1 x_{i1} + \beta_2 x_{i2} + \epsilon_i} \] for \(i = 1, \ldots, n\), where \(x_{i1}\) is food consumption and \(x_{i2}\) is exercise for the \(i^{th}\) patient. Equivalently, using our matrix formulation, we can express the model as \[ \color{red}{\underline{Y} = \mathrm{X}\underline{\beta} + \underline{\epsilon}} \] We can define the relevant quantities using the available data for the multiple linear regression model as follows: \[\begin{align*} \color{red}{ \underline{Y} = \begin{pmatrix} 60.4 \\ 81.1 \\ \vdots \\ 82.4 \\ \end{pmatrix}, \mathrm{X} = \begin{pmatrix} 1 & 2680& 3 \\ 1 & 3280& 1 \\ \vdots& \vdots& \vdots \\ 1 & 2820& 1 \\ \end{pmatrix}, \underline{\beta} = \begin{pmatrix} \beta_0 \\ \beta_1 \\ \beta_2 \\ \end{pmatrix}, \underline{\epsilon} = \begin{pmatrix} \epsilon_1 \\ \epsilon_2 \\ \vdots \\ \epsilon_{24} \\ \end{pmatrix}.} \end{align*}\]

Note that both \(\underline{Y}\) and \(\underline{\epsilon}\) are unchanged by the addition of variables since these are fixed data in the case of the former, and, as yet, unknown errors in the case of the latter - note that since there is a single error term, \(\epsilon_i\), for each observation, \(Y_i\), the errors will always form a vector of length \(n\).

The residuals (which estimate the errors) can also be collected in a vector of length \(n\), but we only know the values of the residuals after model fitting and their values will be different under each model. Both \(\mathrm{X}\) and \(\underline{\beta}\) do change, with the addition of a column and row respectively for each additional covariate.

1.3 Parameter estimation

Having set up the model, how do we obtain estimates for the unknown parameters? One approach is to use maximum likelihood.

In order to set the scene, recall that our key assumptions are normality, zero mean, common variance and independence of the errors, i.e. \(\epsilon_i \sim N(0, \sigma_{\epsilon}^2)\) and \(\textrm{Cov}(\epsilon_i, \epsilon_j) = 0, i \neq j\). Since \(Y_i = \beta_0 + \beta_1 x_{i1} + \ldots + \beta_p x_{ip} + \epsilon_i\), then the normality assumption on \(\epsilon_i\) induces a normal distribution on each \(Y_i\) (and, in turn, the vector \(\underline{Y}\)). The values for the covariates are fixed, as are the true (usually unknown) values of \(\underline{\beta}\) so these just serve as constants when thinking about the distribution of \(Y_i\).

Furthermore, the assumption of independence also carries through and the variance is unaltered (since we are simply adding scalars to a random variable). This tells us that \(Y_i \mid \underline{x}_i, \underline{\beta}, \sigma_{\epsilon}^2 \sim N(\beta_0 + \beta_1 x_{i1} + \ldots + \beta_p x_{ip}, \sigma_{\epsilon}^2)\) for \(i = 1, \ldots, n\). This can be written more compactly as \(Y_i \mid \underline{x}_i, \underline{\beta}, \sigma_{\epsilon}^2 \sim N(\underline{x}_i^T\underline{\beta}, \sigma_{\epsilon}^2)\).

Before going further, recall that the probability density function for a univariate normal random variable, \(Y \sim N(\mu, \sigma^2)\), is: \[ \color{red}{f(y \mid \mu, \sigma^2) = \frac{1}{\sqrt{2\pi\sigma^2}} \exp\left\{-\frac{(y - \mu)^2}{2\sigma^2}\right\}} \]

The likelihood principle instructs us to pick values of the parameters that maximise the likelihood. If observations are independent (as we have just shown they are here), then the likelihood function for all observations is a product of the individual normal densities for each observation of the form: \[ \color{red}{L(\mu, \sigma^2 \mid y_1, \ldots, y_n) = \prod_{i=1}^n f(y_i \mid \mu, \sigma^2).} \] Now, for our multiple linear regression model, we have \[\begin{align*} \color{red}{L(\underline{\beta}, \sigma_{\epsilon}^2 \mid y_1, \ldots, y_n, \mathrm{X})} &\color{red}{= \prod_{i=1}^n f(y_i \mid \underline{x}_i, \underline{\beta}, \sigma_{\epsilon}^2)} \\ &\color{red}{= \prod_{i=1}^n \frac{1}{\sqrt{2\pi\sigma_{\epsilon}^2}} \exp\left\{-\frac{(y_i - \underline{x}_i^T\underline{\beta})^2}{2\sigma_{\epsilon}^2}\right\}} \\ &\color{red}{= (2\pi\sigma_{\epsilon}^2)^{-n/2} \exp\left\{-\frac{1}{2\sigma_{\epsilon}^2}\sum_{i=1}^n (y_i - \underline{x}_i^T\underline{\beta})^2 \right\}} \\ &\color{red}{= (2\pi)^{-n/2} (\sigma_{\epsilon}^2)^{-n/2} \exp\left\{-\frac{1}{2\sigma_{\epsilon}^2}\sum_{i=1}^n (y_i - \underline{x}_i^T\underline{\beta})^2 \right\}} \\ \end{align*}\]

1.3.1 Estimation of \(\underline{\beta}\)

To estimate the parameter vector, \(\underline{\beta}\) we consider the log-likelihood since this is typically easier to work with, being additive as opposed to multiplicative. The log-likelihood is given by \[\begin{equation*} \color{red}{\ell(\underline{\beta}, \sigma_{\epsilon}^2 \mid y_1, \ldots, y_n, \mathrm{X}) = -\frac{n}{2} \ln(2\pi) -\frac{n}{2} \ln (\sigma_{\epsilon}^2) - \frac{1}{2\sigma_{\epsilon}^2} \sum_{i=1}^n (y_i - \underline{x}_i^T\underline{\beta})^2} \end{equation*}\] We could (partially) differentiate with respect to each element of \(\underline{\beta}\) and this would lead to a system of equations (with a lot of structure) known as the normal equations. However, it is easier to consider the matrix formulation of the model by noting that \[ \sum_{i=1}^n (y_i - \underline{x}_i^T\underline{\beta})^2 = (\underline{y} - \mathrm{X}\underline{\beta})^T (\underline{y} - \mathrm{X}\underline{\beta}) \]

This is, in fact, the exact quantity that is minimised using the method of least squares, which can also be used to estimate the parameters for this model (and simple linear regression). Now, ignoring the constant terms (by treating \(\sigma_{\epsilon}^2\) as fixed at this stage), we wish to minimise \[\begin{align*} \color{red}{(\underline{y} - \mathrm{X}\underline{\beta})^T (\underline{y} - \mathrm{X}\underline{\beta})} &\color{red}{= (\underline{y}^T - \underline{\beta}^T\mathrm{X}^T)(\underline{y} - \mathrm{X}\underline{\beta})} \\ &\color{red}{= \underline{y}^T\underline{y} - \underline{y}^T\mathrm{X}\underline{\beta} - \underline{\beta}^T\mathrm{X}^T\underline{y} + \underline{\beta}^T\mathrm{X}^T \mathrm{X}\underline{\beta}} \end{align*}\]

By noting that \(\underline{y}^T\mathrm{X}\underline{\beta}\) is a scalar we can rewrite this as its transpose, i.e. \(\underline{y}^T\mathrm{X}\underline{\beta} = \underline{\beta}^T\mathrm{X}^T\underline{y}\). Hence \[ \color{red}{(\underline{y} - \mathrm{X}\underline{\beta})^T (\underline{y} - \mathrm{X}\underline{\beta}) = \underline{y}^T\underline{y} - 2\underline{\beta}^T\mathrm{X}^T\underline{y} + \underline{\beta}^T\mathrm{X}^T \mathrm{X}\underline{\beta}} \]

We can now differentiate to obtain \[ \color{red}{\frac{\partial \ell}{\partial \underline{\beta}} = - 2\mathrm{X}^T\underline{y} + 2\mathrm{X}^T \mathrm{X}\underline{\beta}} \]

Setting equal to zero and solving for \(\underline{\beta}\) leads to the solution \[\begin{equation} \color{red}{\underline{\hat{\beta}} = (\mathrm{X}^T\mathrm{X})^{-1}\mathrm{X}^T\underline{y}} \end{equation}\]

This estimate exists as long as the inverse exists, i.e. no column of \(\mathrm{X}\) is a linear combination of other columns, i.e. there is no multicollinearity. We have to be careful that no columns of \(\mathrm{X}\) are linearly related as this can be harder to detect and leads to serious issues; we will discuss multicollinearity in more detail in section 1.7.

1.3.2 Estimation of \(\sigma_{\epsilon}^2\)

We can estimate \(\sigma_{\epsilon}^2\) in a similar fashion (now treating \(\underline{\beta}\) as fixed), using maximum likelihood once more. Recall that \[\begin{equation*} \ell(\underline{\beta}, \sigma_{\epsilon}^2 \mid y_1, \ldots, y_n, \mathrm{X}) = -\frac{n}{2} \ln(2\pi) -\frac{n}{2} \ln (\sigma_{\epsilon}^2) - \frac{1}{2\sigma_{\epsilon}^2} \sum_{i=1}^n (y_i - \underline{x}_i^T\underline{\beta})^2. \end{equation*}\]

For ease of calculation we let \(\tau = \sigma_{\epsilon}^2\) and then differentiating with respect to \(\tau\) we obtain \[\begin{equation*} \color{red}{\frac{\partial \ell}{\partial\tau} = -\frac{n}{2\tau} + \frac{\sum (y_i - \underline{x}_i^T\underline{\beta})^2}{2\tau^2}} \end{equation*}\]

Setting the above equal to zero and solving for \(\tau\) we obtain \[\begin{equation} \color{red}{\hat{\tau} = \frac{\sum (y_i - \underline{x}_i^T\underline{\beta})^2}{n} = \frac{\sum (y_i - \hat{y_i})^2}{n}.} \end{equation}\]

However, this is a biased estimate (akin to the sample variance bias problem), so we adjust for the fact that we have estimated the \(p\)-vector \(\underline{\beta}\) by using \[ \color{red}{\hat{\sigma}_{\epsilon}^2 = s^2 = \frac{\sum (y_i - \hat{y_i})^2}{n - p - 1} } \tag{1.2} \]

where \(p\) is the number of explanatory variables in the model. Note that \(\underline{\beta}\) has length \(k = p+1\) typically, with the additional intercept term.

1.3.3 Residuals, fitted values and the ‘hat matrix’

The vector of residuals (which estimate the errors) can be obtained by subtraction after model fitting, namely as ‘observed - fitted’ \[\begin{align*} \color{red}{\underline{\hat{\epsilon}}} &\color{red}{= \underline{y} - \underline{\hat{y}}} \\ &\color{red}{= \underline{y} - \mathrm{X}\underline{\hat{\beta}},} \end{align*}\] where the fitted values are found as \(\underline{\hat{y}} = \mathrm{X}\underline{\hat{\beta}}\).

The hat matrix

We can rewrite the estimate for the errors by substituting in \(\underline{\hat{\beta}} = (\mathrm{X}^T\mathrm{X})^{-1}\mathrm{X}^T\underline{y}\) to obtain \[\begin{align*} \color{red}{\underline{\hat{\epsilon}}} &\color{red}{= \underline{y} - \mathrm{X}(\mathrm{X}^T\mathrm{X})^{-1}\mathrm{X}^T\underline{y}} \\ &\color{red}{= (\mathrm{I} - \mathrm{X}(\mathrm{X}^T\mathrm{X})^{-1}\mathrm{X}^T) \underline{y}} \\ &\color{red}{= (\mathrm{I} - \mathrm{H})\underline{y},} \end{align*}\] where \(\mathrm{H} = \mathrm{X}(\mathrm{X}^T\mathrm{X})^{-1}\mathrm{X}^T\). \(\mathrm{H}\) is known as the ‘hat’ matrix since \[\begin{align*} \color{red}{\underline{\hat{y}}} &\color{red}{= \mathrm{X}\underline{\hat{\beta}}} \\ &\color{red}{= \mathrm{X}(\mathrm{X}^T\mathrm{X})^{-1}\mathrm{X}^T\underline{y}} \\ &\color{red}{= \mathrm{H}\underline{y}.} \end{align*}\]

Hence, multiplying by \(\mathrm{H}\) converts \(\underline{y}\) to \(\underline{\hat{y}}\), i.e. it is the matrix that puts a hat on \(\underline{y}\). The hat matrix is an \(n \times n\) matrix with elements \[\begin{align*} \color{red}{\mathrm{H} = \begin{pmatrix} h_{11}& h_{12}& \ldots& h_{1n} \\ h_{21}& h_{22}& \ldots& h_{2n} \\ \vdots& \vdots& & \vdots \\ h_{n1}& h_{n2}& \ldots& h_{nn} \\ \end{pmatrix}} \;\;\;\; \end{align*}\]

The diagonal values of \(\mathrm{H}\) (i.e. the \(h_{ii}\) values for \(i = 1, \ldots, n\)) are called the leverages (see chapter 2).

1.3.4 Properties of the hat matrix

It turns out that the hat matrix, \(\mathrm{H}\), has some useful properties, which will prove to be handy later. Namely,

  1. \(\mathrm{H}\) is symmetric, whereby \(\mathrm{H}^T = \mathrm{H}\),
  2. \(\mathrm{H}\) is idempotent, i.e \(\mathrm{H}^2 = \mathrm{H}\mathrm{H} = \mathrm{H}\).

Proof

  1. Now \[\begin{align*} \color{red}{\mathrm{H}^T} &\color{red}{=\Bigl(\mathrm{X} (\mathrm{X}^T\mathrm{X})^{-1} \mathrm{X}^T\Bigr)^T} \\ &\color{red}{= \mathrm{X}\left(\mathrm{X}^T\mathrm{X})^{-1}\right)^T \mathrm{X}^T} \\ &\color{red}{= \mathrm{X}\left(\mathrm{X}^T\mathrm{X})^T\right)^{-1} \mathrm{X}^T} \\ &\color{red}{= \mathrm{X} (\mathrm{X}^T\mathrm{X})^{-1} \mathrm{X}^T} \\ &\color{red}{= \mathrm{H}.} \\ \end{align*}\]
  2. We now have \[\begin{align*} \color{red}{\mathrm{H}\mathrm{H}} &\color{red}{= \Bigl(\mathrm{X} (\mathrm{X}^T\mathrm{X})^{-1} \mathrm{X}^T\Bigr) \Bigl(\mathrm{X} (\mathrm{X}^T\mathrm{X})^{-1} \mathrm{X}^T\Bigr)} \\ &\color{red}{= \mathrm{X} (\mathrm{X}^T\mathrm{X})^{-1} \mathrm{X}^T\mathrm{X} (\mathrm{X}^T\mathrm{X})^{-1} \mathrm{X}^T} \\ &\color{red}{= \mathrm{X} (\mathrm{X}^T\mathrm{X})^{-1} \mathrm{X}^T} \\ &\color{red}{= \mathrm{H}.} \\ \end{align*}\]

Example: Multiple linear regresion analysis of bodyweight data

We are now in a position to estimate the parameters for the data on pre-diabetes patients introduced earlier. First we recall that for the multiple linear regression model for the bodyweight data we have the design matrix given by: \[\begin{align*} \mathrm{X} = \begin{pmatrix} 1 & 2680& 3 \\ 1 & 3280& 1 \\ \vdots& \vdots& \vdots \\ 1 & 2820& 1 \\ \end{pmatrix} \end{align*}\] Using matrix algebra we can calculate \[\begin{align*} \mathrm{X}^T\mathrm{X} &= \begin{pmatrix} 24 & 71560 & 19 \\ 71560 & 216577400 & 55380\\ 19 & 55380 & 35\\ \end{pmatrix} \\ \mathrm{X}^T\underline{y} &= \begin{pmatrix} 1837.50 \\ 5570449 \\ 1354.60 \end{pmatrix}. \end{align*}\]

Taking the (3 by 3) matrix inverse we get

\[\begin{align*} (\mathrm{X}^T\mathrm{X})^{-1} = \begin{pmatrix} 3.02 & -9.69 \times 10^{-4} & -0.10 \\ -9.69 \times 10^{-4} & 3.20 \times 10^{-7} & 2.04 \times 10^{-5} \\ -0.10 & 2.04 \times 10^{-5} & 0.05 \\ \end{pmatrix} \end{align*}\]

Hence, the parameter estimates can be found as

\[\begin{eqnarray*} \color{red}{ \underline{\beta} = (\mathrm{X}^T\mathrm{X})^{-1}\mathrm{X}^T\underline{y} = (-2.104, 0.027, -3.278)^T.} \end{eqnarray*}\] Note that we can also calculate \(\underline{\hat{y}} = \mathrm{X}\underline{\hat{\beta}}\) from the above, and subsequently \(\underline{\hat{\varepsilon}} = \underline{y} - \underline{\hat{y}}\). The fitted line for the multiple linear regression model is \[ \color{red}{\textrm{Weight} = -2.014 + 0.027\times \textrm{Consumption} -3.278\times\textrm{Exercise}.} \] This can be interpreted in a similar way to simple linear regression, but with a few caveats:

  • body weight goes up by 0.027 kg for every additional calorie consumed, .

  • body weight decreases by around 3.3 kg as individuals move up an exercise category, .

  • be careful not to interpret the above as exercise being ‘more important’ than consumption due to having a larger coefficient - the scales of the variables are different and we also have no idea (yet!) whether these values are significant.

It may be more meaningful to express the change due to consumption in different units. Also, note that the values of the parameter estimates change with the introduction (or removal) of variables into (from) the model - this is always the case (unless the covariates are independent), no matter how significant (or not) they are. This is important when building a regression model.

We have seen that fitting a multiple linear regression model with two covariates can be achieved `by hand’. However, it is clear that as we look to build more complex models then it may be advantageous to use software - we will see how to do this in section 1.5.

1.4 Expectations, variances and inference

We now consider the properties of the estimators \(\underline{\hat{\beta}}\) and \(\underline{\hat{\epsilon}}\), i.e. is \(\underline{\hat{\beta}}\) unbiased? This will allow us, among other things, to assess the significance (or otherwise) of the parameter estimates. We begin by considering the expectation and variance of \(\underline{\hat{\beta}}\).

1.4.1 Expectation of \(\underline{\hat{\beta}}\)

Now, \[\begin{align*} \color{red}{\mathrm{E}\left[\underline{\hat{\beta}}\right]} &\color{red}{= \mathrm{E}\left[\left(\mathrm{X}^T\mathrm{X}\right)^{-1}\mathrm{X}^T\underline{Y}\right]} \\ &\color{red}{= \left(\mathrm{X}^T\mathrm{X}\right)^{-1}\mathrm{X}^T \mathrm{E}[\underline{Y}]} \\ &\color{red}{= \left(\mathrm{X}^T\mathrm{X}\right)^{-1}\mathrm{X}^T \mathrm{X}\underline{\beta}} \\ &\color{red}{= \underline{\beta}} \end{align*}\]

Hence, \(\underline{\hat{\beta}}\) is an unbiased estimator of \(\underline{\beta}\).

1.4.2 Variance of \(\underline{\hat{\beta}}\)

Before looking at the variance in detail we note \[ \left\{\left(\mathrm{X}^T\mathrm{X}\right)^{-1}\mathrm{X}^T\right\}^T = \mathrm{X}\left\{\left(\mathrm{X}^T\mathrm{X}\right)^{-1}\right\}^T = \mathrm{X}\left(\mathrm{X}^T\mathrm{X}\right)^{-1} \] since \((\mathrm{X}^T\mathrm{X})^{-1}\) is a symmetric matrix.

We are now in a position to look at the variance \[\begin{align*} \color{red}{\mathrm{Var}\left[\underline{\hat{\beta}}\right]} &\color{red}{= \mathrm{Var}\left[\left(\mathrm{X}^T\mathrm{X}\right)^{-1}\mathrm{X}^T\underline{Y}\right]} \\ &\color{red}{= \left(\mathrm{X}^T\mathrm{X}\right)^{-1}\mathrm{X}^T \mathrm{Var}\left[\underline{Y}\right] \left\{\left(\mathrm{X}^T\mathrm{X}\right)^{-1}\mathrm{X}^T\right\}^T} \\ &\color{red}{= \left(\mathrm{X}^T\mathrm{X}\right)^{-1}\mathrm{X}^T \mathrm{I}\sigma_{\epsilon}^2 \mathrm{X}\left(\mathrm{X}^T\mathrm{X}\right)^{-1}} \\ &\color{red}{= \sigma_{\epsilon}^2 \left(\mathrm{X}^T\mathrm{X}\right)^{-1}\mathrm{X}^T \mathrm{X}\left(\mathrm{X}^T\mathrm{X}\right)^{-1}} \\ &\color{red}{= \sigma_{\epsilon}^2 \left(\mathrm{X}^T\mathrm{X}\right)^{-1}} \end{align*}\]

1.4.3 Inference for \(\underline{\hat{\beta}}\)

Since \(\underline{\hat{\beta}}\) consists of linear combinations of the \(Y_i\)’s, which are independent and normally distributed, it has a multivariate normal distribution, namely \(\underline{\hat{\beta}} \sim N_{p+1}\left(\underline{\beta}, \sigma_{\epsilon}^2 (\mathrm{X}^T\mathrm{X})^{-1}\right)\) and each of the individual parameter estimates are univariate normal (due to properties of the multivariate normal distribution). Their (individual) significance can be asssessed via the test statistic \[ \color{red}{\hat{\beta_j}\bigg/\sqrt{v_{jj}s^2} \sim t_{n-p-1}} \] where \(v_{jj}\) is the \((j+1)^{th}\) diagonal element of \(\mathrm{V} = (\mathrm{X}^T\mathrm{X})^{-1}\), and \(s^2\) is our (unbiased) estimate of \(\sigma_{\epsilon}^2\) from equation (1.2). Note the use of the \(t\)-distribution since we must also estimate \(s^2\).

1.4.4 Expectation and variance of the fitted values

The fitted values are calculated as \[ \underline{\hat{Y}} = \mathrm{X}\underline{\hat{\beta}} \] or equivalently as \[ \underline{\hat{Y}} = \mathrm{H}\underline{Y}. \]

Their expectation is

\[\begin{align*} \color{red}{\mathrm{E}\left[\underline{\hat{Y}}\right]} &\color{red}{= \mathrm{E}\left[\mathrm{X}\underline{\hat{\beta}}\right]} \\ &\color{red}{=\mathrm{X} \mathrm{E}[\underline{\hat{\beta}}]} \\ &\color{red}{=\mathrm{X}\underline{\beta}}, \end{align*}\]

with variance given by

\[\begin{align*} \color{red}{\mathrm{Var}\left[\underline{\hat{Y}}\right]} &\color{red}{= \mathrm{Var}\left[\mathrm{H}\underline{Y}\right]} \\ &\color{red}{= \mathrm{H} \mathrm{Var}\left[\underline{Y}\right] \mathrm{H}^T} \\ &\color{red}{= \mathrm{H} \mathrm{I}\sigma_{\epsilon}^2 \mathrm{H}^T} \\ &\color{red}{= \mathrm{H}\mathrm{H}\sigma_{\epsilon}^2} \\ &\color{red}{= \mathrm{H} \sigma_{\epsilon}^2} \end{align*}\] Hence, the variability of the fitted values depends on the hat matrix, \(\mathrm{H}\). We will discuss this further in chapter 2.

1.4.5 Expectation and variance of the residuals

Recall that the residuals are found as \[ \color{red}{\underline{\hat{\epsilon}} = \underline{Y} - \underline{\hat{Y}}.} \] or alternatively as \[ \color{red}{\underline{\hat{\epsilon}} = (\mathrm{I} - \mathrm{H})\underline{Y}.} \]

We can find the expectation and variance as \[\begin{align*} \color{red}{\mathrm{E}[\underline{\hat{\epsilon}}]} & \color{red}{= \mathrm{E}[\underline{Y} - \underline{\hat{Y}}]} \\ &\color{red}{= \mathrm{E}[\mathrm{X}\underline{\beta}] - \mathrm{E}\left[\mathrm{X}\underline{\hat{\beta}}\right]} \\ &\color{red}{= \mathrm{X}\underline{\beta} - \mathrm{X}\underline{\beta}} \\ &\color{red}{= \underline{0},} \end{align*}\] and \[\begin{align*} \color{red}{\mathrm{Var}[\underline{\hat{\epsilon}}]} &\color{red}{= \mathrm{Var}[(\mathrm{I} - \mathrm{H})\underline{Y}]} \\ &\color{red}{= (\mathrm{I} - \mathrm{H})\mathrm{Var}[\underline{Y}](\mathrm{I} - \mathrm{H})^T} \\ &\color{red}{= (\mathrm{I} - \mathrm{H})\sigma_{\epsilon}^2\mathrm{I}(\mathrm{I} - \mathrm{H})} \\ &\color{red}{= \sigma_{\epsilon}^2(\mathrm{I} - \mathrm{H} - \mathrm{H} + \mathrm{H}\mathrm{H})} \\ &\color{red}{= \sigma_{\epsilon}^2(\mathrm{I} - \mathrm{H}).} \end{align*}\]

Note that this implies that, unless all the diagonal values of \(\mathrm{H}\) are equal then the errors have different variances, and that these variances are smaller for larger values of \(h_{ii}\), i.e. higher leverages (see chapter 2).

1.5 Multiple linear regression in R

Once we start to think about large datasets and a large number of parameters, finding the parameter estimates by hand becomes laborious, not to mention the possibility of both data entry and/or numerical errors occurring increases greatly. Happily, we can use R to conduct the analyses instead.

1.5.1 Using data in R

There are various ways of using data with R. Data can be read in manually, i.e.

## Type the data in the console
bodyweight = c(60.4, 81.1, 94.9, 86.4, 90.3, 60.4, 77.8, 85.0, 71.6, 64.6, 75.1, 89.6, 
84.4, 93.0, 61.3, 74.9, 94.3, 61.8, 78.1, 74.8, 59.0, 69.2, 67.1, 82.4)

The majority of the time, in this module and the wider world, the (external) data in the file ExternalData.RData will be read/loaded directly into R, e.g.

## Load in an external dataset
load("ExternalData.RData")

Alternatively, we may sometimes make use of datasets that are internal to R in that they are part of an R package, i.e. for the dataset InternalRDataset:

## Load in an internal dataset
data(InternalRDataset)

To view the available datasets in R we can type data() at the console, or, for datasets attached to a particular package we can use data(library = "Rpackage").

Example: Analysis of bodyweight data using R

To use R for the plots and analysis seen earlier:

## Load the data
load("bodyweight.RData")
## Plots
# Weight versus consumption
plot(Weight ~ Consumption, data = bodyweight, pch = 16)
# Weight versus exercise
plot(Weight ~ Exercise, data = bodyweight, pch = 16)
## Analysis
# Simple linear regression on consumption
fit1 <- lm(Weight ~ Consumption, data = bodyweight)
# Multiple linear regression on consumption & exercise
fit2 <- lm(Weight ~ Consumption + Exercise, data = bodyweight)
We can inspect a model fit using various commands
  1. The summary() command gives an overview of the fit 
    2. summary(fit2)
    ## 
## Call:
## lm(formula = Weight ~ Consumption + Exercise, data = bodyweight)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -5.5751 -2.5704 -0.7894  2.4049 10.9266 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -2.104925   7.017861  -0.300  0.76717    
## Consumption  0.027254   0.002286  11.921 8.23e-11 ***
## Exercise    -3.278296   0.916795  -3.576  0.00178 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.044 on 21 degrees of freedom
## Multiple R-squared:  0.8916, Adjusted R-squared:  0.8813 
## F-statistic:  86.4 on 2 and 21 DF,  p-value: 7.346e-11
    We will consider output of this nature in detail later in the module.

  2. The fitted values and residuals can also be extracted (output - to three decimal places - is suppressed here)

    round(fitted.values(fit2), 3)
round(residuals(fit2), 3)

  3. The variance-covariance matrix for \(\underline{\hat{\beta}}\) is also contained within the fit. For the second model fit we get

    vcov(fit2)
    ##             (Intercept)   Consumption      Exercise
## (Intercept)  49.2503765 -1.584890e-02 -1.6584367474
## Consumption  -0.0158489  5.227021e-06  0.0003330453
## Exercise     -1.6584367  3.330453e-04  0.8405137595

  4. The hat-values that make up the diagonal of the \(\mathrm{H}\) matrix can also be found - again we round to three decimal places:

    round(hatvalues(fit2), 3)
    ##     1     2     3     4     5     6     7     8     9    10    11    12    13 
## 0.294 0.075 0.425 0.079 0.101 0.132 0.095 0.086 0.106 0.087 0.080 0.124 0.078 
##    14    15    16    17    18    19    20    21    22    23    24 
## 0.101 0.210 0.073 0.114 0.128 0.074 0.045 0.197 0.117 0.127 0.051

  5. We can also add a fitted regression line to a scatterplot via the abline() command:

    plot(Weight ~ Consumption, data = bodyweight, pch = 16)
abline(fit1, lty = 2)

1.6 The role of the intercept

The intercept, via the parameter \(\beta_0\), is included as a matter of course when fitting a regression model (the default behaviour in R is to have an intercept present in a model). Why is this the case? What would happen if we removed the intercept?

Suppose we thought that we should fit the model without an intercept, then the multiple linear regression model takes the form \[ \color{red}{Y_i = \beta_0 + \beta_1 x_{i1} + \beta_2 x_{i2} + \ldots \beta_p x_{ip} + \epsilon_i \tag{1.3}} \] for \(i = 1, \ldots, n\). Or, equivalently, in matrix notation \[\begin{equation} \color{red}{\underline{Y} = \mathrm{\tilde{X}}\underline{\beta} + \underline{\epsilon}.} \end{equation}\] where \(\mathrm{\tilde{X}}\) represents the design matrix that does not now have a first column of ’1’s.

Example: Analysis of bodyweight data without an intercept term

Returning to our example on pre-diabetes we would have \[\begin{align*} \color{red}{\underline{Y} = \begin{pmatrix} 60.4 \\ 81.1 \\ \vdots \\ 82.4 \\ \end{pmatrix}, \mathrm{\tilde{X}} = \begin{pmatrix} 2680& 3 \\ 3280& 1 \\ \vdots& \vdots \\ 2820& 1 \\ \end{pmatrix}, \underline{\beta} = \begin{pmatrix} \beta_1 \\ \beta_2 \\ \end{pmatrix}, \underline{\epsilon} = \begin{pmatrix} \epsilon_1 \\ \epsilon_2 \\ \vdots \\ \epsilon_{24} \\ \end{pmatrix}.} \end{align*}\]

Note that \(\underline{Y}\) and \(\underline{\epsilon}\) are unchanged, whereas both the design matrix and \(\underline{\beta}\) are affected by the removal of the intercept term.

Upon fitting we would obtain the fitted model \[ \textrm{Weight} = \hat{\beta}_1\times \textrm{Consumption} + \hat{\beta}_2\times\textrm{Exercise}. \]

This model - and the equivalent model with an intercept term - assumes the relationship between weight and consumption remains the same for all values of calorific consumption and exercise. Moreover, the model without the intercept further assumes that zero calorie intake and zero exercise gives zero body weight!

This may well not be true (or possible), not just here but for many datasets. Forcing a zero intercept can give nonsensical values for the predicted response and it can also severely affect the fit of the regression line, particularly if our estimate of the intercept is significantly different from zero. In the absence of an intercept term, the line of best fit is forced to go through the origin. We will now investigate further with another example.

Example: Analysis of men’s Premier League football data - the role of the intercept

The data in the following example comes from the 2012-13 men’s English Premier League final football table (on Canvas in the file prem.RData). For each team the number of points they achieved (the response - why?), goals they scored, conceded, and their goal difference (scored - conceded) are recorded, alongside how many times they did not concede a goal (a ‘clean sheet’), which will be our primary focus for now. A snapshot of the data are given below:

load("prem.RData")
kable(head(prem, 5))
Position Team Scored Conceded GoalDifference Points CleanSheets
1 Manchester United 86 43 43 89 13
2 Manchester City 66 34 32 78 18
3 Chelsea 75 39 36 75 14
4 Arsenal 72 37 35 73 14
5 Tottenham Hotspur 66 46 20 72 9
Below is a scatterplot of points against clean sheets:
Scatterplot of points against clean sheets for the Premier League 2012/13 data.

Figure 1.2: Scatterplot of points against clean sheets for the Premier League 2012/13 data.

  1. Fit a simple linear regression model with clean sheets as the sole covariate. Overlay the regression line on the scatterplot of the raw data and comment.

    fitprem1 <- lm(Points ~ CleanSheets, data = prem)
plot(Points ~ CleanSheets, data = prem, pch = 16, xlim = c(0, 25), ylim = c(0, 100))
abline(fitprem1, lty = 2, lwd = 1.5, col = "red")
    Scatterplot of points against clean sheets for the Premier League 2012/13 data with overlaid model fit.

    Figure 1.3: Scatterplot of points against clean sheets for the Premier League 2012/13 data with overlaid model fit.

  2. Fit a second model, this time without an intercept and overlay this regression line. What do you observe?
    We can fit the second model and overlay the line using:

    fitprem2 <- lm(Points ~ CleanSheets - 1, data = prem)
abline(fitprem2, lty = 3, lwd = 1.5)
    Scatterplot of points against clean sheets for the Premier League 2012/13 data with two overlaid model fits.

    Figure 1.4: Scatterplot of points against clean sheets for the Premier League 2012/13 data with two overlaid model fits.

We can see that the model without the intercept has a different slope since \(\hat{\beta}_0\) and \(\hat{\beta}_1\) are correlated (see MAS2902). Note also that the line of best fit for the model without the intercept is forced to go through the origin. Both models seem to do a reasonable job - assessing by eye - of capturing the relationship between points and clean sheets. This is not always the case though, as we will see in practical 1.

1.6.1 Interpretability of the intercept and extrapolation

Having established that including an intercept is a sensible thing to do, we now move on to the question of its interpretation. Note that we did not formally interpret the intercept in our previous analysis of the bodyweight data, and this is common practice.

However, if we did wish to say something meaningful about the intercept how would we go about it? We first inspect the fit for our first model from the previous example:

summary(fitprem1)
## 
## Call:
## lm(formula = Points ~ CleanSheets, data = prem)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -16.6526  -7.9816  -0.7842   7.0974  26.8211 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  16.3368     8.0303   2.034 0.056915 .  
## CleanSheets   3.5263     0.7544   4.674 0.000189 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 12.3 on 18 degrees of freedom
## Multiple R-squared:  0.5483, Adjusted R-squared:  0.5232 
## F-statistic: 21.85 on 1 and 18 DF,  p-value: 0.0001888
Scatterplot of points against clean sheets for the Premier League 2012/13 data with overlaid model fit.

Figure 1.5: Scatterplot of points against clean sheets for the Premier League 2012/13 data with overlaid model fit.

We see that the estimate of the intercept, \(\hat{\beta}_0\), is 16.34. This tells us that when clean sheets takes the value zero, then we would expect a team to obtain around 16 points (as points is an integer we round). Here this makes some sense, since no clean sheets would mean a team concedes at least one goal in every match they play. Note, however, that the smallest observed value for this variable is five, so by using the value of zero we are extrapolating beyond the observed range of our data and this can be problematic.

1.6.2 Mean-centering of covariates

Without any data manipulation prior to model fitting we have seen that the estimate for the intercept is interpreted as the value for the response when all of the covariates take the value zero. This, however, might be a scenario that is either not likely (i.e. a weight of zero kg for an adult), or not permissible (amount of a drug administered as part of a treatment) in the context of the data at hand.

One solution to this issue is to scale the covariates via mean-centering. \[ \color{red}{\underline{\tilde{x}}^{(j)} = \underline{x}^{(j)} - \bar{x}^{(j)}} \] where \(j = 1, \ldots, p\), \(\underline{x}^{(j)} = (x_{1j}, x_{2j}, \ldots, x_{nj})\) is the vector of values for the \(j^{th}\) covariate and \(\bar{x}^{(j)}\) is the sample mean for the \(j^{th}\) covariate, for example the mean of the exercise values in the bodyweight data. Note the distinction between \(\underline{x}^{(j)}\) and \(\underline{x}_i\) introduced earlier, which is the vector of values for each individual (or subject).

The intercept has the same interpretation as above, namely the value of the response when the covariates are all simultaneously set to zero, i.e. \(\underline{x}_i = \underline{0}\). However, zero is now the mean value for each covariate, after mean-centering, so the intercept can also now be interpreted as the value of the response when each covariate is at its (own) average value. Furthermore, the value of the intercept turns out to be \(\bar{y}\), the sample mean of the response vector. Recall that in simple linear regression \[ \hat{\beta}_0 = \bar{y} - \hat{\beta}_1 \times \bar{x} \] and this will clearly reduce to \(\hat{\beta}_0 = \bar{y}\) when \(\bar{x} = 0\). This result generalises to the multiple linear regression case.

This tends to give a more intuitive interpretation generally. Mean-centering also removes the correlation between \(\beta_0\) and \(\beta_1, \ldots \beta_p\). We will now see the effect of mean-centering in an example.

Example: Mean-centering (men’s Premier League football data)

Returning to the data from the men’s football Premier League. Below is a scatterplot of points against the raw (solid circles) and mean-centered (triangles) versions of our clean sheets covariate.

Scatterplot of points against the observed clean sheets (solid circles) and their mean-centered counterpart (triangles).

Figure 1.6: Scatterplot of points against the observed clean sheets (solid circles) and their mean-centered counterpart (triangles).

  1. Fit a model using a mean-centered version of clean sheets.
    [Hint: use the scale command in R to perform the mean-centering].
    We fit - and inspect - the model using the R commands

    CleanSheetsScaled <- scale(prem$CleanSheets, 
                           scale = FALSE)
fit_mean_centre <- lm(Points ~ CleanSheetsScaled, 
                      data = prem)
summary(fit_mean_centre)  
    ## 
## Call:
## lm(formula = Points ~ CleanSheetsScaled, data = prem)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -16.6526  -7.9816  -0.7842   7.0974  26.8211 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)        51.6000     2.7513  18.755 2.91e-13 ***
## CleanSheetsScaled   3.5263     0.7544   4.674 0.000189 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 12.3 on 18 degrees of freedom
## Multiple R-squared:  0.5483, Adjusted R-squared:  0.5232 
## F-statistic: 21.85 on 1 and 18 DF,  p-value: 0.0001888
  2. Overlay the lines of best fit for the models using the raw and mean-centered covariates. What do you notice?
    Scatterplot of points against the observed clean sheets (solid circles) and their mean-centered counterpart (triangles) with overlaid lines of best fit.

    Figure 1.7: Scatterplot of points against the observed clean sheets (solid circles) and their mean-centered counterpart (triangles) with overlaid lines of best fit.

    3. From the summary, we see that the estimate of the slope is exactly the same as before, i.e. \(\hat{\beta}_1 = 3.53\) so the line has the same gradient, but the intercept is different. The intercept estimate is \(\hat{\beta}_0 = 51.60\) (recall, it was around 16 earlier) which suggests that a team with the number of clean sheets will obtain around 52 points (nearest integer, as before). This interpretation is cleaner than our earlier interpretation using the raw rather than mean-centered covariate.

Although we have illustrated the role of the intercept using simple linear regression, the same ideas hold in the multiple linear regression model. We now return to the issue of multicollinearity.

1.7 Properties of \(\left(\mathrm{X}^T\mathrm{X}\right)^{-1}\): multicollinearity

We saw earlier that both the estimator of \(\underline{\beta}\) and its variance depend on the quantity \(\left(\mathrm{X}^T\mathrm{X}\right)^{-1}\). As such, this quantity plays a critical part in fitting a regression model and in determining the significance (or otherwise) of estimated parameters. We will now consider a situation known as multicollinearity that leads to problems with taking the inverse of \(\mathrm{X}^T\mathrm{X}\).

Example: Multicollinearity in men’s Premier League football data

Returning to the Premier League football data, a sports data analyst sets out to fit the following model: \[ \textrm{Points}_i = \beta_0 + \beta_1 \textrm{Goal difference}_i + \beta_2 \textrm{Scored}_i + \beta_3 \textrm{Conceded}_i + \epsilon_i \]

  1. Construct the design matrix \(\mathrm{X}\), and hence calculate \(\mathrm{X}^T\mathrm{X}\) and \((\mathrm{X}^T\mathrm{X})^{-1}\).

    X <- cbind(1, prem$GoalDifference, prem$Scored, 
           prem$Conceded)
XTX <- crossprod(X)
solve(XTX)
    The last line fails, we cannot invert the matrix as it is singular.

  2. Can you spot an obvious problem with this model?

    The problem here is that one of the variables is a linear combination of the others, namely goal difference which is defined as ‘scored’ - ‘conceded’. This means they are collinear and that we will have problems inverting \(\mathrm{X}^T\mathrm{X}\).

  3. Fit the model in R and inspect the fit - what do you notice?

    Implementing the model we get

    summary(lm(Points ~ GoalDifference + Scored + Conceded, data = prem))
    ## 
## Call:
## lm(formula = Points ~ GoalDifference + Scored + Conceded, data = prem)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -11.0387  -2.7741   0.2508   3.7633   5.7760 
## 
## Coefficients: (1 not defined because of singularities)
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     38.3511     9.6368   3.980 0.000969 ***
## GoalDifference   0.5710     0.1097   5.206 7.13e-05 ***
## Scored           0.2493     0.1803   1.382 0.184780    
## Conceded             NA         NA      NA       NA    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.471 on 17 degrees of freedom
## Multiple R-squared:  0.9437, Adjusted R-squared:  0.937 
## F-statistic: 142.4 on 2 and 17 DF,  p-value: 2.404e-11
    There is no parameter estimate or standard error for goals conceded, which seems to have been removed from the model.

Here it was clear what was driving the multicollinearity, and the issue could be easily spotted, and resolved. Sometimes, however, the problem is more subtle and we will consider this scenario further in practical 1. In the next chapter we will investigate whether our model conforms to assumptions and/or has any unusual observations that warrant further investigation.