---
subtitle: "TMA4315 Generalized linear models H2018"
title: "Module 2: MULTIPLE LINEAR REGRESSION Week 1"
author: "Mette Langaas, Department of Mathematical Sciences, NTNU -- with contributions from Øyvind Bakke and Ingeborg Hem"
date: "30.08 and 06.09 [PL], 31.08 and 07.09 [IL]"
output: #3rd letter intentation hierarchy
 beamer_presentation:
   keep_tex: yes
   fig_caption: false
   latex_engine: xelatex
---

```{r setup, include=FALSE}
showsol<-TRUE
library(knitr)
opts_chunk$set(tidy.opts=list(width.cutoff=68),results="hold",tidy=TRUE,warning=FALSE,error=FALSE,message=FALSE)
knitr::opts_chunk$set(echo = TRUE)
```

# Overview

## Learning material

* Textbook: Chapter 2.2, 3 and B.4.
* [Classnotes 30.08.2018](https://www.math.ntnu.no/emner/TMA4315/2018h/TMA4315M2H20180830.pdf)
* [Classnotes 06.09.2018](https://www.math.ntnu.no/emner/TMA4315/2018h/TMA4315M2H20180906.pdf)


---

## Topics
### [First week](#firstweek)

* Aim of multiple linear regression.
* Define and understand the multiple linear regression model - traditional and GLM way
* parameter estimation with maximum likelihood (and least squares)
* likelihood, score vector and Hessian (observed Fisher information matrix)
* big data implementation (if time)
* properties of parameter estimators
* assessing model fit (diagnostic), residuals, QQ-plots
* design matrix: how to code categorical covariates (dummy or effect coding), and how to handle interactions

---

### [Second week](#secondweek)

* What did we do last week?
* Statistical inference for parameter estimates 
    + confidence intervals, 
    + prediction intervals, 
    + hypothesis test, 
    + linear hypotheses
* analysis of variance decompositions and $R^2$, sequential ANOVA table
* DEVIANCE???
* model selection with AIC and variants


---

# Aim of multiple linear regression

1. **Explanation** Construct a model to help understand the relationship between a response and one or several explanatory variables.

2. **Prediction** Construct a model to predict the response from a set of (one or several) explanatory variables.

---

### Munich rent index

Munich, 1999: 3082 observations on 9 variables.

Interested in `rent`. Govt wants to know what a typical rent is in different areas of Munich

- what variables explain it
- how well can we predict it?

More information in Fahrmeir et. al., (2021) page 6.

---


### Munich rent index Variables

* `rent`: the net rent per month (in Euro).
* `rentsqm`: the net rent per month per square meter (in Euro).
* `area`: Living area in square meters.
* `yearc`: year of construction.
* `location`: quality of location: a factor indicating whether the location is average location, 1, good location, 2, and top location, 3.
* `bath`: quality of bathroom: a a factor indicating whether the bath facilities are standard, 0, or premium, 1.
* `kitchen`: Quality of kitchen: 0 standard 1 premium.
* `cheating`: central heating: a factor 0 without central heating, 1 with central heating.
* `district`: District in Munich.

---

```{r,echo=FALSE,eval=TRUE,warnings=FALSE,messages=FALSE,error=FALSE}
library("gamlss.data")
library(GGally)
ggpairs(rent99,lower = list(combo = wrap(ggally_facethist, binwidth = 0.5)))
```

---

**Interesting questions**

1. Is there a relationship between `rent` and `area`?
2. How strong is this relationship?
3. Is the relationship linear?
4. Are also other variables associated with `rent`? 
5. How well can we predict the rent of an apartment?
6. Is the effect of `area` the same on `rent` for apartments at average, good and top `location`? (interaction)

---


# Notation

${\bf Y}: (n \times 1)$ vector of responses (random variable) [e.g. one of the following: rent, rent pr sqm, weight of baby, ph of lake, volume of tree]

${\bf X}: (n \times p)$ design matrix [e.g. location of flat, gestation age of baby, chemical measurement of the lake, height of tree]

$\boldsymbol{\beta}: (p \times 1)$ vector of regression parameters (intercept included, so $p=k+1$)

$\boldsymbol{\varepsilon}: (n\times 1)$ vector of random errors. Used in "traditional way".

We assume that pairs $({\bf x}_i^T,y_i)$ ($i=1,...,n$) are measured from sampling units. That is, the observation pair $({\bf x}_1^T,y_1)$ is independent from $({\bf x}_2^T,y_2)$, and so on.


----

### Hands on: Munich rent index — response and covariates

From the list of variable and the statement of the questions, answer these questions:

- <span class="question">What can be response, and what covariates? (using what you know about rents)</span>
- <span class="question">What type of response(s) do we have? (continuous, categorical, nominal, ordinal, discrete, factors, ...). </span>
- <span class="question">What types of covariates? (continuous, categorical, nominal, ordinal, discrete, factors, ...) </span>
- <span class="question">Explain what the elements of `model.matrix` will be (Hint: coding of `location`) </span>


```{r , echo=FALSE, eval=FALSE, warning=FALSE, message=FALSE}
cat("Answers \n")
cat("a) rent or rentsqm for response and all other as covariates.\n")
cat("b) rent, rentsqm (continuous). \n")
cat("c) area (continuous), yearc (ordinal or continuous), location (categorical - ordinal), bath, kichen, cheating (categorical, nominal?), district (categorical - nominal).\n")
cat("d) our design matrix, contains numerical values for the continuous covariates, and dummy variable coding for categorical covariates (more about this later).\n")
```

```{r, echo=FALSE}
library("gamlss.data")
ds=rent99
# colnames(ds)
# summary(ds)
# dim(ds)
# head(ds)
# str(ds$location)
# contrasts(ds$location)

X=model.matrix(rentsqm~area+yearc+location+bath+kitchen+cheating+district,data=ds)
# head(X)
```

---

# Regression Models

---

# Model

## The traditional way

$${\bf Y}=X \boldsymbol{\beta}+\boldsymbol{\varepsilon}$$
is called a classical linear model if the following is true:

1. $\text{E}(\boldsymbol{\varepsilon})=\bf{0}$.

2. $\text{Cov}(\boldsymbol{\varepsilon})=\text{E}(\boldsymbol{\varepsilon}\boldsymbol{\varepsilon}^T)=\sigma^2\bf{I}$.

3. The design matrix has full rank, $\text{rank}({\bf X})=k+1=p$.

The classical _normal_ linear regression model is obtained if additionally

4. $\varepsilon\sim N_n({\bf{0}}, \sigma^2 \bf{I})$ holds. 

For random covariates these assumptions are to be understood conditionally on $\bf{X}$.

---

## The GLM way

Independent pairs $(Y_i, {\bf x}_i)$ for $i=1,\ldots,n$.

1. Random component: $Y_i \sim N(\mu_i, \sigma^2)$.
2. Systematic component: $\eta_i={\bf x}_i^T \boldsymbol{\beta}$.
3. Link function: linking the random and systematic component (linear predictor): Identity link and response function.
$\mu_i=\eta_i$.

---

### Questions

* <span class="question">Compare the traditional and GLM way. Have we made the same assumptions for both?
* <span class="question">What is the connection between each ${\bf x}_i$ and the design matrix? 
* <span class="question">What is "full rank"? Why is this needed? Example of rank less than $p$?
* <span class="question">Why do you think we move from traditional to GLM way? Could we not just let $\varepsilon$ be from binomial, Poisson, etc. distribution?

---

# Parameter estimation

In multiple linear regression there are two popular methods for estimating the regression parameters in $\beta$:

- maximum likelihood and 
- least squares. 

These two methods give the same estimator when we assume the normal linear regression model.

---

## Likelihood theory

### Likelihood $L(\beta)$

We assume that pairs of covariates and response are measured independently of each other: $({\bf x}_i,Y_i)$, and $Y_i$ follows the distribution specified above, and ${\bf x}_i$ is fixed.

$$L(\beta)=\prod_{i=1}^n L_i(\beta)=\prod_{i=1}^n f(y_i; \beta)$$
**Q**: <span class="question">fill in with the normal density for $f$ and the multiple linear regression model.

---

### Loglikelihood $l(\beta)$

We work with the log-likelihood because this makes the mathematics simpler

The main aim with the likelihood is to maximize it to find the maximum likelihood estimate, and since the log is a monotone function the maximum of the log-likelihood will be in the same place as the maximum of the likelihood.

$$
l(\beta)=\ln L(\beta)=\sum_{i=1}^n \ln L_i(\beta)=\sum_{i=1}^n l_i(\beta)
$$

Observe that the log-likelihood is a sum of individual contributions for each observation pair $i$.

**Q**: <span class="question">fill in with the normal density for $f$ and the multiple linear regression model.

---

### Rules for derivatives of vector

Härdle and Simes (2015), page 65.

* Let ${\beta}$ be a $p$-dimensional column vector of interest, 
* and let $\frac{\partial}{\partial {\beta}}$ denote the $p$-dimensional vector with partial derivatives wrt the $p$ elements of ${\beta}$. 
* Let ${\bf d}$ be a $p$-dimensional column vector of constants and 
* ${\bf D}$ be a $p\times p$ symmetric matrix of constants.

---

### Rules for derivatives of vector II

**Rule 1:**

$$\frac{\partial}{\partial \boldsymbol{\beta}}({\bf d}^T\boldsymbol{\beta})=\frac{\partial}{\partial \boldsymbol{\beta}}(\sum_{j=1}^p d_j \beta_j)={\bf d}$$
**Rule 2:**
$$\frac{\partial}{\partial \boldsymbol{\beta}}(\boldsymbol{\beta}^T{\bf D}\boldsymbol{\beta})=\frac{\partial}{\partial \boldsymbol{\beta}}(\sum_{j=1}^p \sum_{k=1}^p \beta_j d_{jk} \beta_k)=2{\bf D}\boldsymbol{\beta}$$

**Rule 3:**
The Hessian of the quadratic form $\boldsymbol{\beta}^T{\bf D}\boldsymbol{\beta}$ is

$$\frac{\partial^2 \boldsymbol{\beta}^T{\bf D}\boldsymbol{\beta}}{\partial \boldsymbol{\beta}\partial \boldsymbol{\beta}^T}= 2{\bf D}$$

---

### Score function $s(\beta)$

The score function is a $p\times 1$ vector, $s(\beta)$, with the partial derivatives of the log-likelihood with respect to the $p$ elements of the $\beta$ vector.

<!-- **Q**: Write down the rules for derivatives: chain rule, product rule, fraction rule, and in particular derivative of $\ln(x)$, $\exp(x)$ and $\frac{1}{x}$, you will need them now. (Yes, known from the R1 course in vgs.) -->
<!-- A: Chain rule: $\frac{d f(u(x))}{du}=\frac{df}{du}\cdot \frac{du}{dx}$, product rule: $(u\cdot v)'=u'\cdot v+u\cdot v'$, fraction rule: $(\frac{u}{v})'=\frac{u' \cdot v - u\cdot v'}{v^2}$, $\frac{d \ln(x)}{dx}=\frac{1}{x}$, $\frac{d\exp(x)}{dx}=\exp(x)$ and $\frac{d(\frac{1}{x})}{dx}=-\frac{1}{x^2}$. -->

$$
s(\beta)=\frac{\partial l(\beta)}{\partial \beta}=
\sum_{i=1}^n \frac{\partial l_i(\beta)}{\partial \beta}=
\sum_{i=1}^n s_i(\beta)
$$

Again, observe that the score function is a sum of individual contributions for each observation pair $i$.

**Q**: <span class="question">fill in for the multiple linear regression model.

---

To find the maximum likelihood estimate $\hat{\beta}$ we solve the set of $p$ equations:
$$s(\hat{\beta})=0$$

---

**Q**: <span class="question">fill in for the multiple linear regression model. Specify what the _normal equations_ are.

For the normal linear regression model, these equations $s(\hat{\beta})=0$ have a solution to be written on closed form. 

---

### Least squares and maximum likelihood (ML) estimator for ${\boldsymbol \beta}$:

$$
\hat\beta=({\bf X}^T{\bf X})^{-1} {\bf X}^T {\bf Y}
$$

**Q**: <span class="question">How can you see that least squares and ML gives the same estimator?

---

### Looking ahead: Hessian and Fisher information

But, for other distribution than the normal we get a set of non-linear equations when we look at $s(\hat{\beta})=0$, and then we will use the Newton-Raphson or Fisher Scoring iterative methods.


---

**Observed Fisher information matrix $H(\beta)$**

$$
H(\beta) = -\frac{\partial^2l(\beta)}{\partial\beta\partial\beta^T} = -\frac{\partial s(\beta)}{\partial\beta^T}
$$
- May be considered as a _local measure of information_ that the likelihood contains. 
- The higher the curvature of the log-likelihood near its maximum the more information is providd by the likelihood about the unknown parameter. 


**Q:** <span class="question">Calculate this for the multiple linear regression model. What is the dimension of $H(\beta)$?

---

In addition we also use the _expected Fisher information matrix $F(\beta)$_ which we may find in two ways, one is by taking the mean of the observed Fisher information matrix:

$$
F(\beta) = E\left( -\frac{\partial^2l(\beta)}{\partial\beta\partial\beta^T} \right).
$$

**Q:** <span class="question">Calculate this for the multiple linear regression model. What is the dimension of $F(\beta)$?

More on these later: for the moment, note the relationship between $\frac{\partial^2l(\beta)}{\partial\beta\partial\beta^T}$ and $Var(\hat{\beta})$

---

### Hands on: Munich rent index parameter estimates

Explain what the values under `Estimate` mean in practice.

```{r}
fit=lm(rentsqm~area+yearc+location+bath+kitchen+cheating,data=ds)
summary(fit)$coefficients
```


---


## Restricted maximum likelihood estimator for ${\boldsymbol \sigma}^2$

$$
\hat{\sigma}^2=\frac{1}{n-p} ({\bf Y}-{\bf X}{\hat{\beta}})^T ({\bf Y} - {\bf X}{\boldsymbol{\hat{\beta}}})=\frac{\text{SSE}}{n-p}
$$

The regression parameters $\beta$ are therefore our prime focus.

We will look at the parameter $\sigma^2$ as a nuisance parameter = parameter that is not of interest to us. 

---

To perform inference we need an estimator for $\sigma^2$.

The maximum likelihood estimator for $\sigma^2$ is $\frac{\text{SSE}}{n}$, which is found from maximizing the likelihood inserted our estimate of $\hat{\beta}$

$$
L(\hat{\beta},{\sigma^2})=(\frac{1}{2\pi})^{n/2}(\frac{1}{\sigma^2})^{n/2}\exp(-\frac{1}{2\sigma^2} ({\bf y}-{\bf X}\hat{\beta})^T({\bf y}-{\bf X}\hat{\beta}))
$$

$$
\begin{aligned}
l(\hat{\beta},{\sigma^2}) &= \text{ln}(L(\hat{\beta},{\sigma^2})) \\ 
&= -\frac{n}{2}\text{ln}(2\pi)-\frac{n}{2}\text{ln}\sigma^2-\frac{1}{2\sigma^2} ({\bf y}-{\bf X}\hat{\beta})^T({\bf y}-{\bf X}\hat{\beta})
\end{aligned}
$$


The score vector with respect to $\sigma^2$ is

$$
\frac{\partial l}{\partial \sigma^2}=0-\frac{n}{2\sigma^2}+\frac{1}{2\sigma^4}({\bf y}-{\bf X}\hat{\beta})^T({\bf y}-{\bf X}\hat{\beta})
$$

Solving $\frac{\partial l}{\partial \sigma^2}=0$ gives us the estimator

$$
\frac{1}{n}({\bf y}-{\bf X}\hat{\beta})^T({\bf y}-{\bf X}\hat{\beta})=\frac{\text{SSE}}{n}
$$ 

But, this estimator is biased.

<!---
To prove this you may use the trace-formula, that is $\text{E}({\bf Y}^T {\bf A}{\bf Y})=\text{tr}({\bf A}\text{Cov}({\bf Y}))+\text{E}({\bf Y})^T{\bf A}\text{E}({\bf Y})$, and we use that $\text{SSE}={\bf Y}^T ({\bf I}-{\bf H}){\bf Y}$. This was done in [class notes from TMA4267 - lecture 10](https://www.math.ntnu.no/emner/TMA4267/2017v/L10classnotes20170217.pdf)
--->

---

When an unbiased version is preferred, it is found using _restricted maximum likelihood_ (REML). We will look into REML-estimation in Module 7. In our case the (unbiased) REML estimate is

$$ \hat{\sigma}^2=\frac{1}{n-p}({\bf y}-{\bf X}\hat{\beta})^T({\bf y}-{\bf X}\hat{\beta})=\frac{\text{SSE}}{n-p}$$

The restricted maximum likelihood estimate is used in `lm`.

**Q:** <span class="question">What does it mean that the REML estimate is unbiased?</span>

<!---
  cat("Answer: \n")
cat("Residual standard error: 2.031 is sigmahat")
--->



<!-- ## Orthogonal parameters -->

<!-- **Q**: Why can we find $\hat{\boldsymbol{\beta}}$ first (from maximizing the likelihood wrt $\boldsymbol{\beta}$), and then replace $\boldsymbol{\beta}$ by $\hat{\boldsymbol{\beta}}$ in the likelihood and then maximize wrt $\sigma^2$? Why do we not have to maximize wrt both $\boldsymbol{\beta}$ and $\sigma^2$ simultaneously? That is, solve -->

<!-- $$ \begin{pmatrix} -->
<!-- \sum_{i=1}^n \frac{\partial l_i(\boldsymbol{\beta},\sigma^2)}{\partial \boldsymbol{\beta}} \\ \sum_{i=1}^n \frac{\partial l_i(\boldsymbol{\beta},\sigma^2)}{\partial \sigma^2} -->
<!-- \end{pmatrix}= \begin{pmatrix} {\bf 0} \\ 0 \end{pmatrix} -->
<!-- $$  -->
<!-- for $\boldsymbol{\beta}$ and $\sigma^2$ at the same time? -->

<!-- See class notes. -->

# Properties for the normal linear model

To be able to do inference we need to know about the properties of our parameter estimators in the (normal) linear model.

* Least squares and maximum likelihood estimator for ${\boldsymbol \beta}$:
\[ \hat{\beta}=({\bf X}^T{\bf X})^{-1} {\bf X}^T {\bf Y}\]
with $\hat{\beta}\sim N_{p}(\beta,\sigma^2({\bf X}^T{\bf X})^{-1})$.

* Restricted maximum likelihood estimator for ${\boldsymbol \sigma}^2$:
\[ \hat{\sigma}^2=\frac{1}{n-p}({\bf Y}-{\bf X}\hat{\beta})^T({\bf Y}-{\bf X}\hat{\beta})=\frac{\text{SSE}}{n-p}\]
with $\frac{(n-p)\hat{\sigma}^2}{\sigma^2} \sim \chi^2_{n-p}$.

* Statistic for inference about $\beta_j$, $c_{jj}$ is diagonal element $j$ of $({\bf X}^T{\bf X})^{-1}$.
\[ T_j=\frac{\hat{\beta}_j-\beta_j}{\sqrt{c_{jj}}\hat{\sigma}}\sim t_{n-p}\]
This requires that $\hat{\beta}_j$ and $\hat{\sigma}$ are independent (see below).

---

However, when we work with _large samples_ then $n-p$ becomes large and the $t$ distribution goes to a normal distribution, so we may use the standard normal in place of the $t_{n-p}$.

**Asymptotically** we have:
$$\hat{\beta}\sim N_{p}(\beta,\tilde{\sigma}^2({\bf X}^T{\bf X})^{-1})$$.
and
$$ T_j=\frac{\hat{\beta}_j-\beta_j}{\sqrt{c_{jj}}\tilde{\sigma}}\sim N(0,1)$$
where $\tilde{\sigma}^2=\frac{\text{SSE}}{n}$ (the ML estimator).

**Q:** <span class="question">Pointing forwards: do you see any connection between the covariance matrix of $\hat{\boldsymbol{\beta}}$ and the Fisher information?

----

## Are $\hat{\boldsymbol{\beta}}$ and SSE are independent? (optional)

Can be proven using knowledge from TMA4267 on properties of the multivariate normal distribution.

Independence:  Let ${\bf X}_{(p \times 1)}$ be a random vector from
    $N_p({\boldsymbol \mu},{\boldsymbol \Sigma})$. Then ${\bf A}{\bf X}$ and ${\bf B}{\bf X}$ are independent 
iff ${\bf A}{\boldsymbol \Sigma}{\bf B}^T={\bf 0}$.

We have: 

* ${\bf Y}\sim N_n({\bf X}{\boldsymbol \beta},\sigma^2{\bf I})$

* ${\bf AY}=\hat{{\boldsymbol \beta}}=({\bf X}^T{\bf X})^{-1}{\bf X}^T {\bf Y}$, and

* ${\bf BY}=({\bf I}-{\bf H}){\bf Y}$.

* Now ${\bf A}\sigma^2{\bf I}{\bf B}^T=\sigma^2 {\bf A}{\bf B}^T=\sigma^2 ({\bf X}^T{\bf X})^{-1}{\bf X}^T ({\bf I}-{\bf H})={\bf 0}$

* since ${\bf X}({\bf I}-{\bf H})={\bf X}-{\bf HX}={\bf X}-{\bf X}={\bf 0}$.

* We conclude that $\hat{{\boldsymbol \beta}}$ is independent of $({\bf I}-{\bf H}){\bf Y}$,

* and, since SSE=function of $({\bf I}-{\bf H}){\bf Y}$: SSE=${\bf Y}^T({\bf I}-{\bf H}){\bf Y}$,

* then $\hat{{\boldsymbol \beta}}$ and SSE are independent, and the result with $T_j$ being t-distributed with $n-p$ degrees of freedom is correct.

Remark: a similar result will exist for GLMs, using the concept of _orthogonal parameters_. 

---

# Checking model assumptions 

---

# Checking model assumptions 

In the normal linear model we have made the following assumptions.

1. Linearity of covariates: ${\bf Y}={\bf X}\beta+\varepsilon$. Problem: non-linear relationship?

2. Homoscedastic error variance: $\text{Cov}(\varepsilon)=\sigma^2 {\bf I}$.

3. Uncorrelated errors: $\text{Cov}(\varepsilon_i,\varepsilon_j)=0$.

4. Additivity: ${\bf Y}={\bf X}\beta+\varepsilon$

5. Assumption of normality: $\varepsilon \sim N_n({\bf 0},\sigma^2{\bf I})$

The same assumptions are made when we do things the GLM way for the normal linear model. 

In addtion the following might cause problems:

* Outliers
* High leverage points
* Collinearity

---

## Residuals

If we assume the normal linear model then we know that the  residuals ($n\times 1$ vector)
$$\hat{\varepsilon}={\bf Y}-{\bf \hat{Y}}=({\bf I}-{\bf H}){\bf Y}$$
has a normal (singular) distribution with mean $\text{E}(\hat{\varepsilon})={\bf 0}$ and covariance matrix $\text{Cov}(\hat{ \varepsilon})=\sigma^2({\bf I}-{\bf H})$ where ${\bf H}={\bf X}({\bf X}^T{\bf X})^{-1}{\bf X}^T$.

This means that the residuals (possibly) have different variance, and may also be correlated.

---

Our best guess for the error $\varepsilon$ is the residual vector $\hat{\varepsilon}$, and we may think of the residuals as _predictions of the errors_. 
Be aware: don't mix errors (the unobserved) with the residuals ("observed"). 

But, we see that the residuals are not independent and may have different variance, therefore we will soon define variants of the residuals that we may use to assess model assumptions after a data set is fitted.

**Q:** <span class="question">How can we say that the residuals can have different variance and may be correlated? Why is that a problem?

---

We would like to check the model assumptions 

- we see that they are all connected to the error terms. 

But, but we have not observed the error terms $\varepsilon$

However, we have made "predictions" of the errors - our residuals. And, we want to use our residuals to check the model assumptions. 

---

We want to check that our errors are 

- independent, 
- homoscedastic (same variance for each observation), 
- not dependent on our covariates

We want to use the residuals (observed) in place of the errors (unobserved). 

It would have been great if the residuals have these properties when the underlying errors have. 

Enter _standardized_ or _studentized residuals_.

---

### Standardized residuals:

$$r_i=\frac{\hat{\varepsilon}_i}{\hat{\sigma}\sqrt{1-h_{ii}}}$$
where $h_ii$ is the $i$th diagonal element of the hat matrix ${\bf H}$.

In R you can get the standardized residuals from an `lm`-object (named `fit`) by `rstandard(fit)`.

---

### Studentized residuals:

$$r^*_i=\frac{\hat{\varepsilon}_i}{\hat{\sigma}_{(i)}\sqrt{1-h_{ii}}}$$
where $\hat{\sigma}_{(i)}$ is the estimated error variance in a model with observation number $i$ omitted. This seems like a lot of work, but it can be shown that it is possible to calculated the studentized residuals directly from the standardized residuals:
$$r^*_i=r_i (\frac{n-p-1}{n-p-r_i^2})^{1/2}$$

In R you can get the studentized residuals from an `lm`-object (named `fit`) by `rstudent(fit)`.

---

### Plotting residuals - and what to do when assumptions are violated?

Some important plots

1. Plot the residuals, $r^*_i$ against the predicted values, $\hat{y}_i$.

- Dependence of the residuals on the predicted value: wrong regression model?
                
- Nonconstant variance: transformation or weighted least squares is needed?

2. Plot the residuals, $r^*_i$, against predictor variable or functions of predictor variables. 

- Trend suggest that transformation of the predictors or more terms are needed in the regression.

---

3. Assessing normality of errors: QQ-plots and histograms of residuals. 

Tests for normality can be used, but they can be useless: for small sample sizes the test is not powerful and for large sample sizes even very small deviances from normality will be labelled as significant. 

4. Plot the residuals, $r^*_i$, versus time or collection order (if possible). Look for dependence or autocorrelation.

Residuals can be used to check model assumptions, and also to _discover outliers_.

---

### Diagnostic plots in R

```{r echo = FALSE}
library(ggplot2)
```


For simplicity we use the Munich rent index with `rent` as response and only `area` as the only covariate.

\footnotesize

```{r tidy = TRUE, echo=showsol}
fit <- lm(rent ~ area, data=rent99) # Run a regression analysis
format(head(fortify(fit)), digits = 4L)
```

(`ggplot2::fortify.lm` creates a dataframe from an `lm`-object, which can be used to plot diagnostic plots. `ggplot` will do this automatically when asked to plot)

---

#### Residuals vs fitted values

A plot with the fitted values of the model on the x-axis and the residuals on the y-axis shows if the residuals have non-linear patterns. The plot can be used to test the assumption of a linear relationship between the response and the covariates. If the residuals are spread around a horizontal line with no distinct patterns, it is a good indication on no non-linear relationships, and a good model.

Does this look like a good plot for this data set?

---

```{r tidy = TRUE,echo=showsol, fig.height=4}
ggplot(fit, aes(.fitted, .stdresid)) +
  geom_point(pch = 21) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  geom_smooth(se = FALSE, col = "red", size = 0.5, method = "loess") +
  labs(x = "Fitted values", y = "Standardized residuals", title = "Fitted values vs standardized residuals", subtitle = deparse(fit$call))

```

<!-- (Ok linear assumption, but not constant spread.) -->

----

#### Normal Q-Q

This plot shows if the residuals are Gaussian (normally) distributed. If they follow a straigt line it is an indication that they are, and else they are probably not.

```{r tidy = TRUE,echo=showsol, fig.height=4}
ggplot(fit, aes(sample = .stdresid)) +
  stat_qq(pch = 19) +
  geom_abline(intercept = 0, slope = 1, linetype = "dotted") +
  labs(x = "Theoretical quantiles", y = "Standardized residuals", title = "Normal Q-Q", subtitle = deparse(fit$call))
```

---


#### Scale-location

This is also called spread-location plot. It shows if the residuals are spread equally along the ranges of predictors. Can be used to check the assumption of equal variance (homoscedasticity). A good plot is one with a horizontal line with randomly spread points.

Is this plot good for your data?

---

```{r, fig.height=4}
ggplot(fit, aes(.fitted, sqrt(abs(.stdresid)))) +
  geom_point() +
  geom_smooth(se = FALSE, col = "red", size = 0.5, method = "loess") +
  labs(x = "Fitted values", y = expression(sqrt("Standardized residuals")), title = "Scale-location", subtitle = deparse(fit$call))

```

<!-- (Confirms our observation of not constant variance.) -->

----

#### Residual vs Leverage

This plot can reveal influential outliers. 

Not all outliers are influential in linear regression; the results might not change if they are removed from the data set

Influential outliers can be seen as observations that does not get along with the trend in the majority of the observations. 

---

Cook's distance is the Euclidean distance between the $\mathbf{\hat{y}}$ (the fitted values) and $\mathbf{\hat{y}}_{(i)}$ (the fitted values calculated when the $i$-th observation is omitted from the regression). 

This is then a measure on how much the model is influences by observation $i$. 

The distance is scaled, and a rule of thumb is to examine observations with Cook's distance larger than 1, and give some attention to those with Cook's distance above 0.5.

Leverage is defined as the diagonal elements of the hat matrix, i.e., the leverage of the $i$-th data point is $h_{ii}$ on the diagonal of $\mathbf{H = X(X^TX)^{-1}X^T}$. A large leverage indicated that the observation ($i$) has a large influence on the estimation results, and that the covariate values ($\mathbf{x}_i$) are unusual.

---

```{r, tidy = TRUE,echo=showsol, fig.height=4}
ggplot(fit, aes(.hat, .stdresid)) +
  geom_smooth(se = FALSE, col = "red", size = 0.5, method = "loess") +
  geom_point(aes(size = .cooksd)) +
  scale_size_continuous("Cook's dist.") +
  labs(x = "Leverage", y = "Standardized residuals", title = "Residuals vs Leverage", subtitle = deparse(fit$call))
```

(Some observations does not fit our model, but if we fit a more complex model this may change.)

---

# Categorical covariates - dummy and effect coding

(read for yourself - topic of ILw1)



# Interactions 

(if we have time)

To illustrate how interactions between covariates can be included we use the `ozone` data set from the `ElemStatLearn` library. This data set is measurements from 1973 in New York and contains 111 observations of the following variables:

* `ozone` : ozone concentration (ppm)
* `radiation` : solar radiation (langleys)
* `temperature` : daily maximum temperature (F)
* `wind` : wind speed (mph)

---

# Ozone

We start by fitting a multiple linear regression model to the data, with `ozone` as our response variable and `temperature` and `wind` as covariates.

```{r,echo=FALSE,eval=TRUE,messages=FALSE,warnings=FALSE,error=FALSE,results="hold"}
#  The ElemStatLearn package isn't on CRAN, so download the data directly
# library(ElemStatLearn)
ozone <- read.table("https://hastie.su.domains/ElemStatLearn/datasets/ozone.data",
                    header=TRUE)

library(knitr)
library(kableExtra)
kable(head(ozone))
ozone.lm = lm(ozone~temperature+wind, data=ozone)
```

---

\small
```{r,echo=FALSE,eval=TRUE,messages=FALSE,warnings=FALSE,error=FALSE,results="hold"}
summary(ozone.lm)
```
\normalsize

---

The model can be written as:
$$Y = \beta_0 + \beta_1 x_t + \beta_2 x_w + \varepsilon$$
In this model we have assumed that increasing the value of one covariate is independent of the other covariates. For example: by increasing the `temperature` by one-unit always increases the response value by $\beta_2 \approx 1.651$, regardless of the value of `wind`.

---

However, one might think that the covariate `wind` (wind speed) might act differently upon `ozone` for different values of `temperature` and vice verse.
$$
\begin{aligned} Y &= \beta_0 +  \beta_1 x_t + \beta_2 x_w + \beta_3\cdot(x_t  \cdot x_w) +\varepsilon \\ 
&= \beta_0 +  (\beta_1 + \beta_3 x_w) \cdot x_t + \beta_2 x_w + \varepsilon \\ 
&= \beta_0 + \beta_1 x_t + (\beta_2 + \beta_3 x_t) \cdot x_w + \varepsilon \end{aligned}.
$$
We fit this model in `R`. An interaction term can be included in the model using the `*` symbol.

**Q:** <span class="question">Look at the `summary` below. Is this a better model than without the interaction term? It the term significant?

---

\footnotesize
```{r}
ozone.int = lm(ozone~temperature+wind+ temperature*wind, data=ozone)
summary(ozone.int)
```
\normalsize

---

Below we see that the interaction term is highly significant. The $p$-value is very small, so that there is strong evidence that $\beta_3 \neq 0$. Furthermore, $R^2_{\text{adj}}$ has increased, indicating that more of the variability in the data has been explained by the model (than without the interaction).


---

*Interpretation of the interaction term:*

* If we now increase the `temperature` by $10^{\circ}$ F, the increase in `wind` speed will be $$(\hat \beta_1+\hat \beta_3 \cdot x_w) \cdot 10 = (4.0 -0.22 \cdot x_w) \cdot 10 = 40-2.2 x_w \text{ units}.$$

* If we increase the `wind` speed by 10 mph, the increase in `temperature` will be $$(\hat \beta_2 + \hat \beta_3 \cdot x_t) \cdot 10 = (14 -0.22 \cdot x_t) \cdot 10 = 140-2.2 x_t \text{ units}.$$

---

## The hierarchical principle

It is possible that the interaction term is highly significant, but the main effects are not.

In our `ozone.int` model above: the main effects are `temperature` and `wind`. The hierarchical principle states that if we include an interaction term in our model, the main effects are also to be included, even if they are not significant. This means that if the coefficients $\hat \beta_1$ or $\hat \beta_2$ would be insignificant, while the coefficient $\hat \beta_3$ is significant, $\hat \beta_1$ and $\hat \beta_2$ should still be included in the model.


---

## The hierarchical principle: why?

A model with interaction terms, but without the main effects is hard to interpret.

Removing a main effect is the same as setting $\beta_1=0$

$$
Y = \beta_0 +  \beta_1 x_t + \beta_2 x_w + \beta_3\cdot(x_t  \cdot x_w) +\varepsilon
$$

i.e. saying the slope of the $x_t$ effect is 0 when $x_w = 0$.

---

## Interactions between qualitative (discrete) and quantitative (continuous) covariates

We create a new variable `temp.cat` which is a `temperature` as a qualitative covariate with two levels and fit the model:
$$\begin{aligned}y&=\beta_0 + \beta_1 x_w + \begin{cases} \beta_2 + \beta_3  x_w  &\text{ if temperature="low"}\\ 0 &\text{ if temperature = "high"}\end{cases} \\\\ &= \begin{cases} (\beta_0 + \beta_2) + (\beta_1 + \beta_3) \cdot x_w &\text{ if temperature="low"}\\ \beta_0 + \beta_1 x_w &\text{ if temperature="high""} \end{cases} \end{aligned}$$

---

# Ozone: make temperature categorical

```{r,echo=showsol}
temp.cat = ifelse(ozone$temperature < mean(ozone$temperature), "low", "high")
ozone2 = cbind(ozone, temp.cat)
kable(head(ozone2))
```

---

# Model with interaction

```{r,echo=showsol}
ozone.int2 = lm(ozone~wind + temp.cat+ temp.cat*wind, data=ozone2)
summary(ozone.int2)$coefficients
```


---

```{r,echo=showsol, fig.height=4}
interceptlow = coef(ozone.int2)[1]+coef(ozone.int2)[3]
slopelow = coef(ozone.int2)[2]+coef(ozone.int2)[4]
intercepthigh = coef(ozone.int2)[1]
slopehigh = coef(ozone.int2)[2]
ggplot(ozone)+geom_line(aes(y=interceptlow+slopelow*wind, x = wind), col="blue")+geom_line(aes(y=intercepthigh+slopehigh*wind, x=wind), col="red")+ylab("ozone")
```



