---
subtitle: "TMA4268 Statistical Learning V2018"
title: "Module 3: LINEAR REGRESSION"
author: "Mette Langaas and Julia Debik, Department of Mathematical Sciences, NTNU"
date: "week 4 2018 (Version 23.05.2018)"
output: 
  # html_document:
  #   toc: true
  #   toc_float: true
#  beamer_presentation:
###    incremental: true # or >* for one at a time
#  slidy_presentation:
#    font_adjustment: +1  
  pdf_document:
   # toc: true
   # toc_depth: 2
editor_options: 
  chunk_output_type: console
fig.width: 6.5
fig.height: 6.5
---

Last changes: (23.05 typo p-q from Rasmus, 22.05 typos from Elisabeth and html toc restructuring.)

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

# Introduction

###Learning material for this module:

* James et al (2013): An Introduction to Statistical Learning. Chapter 3.  

Remark: there are no classnotes for this module since I wrote on the slides in class. Screencast from the lectures are available under Course Modules (column with "Dates") on Blackboard.

---

<!-- ## Will be added after the lecture(s): -->

<!-- * [Module3PartA: Classnotes 22.01.2018](https://www.math.ntnu.no/emner/TMA4268/2018v/notes/TMA4268ClassnotesM3L1.pdf).  -->
<!-- * [Module3PartB: Classnotes 24.01.2018](https://www.math.ntnu.no/emner/TMA4268/2018v/notes/TMA4268ClassnotesM3L2.pdf).  -->

### Move to:

* [Part A: Simple linear regression and introduction to multiple linear regression ](#PartA)
* [Part B: Multiple linear regression - continued](#PartB)
* [Recommended exercises](#RecEx)
* [Further reading](#further)
* [Packages to install before knitting this R Markdown file](#Rpackages)

---

## Topics in this module

**Part A: Simple linear regression and introduction to multiple linear regression**

* Aim of linear regression
* Simple linear regression 
    + model
    + parameter estimation
    + confidence intervals
    + single hypothesis testing: set-up and $p$-values
    + model fit
* Multiple linear regression
    + model
    + parameter estimation with least squares 
    + properties of parameter estimators

---

**Part B: Multiple linear regression - continued**

* So far: simple linear regression - and multiple linear regression - model, estimators
* Statistical inference 
    + estimators and properties
    + confidence intervals and hypothesis tests
    + significance of regression: F-test
    + prediction and prediction intervals 
* Model assessement and selection
    + $R^2$
    + subset selection 
    + diagnostic plots - studentized residuals and leverages
    
---

* Special cases
    + qualitative predictors: dummy coding 
    + non-additivity: including interactions
    + nonlinearity: polynomials
    + multiple linear regression vs. KNN-regression 
* Optional material
    + projection matrices and geometry of least squares
 
**Recommended exercises**

**Further reading**

**Rpackages to install**

---

# <a id="PartA"> Part A: Simple linear regression and introduction to multiple linear regression</a>

## Aim of linear regression

1. Construct a model to help understand the relationship between _one response_ and _one or several explanatory variables_. [Correlation, or cause and effect?]

2. Construct a model to predict the _reponse_ from a set of (one or several) _explanatory variables_. [More or less "black box"]

Is linear regression dull? Maybe, but very useful and widely used. Important to _understand_ because many learning methods can be seen as generalization of linear regression.

Linear regression is a supervised and parametric method.

---

## Motivating example: Munich rent index

Munich, 1999: 3082 observations on 9 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.

More information in Fahrmeir et. al., (2013) page 5 (textbook used in TMA4267 Linear statistical models).

---

```{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 appartment?
6. Is the effect of `area` the same on `rent` for appartments at average, good and top `location`? (interaction)

---

# Simple Linear Regression

* One quantitative response $Y$ is modelled with
* from _one covariate_ $x$ (=simple), 
* and the relationship between $Y$ and $x$ is assumed to be _linear_. 

---

## Model and assumptions
$$Y = f(x)+ \varepsilon= \beta_0 + \beta_1 x + \varepsilon$$

 * $Y$ is a _quantitative_ (continuous) response variable.
 * $\beta_0$ is the _intercept_. It is the value of $Y$ when $x=0$.
 * $\beta_1$ is the _regression coefficient_ or _regression parameter_ (slope). The slope represents the average increase in $Y$ given a one-unit increase in $x$.
 * $x$ is the covariate. It can be continuous or discrete.
 * $\varepsilon$ is the _error term_ also called the measurement noise. It replaces all the unobserved covariates that influence $Y$.
      + $\text{E}(\varepsilon)=0$ 
      + $\text{Var}(\varepsilon)=\sigma_{\varepsilon}^2$ (for all $x$).
      + Often we also assume that $\varepsilon$ is normally distributed. 
 
Here the _model parameters_ $\beta_0$, $\beta_1$ and $\sigma_{\varepsilon}^2$ are unknown and have to be estimated from _observed data_.
 
---

## Estimation and prediction

* _Fitting_ a model is to find the best estimates for the model parameters.
* This is based on data: independent pairs $(x_i,Y_i)$, $i=1,\ldots,n$.
* The _least squares method_ is used.
* We can predict the value of the response for a (new) observation of the covariate at $x_0$.
$$\hat{y} = \hat{\beta}_0 + \hat{\beta_1}x_0.$$

We use the $\hat{}$ symbol for estimated or predicted values.

---

### Example continued: Munich rent index

Assume a linear relationship between the monthly net rent (`rent`) and the living area (`area`) in square meters. 

In R: fit simple linear regression model using the `lm` function to the data,  and plot the data together with the _fitted model_.

We see that the model fits the data quite well. It captures the essence. It looks that a linear relationship between the response `rent` and covariate `area` is a good approximation. 

---

```{r,eval=TRUE,echo=FALSE}
library(ggplot2)
munich1.lm = lm(rent~area, data=rent99)
munich1.beta0 = coef(munich1.lm)[1]
munich1.beta1 = coef(munich1.lm)[2]
ggplot(rent99, aes(x=area, y=rent))+geom_point(size=0.7)+
  geom_line(aes(x=area, y=munich1.beta0+area*munich1.beta1), 
            col="blue", size=0.8)
```

---

**Q:** 

* The blue line gives the estimated model. Explain what the line means in practice. 
* How does this relate to the _true_ (population) model? 
* The spread around the fit does not look constant. What does this mean? Is that a problem?

```{r, echo=FALSE, results='asis', eval=showsolA}
cat("
**A**: 

* If we compare two appartment where one has area 1 sq metre more then the other, we would expect the rent to be $\\hat{\\beta}$ Euro higher.
* The true model: tja, might also be approximately linear on the population level (all appartments in Munich), and it might be different from the blue line. To give the population level `area` vs `rent` we need to plot the relationship between all units in the population - and if we assume the relationship is linear that line might deviate from our blue line.
* Yes, in linear regression we assume the variance of $\\varepsilon$ to be constant - the same for all `area`s.
")
```

---

## Parameter estimates: least squares

Assume that we have a data set consisting of $n$ observation pairs $(x_i, Y_i)$ for $i = 1,2,..., n$. 
We usually assume the _pairs are independent_.

Further, assume we fit a linear model to our data. For a single pair of observations $(x_i, y_i)$ it can be written as $\hat{Y}_i = \hat{\beta}_0 + \hat{\beta}_1 x_i$. 

The $i$-th *residual* of the model is the difference between the $i$-th *observed* response value and the $i$-th *predicted* value, and is written as: 
$$e_i = Y_i - \hat{Y}_i.$$

We may regard the residuals as _predictions_ (not estimates) of the error terms $\varepsilon_i$.

Remark: the error terms are random variables and can not be estimated - they can be predicted.

---

The _residual sum of squares (RSS)_ is the squared sum of all residuals 
$$\begin{aligned} \text{RSS} &= e_1^2+e_2^2+...+e_n^2 \\ &= (Y_1 - \hat{\beta}_0 - \hat{\beta}_1 x_1)^2  + (Y_2 -\hat{\beta}_0 - \hat{\beta}_1 x_2)^2+ ... + (Y_n -\hat{\beta}_0 - \hat{\beta}_1 x_n)^2  \end{aligned}$$

We find the parameter estimates for $\beta_0$ and $\beta_1$ by minimizing the RSS.


**Least squares estimators:**

$$\hat{\beta}_0 = \bar{Y}-\hat{\beta}_1 \bar{x}$$
and
$$\hat{\beta}_1 = \frac{\sum_{i=1}^n(x_i-\bar{x})(Y_i-\bar{Y})}{\sum_{i=1}^n(x_i-\bar{x})^2},$$
where 
$\bar{Y} = \frac{1}{n} \sum_{i=1}^n Y_i$ and $\bar{x} = \frac{1}{n}\sum_{i=1}^n x_i$ are the sample means.

---

### Example continued
The coefficient estimates from our simple linear model fit to the Munich rent index data can be printed using the `coef` function in R:
```{r}
coef(munich1.lm)
```

**Q**: Explain what these two values mean.

---

## Accuracy of the parameter estimates

(We will _derive a general version_ of these formulas for multiple linear regression, because without matrix notation this is very cumbersome - e.g. ["Egenskaper til estimatorene" for the TMA4240/TMA4245 Thematic pages](https://wiki.math.ntnu.no/tma4245/tema/begreper/regression).)

The standard errors of the estimates are given by the following formulas:
$$\text{Var}(\hat{\beta}_0)=\text{SE}(\hat{\beta}_0)^2 = \sigma_{\varepsilon}^2 \Big [ \frac{1}{n} + \frac{\bar{x}^2}{\sum_{i=1}^n (x_i -\bar{x})^2} \Big]$$
and
$$\text{Var}(\hat{\beta}_1)=\text{SE}(\hat{\beta}_1)^2 = \frac{\sigma_{\varepsilon}^2}{\sum_{i=1}^n (x_i-\bar{x})^2}.$$
In addition $\text{Cov}(\hat{\beta_0},\hat{\beta_1})$ is in general different from zero.

We see how the estimated coefficient will vary when our experiment is repeated.

Remark: SE=standard error vs SD=standard deviation - since we write SE of a random variable this means the standard deviation of this variable (square root of variance) - so for us SD and SE will be the same. 

In general much confusion about SD and SE in articles: <https://www.ncbi.nlm.nih.gov/pmc/articles/PMC1255808/>

---

**Q: The best design**

Assume that we know the value of $\sigma_{\varepsilon}^2$. Now, observe that

$$\text{SE}(\hat{\beta}_1)^2 = \frac{\sigma_{\varepsilon}^2}{\sum_{i=1}^n (x_i-\bar{x})^2}$$ is only dependent on the _design_ of the $x_i$'s. 

* Would we like the $\text{SE}(\hat{\beta}_1)^2$ large or small? Why?
* If it is possible for us to choose the $x_i$'s, which strategy should we use to choose them? 
* Assume $x$ can take values from 1 to 10 and we choose $n=10$ values. Which is the best design?
    + evenly in a grid: $[1,2,3,4,5,6,7,8,9,10]$.
    + only lower and upper value: $[1,1,1,1,1,10,10,10,10,10]$.
    + randomly drawn from a uniform distribution on $[1,10]$.

---

```{r,echo=TRUE,eval=TRUE}
x1=seq(1:10)
x2=c(rep(1,5),rep(10,5))
x3=runif(10,1,10)

sd1=sqrt(1/sum((x1-mean(x1))^2))
sd2=sqrt(1/sum((x2-mean(x2))^2))
sd3=sqrt(1/sum((x3-mean(x3))^2))

print(c(sd1,sd2,sd3))
```

**A**: the second design - all observations at extremes.

---

## Residual standard error

Remember - residual sum of squares: $\text{RSS}=\sum_{i=1}^n (Y_i-\hat{\beta}_0-\hat{\beta_1}x_{i})^2$.

The residual standard error, RSE, is given by
$$\text{RSE}  =\sqrt{\frac{1}{n-2} \text{RSS}} = \sqrt{\frac{1}{n-2}\sum_{i=1}^n (Y_i -\hat{Y}_i)^2}$$
and is an estimate of $\sigma_{\varepsilon}$, that is, the standard deviation of the error term $\varepsilon$. This is the socalled _restricted maximum likelihood estimator_, and is unbiased.

This is related to the amount the response variables deviate from the estimated regression line. Recall that we will always have observations with noise.

RSE shows the lack of fit of the model to the data. It is measured in units of $Y$, hence is value may be hard to interpret. 

---

If we assume that the simple linear regression is a good model, observation pairs $(x_i,Y_i)$ are independent for $i=1,\ldots,n$ and that $\varepsilon_i\sim N(0,\sigma_{\varepsilon}^2)$ then it can be shown that

$$\frac{\text{RSE}^2(n-2)}{\sigma_{\varepsilon}^2}= \frac{\sum_{i=1}^n (Y_i -\hat{Y}_i)^2}{\sigma_{\varepsilon}^2}\sim \chi^2_{n-2}$$

### Example continued
In R we can get a summary of our fitted linear model, by calling the `summary` function. In this outprint we see the estimated coefficient values in the first column, and the estimated standard errors in the second column. Thus $\hat{\text{SE}}(\hat{\beta}_0) = 8.6135$ and $\hat{\text{SE}}(\hat{\beta}_1) = 4.8215$.

---

\small
```{r}
summary(munich1.lm)
```
\normalsize

---

To illustrate this point further, we fit four models to our Munich rent index data set. Each of the models has been fit using a random sample of a fraction of our observations (1/4). We see how our fitted line changes given a "new" set of observations.

```{r, echo=FALSE}
set.seed(200)

munich.n1 = dim(rent99)[1]
munich.h = sample(munich.n1, munich.n1/2)
munich.1hds = rent99[munich.h, ]
munich.2hds = rent99[-munich.h, ]
munich.n2 = dim(munich.1hds)[1]
munich.q = sample(munich.n2, munich.n2/2)
munich.s1 = munich.1hds[munich.q,]
munich.s2 = munich.1hds[-munich.q,]
munich.s3 = munich.2hds[munich.q,]
munich.s4 = munich.2hds[-munich.q,]

munich.lm.s1 = lm(rent~area, data=munich.s1)
munich.lm.s2 = lm(rent~area, data=munich.s2)
munich.lm.s3 = lm(rent~area, data=munich.s3)
munich.lm.s4 = lm(rent~area, data=munich.s4)

ggplot(rent99, aes(x=area, y=rent))+geom_point(alpha=0.8, size=0.5)+
  geom_line(aes(x=area, y=coef(munich.lm.s1)[1]+area*coef(munich.lm.s1)[2]), 
            col="navy")+ geom_line(aes(x=area, y=coef(munich.lm.s2)[1] + coef(munich.lm.s2)[2]*area), 
            col="blue")+ geom_line(aes(x=area, y=coef(munich.lm.s3)[1]+coef(munich.lm.s3)[2]*area), 
            col="purple")+geom_line(aes(x=area, y=coef(munich.lm.s4)[1]+coef(munich.lm.s4)[2]*area), col="red")
```

---

## Distribution of parameter estimators
(We will _derive a general version_ for multiple linear regression.)

* Let $\hat{\beta_j}$ be the least squares estimator for $\beta_j$, and that 
* $\text{Var}(\hat{\beta_j})$ and $\widehat{\text{Var}}(\hat{\beta_j})$ is as given above. 

Then $\hat{\beta}_j \sim N(\beta_j,\text{Var}(\hat{\beta}_j))$ and 
$$T_j =\frac{\hat{\beta}_j-\beta_j}{\sqrt{\widehat{\text{Var}}(\hat{\beta}_j)}} \sim t_{n-2}$$
The $t$-distribution with $n-2$ degrees of freedom.

See [Rintermediate-sol.html](https://wwww.math.ntnu.no/emner/TMA4268/2018v/1Intro/Rintermediate-sol.html) for more on the $t$-distribution in R.

For the slope in the simple linear regression this is shown, together with inference for $\beta_1$, in this [video from TMA4240/TMA4245 Statistics (in Norwegian - but formulas are still understood in English).](https://mediasite.ntnu.no/Mediasite/Play/2e9a209c58874e75bd47e3c5e0b7b4e81d?catalog=0fce6173-7a98-4db7-84b7-50cba3a3a341)

---

```{r,tidy=FALSE, echo=FALSE}
require(ggplot2)
df=10
ggplot(data.frame(x=c(-4,4)), aes(x))+ 
  theme(axis.title.x=element_blank(),axis.title.y=element_blank(),axis.text.x=element_blank(),axis.text.y=element_blank(),axis.ticks = element_blank())+
  stat_function(fun=function(x) dt(x,df), geom="line", colour="blue")+
  geom_vline(xintercept=qt(0.975,df),color="red")+
  geom_vline(xintercept=-qt(0.975,df),color="red")+
  geom_hline(yintercept=0)+
  annotate("text",x=qt(0.975,df)+0.6,y=0.03,label="alpha/2",parse=TRUE,size=5)+
  annotate("text",x=-qt(0.975,df)-0.6,y=0.03,label="alpha/2",parse=TRUE,size=5)+
  annotate("text",x=0,y=0.03,label="1-alpha",parse=TRUE,size=5)
  
``` 

The figure shows the $t_{10}$ distribution with $\alpha=0.05$. 

---

## Confidence intervals

The $t$-distribution can be used to create confidence intervals for the regression parameters. The lower and upper limits of a 95\% confidence intervals for $\beta_j$ are
$$\hat{\beta}_j \pm t_{\alpha/2,n-2} \cdot\text{SE} (\hat{\beta}_j) \quad j=0, 1.$$
The interpretation of this confidence is that: (before we have contructed the interval) there is a 95\% chance that the interval 
will contain the *true* value of $\beta_j$.

If $n$ is large, the normal approximation to the $t$-distribution can be used (and is used in the textbook).

**Q:** Calculate the confidence intervals for the slope parameter in the `munich1.lm` model by finding the numbers you need from the `summary` output. Here $t_{0.025,n-2}$=`r round(qt(0.975,3080),3)`.

---

We can find confidence intervals in `R` using the `confint` function on `lm` object:
```{r}
confint(munich1.lm)
```

---

Based on the joint distribution of the intercept and slope it is possible to find the distribution for the linear predictor $\hat{\beta}_0+\hat{\beta}_1 x$, and then confidence intervals for $\beta_0+\beta_1 x$. 

```{r, echo=FALSE}
library(ggpubr)
munich4.lm = lm(rent~area, data=rent99[1:100,])
const=qt(0.975,98)
munich4.beta0.min = coef(munich4.lm)[1]-const*64.285
munich4.beta0.max = coef(munich4.lm)[1]+const*64.285

munich4.beta1.min = coef(munich4.lm)[2]-const*1.310
munich4.beta1.max = coef(munich4.lm)[2]+const*1.310

p1=ggplot(rent99, aes(x=area, y=rent))+geom_point()+
  geom_smooth(method=lm, se=TRUE)
p2=ggplot(rent99[sample(1:dim(rent99)[1],100),], aes(x=area, y=rent))+geom_point()+
  geom_smooth(method=lm, se=TRUE)
ggarrange(p1,p2)
#geom_line(aes(x=area, y=munich4.beta0.min+munich4.beta1.max*area ),col="red")+
#  geom_line(aes(x=area, y=munich4.beta0.max+munich4.beta1.min*area), col="red")
  
```

The figures show the fitted line and the confidence interval for the (true) regression line based on all (left) and the first 100 observations (right) in the Munich rent index data.

---

## Single hypothesis testing set-up

In single hypothesis testing we are interesting in testing one null hypothesis against an alternative hypothesis. In linear regression the hypothesis is often about a regression parameter $\beta_j$:
$$H_0: \beta_j=0 \text{ vs. } H_1: \beta_j\neq 0$$

<!-- Remark: we implicitly say that our test is done given that the other variables are present in the model, that is, the other $\beta_i$s ($j\neq i$) are not zero. -->

###Two types of errors:

* "Reject $H_0$ when $H_0$ is true"="false positives" = "type I error" ="miscarriage of justice".
These are our _fake news_, which are very important for us to avoid.

* "Fail to reject $H_0$ when $H_1$ is true (and $H_0$ is false)"="false negatives" = "type II error"= "guilty criminal go free".

---

We choose to reject $H_0$ at some significance level $\alpha$ if the $p$-value of the test (see below) is smaller than the chosen significance level.
We say that : Type I error is "controlled" at significance level $\alpha$, which means that the probability of miscarriage of justice (Type I error) does not exceed $\alpha$.

**Q**: Draw a 2 by 2 table showing the connection between 

* "truth" ($H_0$ true or $H_0$ false) - rows in the table, and 
* "action" (reject $H_0$ and accept $H_0$) - columns in the table,

and place the two types of errors in the correct position within the table. 

What else should be written in the last two cells?

---

### Hypothesis test on $\beta_j$ (t-test)

In linear regression models our test statistic for testing $H_0: \beta_j=0$ is 
$$T_0=\frac{\hat{\beta}_j-0}{\sqrt{c_{jj}}\hat{\sigma}_{\varepsilon}}\sim t_{n-2}$$
where $c_{jj}\hat{\sigma}_{\varepsilon}^2=\widehat{\text{Var}}(\hat{\beta}_j)$.

Inserted observed values (and estimates) we have $t_0$.

We would in a two-sided setting reject $H_0$ for large values of $\text{abs}(t_0)$. We may rely on calculating a $p$-value. 

---

### The p-value

A p-value is a test statistic satisfying $0 \leq p({\bf Y}) \leq 1$ for every vector of observations ${\bf Y}$. 

* Small values give evidence that $H_1$ is true. 

* In single hypothesis testing, if the p-value is less than the chosen significance level (chosen upper limit for the probability of committing a type I error), then we reject the null hypothesis, $H_0$. The chosen significance level is often referred to as $\alpha$.

* A p-value is _valid_ if
$$ P(p({\bf Y}) \leq \alpha) \leq \alpha$$
for all $\alpha$, $0 \leq \alpha \leq 1$, whenever $H_0$ is true, that is, if the $p$-value is valid, rejection on the basis of the $p$-value ensures that the probability of type I error does not exceed $\alpha$.

* If $P(p({\bf Y}) \leq \alpha) = \alpha$ for all $\alpha$, $0 \leq \alpha \leq 1$, the $p$-value is called an _exact_ p-value.  

---

In our linear regression we use the $t$-distibution to calculate p-values for our two-sided test situation $H_0: \beta_j=0$ vs. $H_1: \beta_j \neq 0$. Assume we have observed that our test statistic $T_0$ takes the numerical value $t_0$. Since the $t$-distribution is symmetric around $0$ we have

$$p\text{-value}=P(T_0>\text{abs}(t_0))+P(T_0<-\text{abs}(t_0))=2\cdot P(T_0>\text{abs}(t_0)).$$

We reject $H_0$ if our calculated $p$-value is below our chosen signficance level. We often choose as significance level $\alpha=0.05$.

---

**Q:** Comment on the $p$-values listed in the summary output from fitting the simple linear regression of `area` and `rent`. Then, pinpoint $\hat{\sigma}_{\varepsilon}=RSE$. Where are the entries in the output that you do not (yet) know what is?

\footnotesize
```{r,echo=FALSE}
summary(lm(rent~area,data=rent99))
```
\normalsize

---

## Model accuracy

### Coefficient of determination, $R^2$

The total sum of squares is defined as $$\text{TSS} = \sum_{i=1}^n (y_i - \bar{y})^2,$$ and is proportional to the estimated variance of the response variable.

The $R^2$ statisitic is the fraction of variance explained by the model and is given by
$$R^2 = \frac{\text{TSS}-\text{RSS}}{TSS}= 1-\frac{\text{RSS}}{\text{TSS}}=1-\frac{\sum_{i=1}^n(y_i-\hat{y}_i)^2}{\sum_{i=1}^n(y_i-\bar{y}_i)^2}.$$ 
The value is between 0 and 1, and we want an as high $R^2$ statistic as possible. This statistic is independent of the scale of $Y$.

For a simple linear regression model the squared correlation $r^2$ is equal to $R^2$.

---

**Q**: Look back at the `summary` outprint for the `munich1.lm` model. What is the value of the $R^2$ statistic? Based on this value, would you conclude that this model gives a good fit to the data?

---

# Multiple Linear Regression

Back to the Munich rent index data, but now we want to include more than one covariate. Suggestions?

* `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.

---


## Model

We assume, for observation $i$:
$$Y_i= \beta_0 + \beta_{1}  x_{i1} + \beta_2 x_{i2} + ... + \beta_p x_{ip} + \varepsilon_i,$$ where $i=1,2,...,n$. 

Here: what is $x_{ij}$? 

The model can be written in matrix form:
$${\bf Y}={\bf X} \boldsymbol{\beta}+{\boldsymbol{\varepsilon}}.$$ 

$\oplus$: write out in detail!

---

## Notation

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

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

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

${\boldsymbol \varepsilon}: (n\times 1)$ vector of random errors.

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)$.


---


### Example continued

Assume that we have `rent` as response and `area` and `bath` as covariates. 

* What is $\bf{Y}$, $\bf{X}$, $\boldsymbol{\beta}$ and $\boldsymbol{\varepsilon}$?
* Which of these are known/unknown, observed/unobserved?

Remember: $n=3802$ and `bath=0` gives standard quality and `bath=1` premium and the model is
$${\bf Y=X \boldsymbol\beta}+{\boldsymbol \varepsilon}.$$ 

---

\small
```{r}
fit=lm(rent~area+bath,data=rent99)
head(model.matrix(fit))
head(rent99$rent)
```
\normalsize

---

## Classical linear model
Assumptions:

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})=p+1$. (We assume $n>>(p+1)$.)

The classical _normal_ linear regression model is obtained if additionally

4. $\boldsymbol\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 interpretation of the coefficients $\beta_j$ is now as following: holding all other covariates fixed, what is the average effect on $Y$ of a one-unit increase in the $j$th covariate.

---

## Distribution of the response vector
Assume that 
${\bf Y=X \boldsymbol\beta}+{\boldsymbol \varepsilon}$
and 
$\boldsymbol\varepsilon\sim N_n({\bf 0},\sigma^2 {\bf I})$.

**Q:** 

* Find the mean $\text{E}(\bf Y)$ and 
* the covariance matrix $\text{Cov}(\bf Y)$.
* What is then the distribution of $\bf Y$?

**A**:
$$ {\bf Y} \sim N_{n}({\bf X} {\boldsymbol\beta},\sigma^2 {\bf I})$$

---

## Parameter estimation

In multiple linear regression parameters in $\beta$ are estimated with maximum likelihood and least squares. These two methods give the same estimator when we assume the normal linear regression model. 

**Least squares and maximum likelihood estimator for ${\bf \beta}$:**
$$ \hat{\boldsymbol\beta}=({\bf X}^T{\bf X})^{-1} {\bf X}^T {\bf Y}$$

The estimator is found by minimizing the RSS for a multiple linear regression model:
$$\begin{aligned} \text{RSS} &=\sum_{i=1}^n (y_i - \hat y_i)^2 = \sum_{i=1}^n (y_i - \hat \beta_0 - \hat \beta_1 x_{i1} - \hat \beta_2 x_{i2} -...-\hat \beta_p x_{ip} )^2 \\
&= \sum_{i=1}^n (y_i-{\bf x}_i^T \boldsymbol \beta)^2=(\bf{Y}-\bf{X}\hat{\boldsymbol{\beta}})^T(\bf{Y}-\bf{X}\hat{\boldsymbol{\beta}})\end{aligned}$$
The estimator is found by solving the system of (p+1) equation :$\frac{\partial \text{RSS}}{\partial \boldsymbol \beta}={\bf 0}$, see [LeastSquaresMLR.pdf](https://www.math.ntnu.no/emner/TMA4268/2018v/notes/LeastSquaresMLR.pdf) for a derivation (from TMA4267V2017).

---

### Example continued
Write down the model and explain what the values under `Estimate` mean in practice.

\footnotesize
```{r,echo=showsolB}
munich3.lm=lm(rentsqm~area+yearc+location+bath+kitchen+cheating,data=rent99)
summary(munich3.lm)
```
\normalsize

---

Reproduce the values under `Estimate` by calculating without the use of `lm`.

```{r}
X=model.matrix(rentsqm~area+yearc+location+bath+kitchen+cheating,data=rent99)
Y=rent99$rentsqm
betahat=solve(t(X)%*%X)%*%t(X)%*%Y
print(betahat)
```

---

## Distribution of the regression parameter estimator

1. We assumed that 
${\bf Y=X \boldsymbol\beta}+{\boldsymbol \varepsilon}$ and 
$\boldsymbol\varepsilon\sim N_n({\bf 0},\sigma^2 {\bf I})$, leading to 
$$ {\bf Y} \sim N_{n}(\bf X \boldsymbol \beta,\sigma^2 {\bf I}).$$
2. Then we "found" that an estimator for ${\bf \beta}$ is
$$ \hat{\boldsymbol\beta}=({\bf X}^T{\bf X})^{-1} {\bf X}^T {\bf Y}.$$
3. From Module 2: Part B: $\mathbf{Y}_{(n\times 1)}$ with mean vector
    $\mathbf{\mu}$ and variance-covariance matrix $\Sigma$, then $\mathbf{Z}=\mathbf{C}\mathbf{Y}$
has $\text{E}(\mathbf{Z})=\mathbf{C}\mathbf{\mu}$ and $\text{Cov}(\mathbf{Z})=
   \mathbf{C}\Sigma\mathbf{C}^T$.
4. Also Module 2: Part B: If ${\bf Y}$ is multivariate normal, then also $\mathbf{C}\mathbf{Y}$ is multivariate normal. 

**Q: Homework**

* Find the mean $\text{E}(\hat{\boldsymbol\beta})$ and the covariance matrix $\text{Cov}(\hat{\boldsymbol\beta})$.
* What is then the distribution of $\hat{\boldsymbol\beta}$?

---

# <a id="PartB"> Part B: Multiple linear regression - continued!</a>

## Data sets
### Munich rent index (as in Part A)

Munich, 1999: 3082 observations on 9 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.

----

### Ozone

New York, 1973: 111 observations of

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

`ozone` as our response variable and `temperature`, `wind` and `radiation as covariates.

First 6 observations in data set printed.

---

```{r,echo=FALSE,eval=TRUE,messages=FALSE,warnings=FALSE,error=FALSE,results="hold"}
library(ElemStatLearn)
library(knitr)
library(kableExtra)
data(ozne)
kable(head(ozone))
ozone.lm = lm(ozone~temperature+wind+radiation, data=ozone)
```

----

## Plan for part B

* So far: simple linear regression - and multiple linear regression - model, estimators
* Statistical inference 
    + estimators and properties
    + confidence intervals and hypothesis tests
    + significance of regression: F-test
    + prediction and prediction intervals 
* Model assessement and selection
    + $R^2$
    + subset selection 
    + diagnostic plots - studentized residuals and leverages
    
---

* Special cases
    + qualitative predictors: dummy coding 
    + non-additivity: including interactions
    + nonlinearity: polynomials
    + multiple linear regression vs. KNN-regression 
* Optional material
    + projection matrices and geometry of least squares

* Summing-up with team Kahoot!

---

# Multiple linear regression (continued)

## Model and parameter estimator (summing up)

Model: $${\bf Y}=\bf X \boldsymbol{\beta}+{\boldsymbol{\varepsilon}}$$
with full rank design matrix ${\bf X}$. Classical _normal_ linear regression model when
$$\boldsymbol{\varepsilon}\sim N_n({\bf 0},\sigma^2{\bf I}).$$

In Part A we found that :
$$ {\bf Y} \sim N_{n}({\bf X} {\boldsymbol\beta},\sigma^2 {\bf I})$$

Parameter of interest is $\boldsymbol{\beta}$ and $\sigma^2$ is a nuisance (=parameter not of interest). 

Using the least squares (and maximum likelihood) method the estimator for $\boldsymbol\beta$ is
$$ \hat\beta=({\bf X}^T{\bf X})^{-1} {\bf X}^T {\bf Y}$$

---

\small
```{r}
head(model.matrix(ozone.lm))
head(ozone$ozone)
```
\normalsize

---

How does this compare to simple linear regression?
$$\hat{\beta}_0 = \bar{Y}-\hat{\beta}_1 \bar{x} \text{ and } \hat{\beta}_1 = \frac{\sum_{i=1}^n(x_i-\bar{x})(Y_i-\bar{Y})}{\sum_{i=1}^n(x_i-\bar{x})^2},$$

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

Often we use centered data (and also scaled) to ease interpretation. In design of experiments often orthogonal columns of the design matrix is chosen to get ${\bf X}^T{\bf X}$ to be diagonal, which leads to easier interpretation and identifiability.

---

## Distribution of the regression parameter estimator (summing up)
$$ \hat{\boldsymbol\beta}=({\bf X}^T{\bf X})^{-1} {\bf X}^T {\bf Y}$$

This can be written as $\hat{\boldsymbol\beta}={\bf C}{\bf Y}$ where 

* ${\bf C}=({\bf X}^T{\bf X})^{-1} {\bf X}^T$
* ${\bf Y} \sim N_{n}({\bf X} {\boldsymbol\beta},\sigma^2 {\bf I})$. 

Therefore

* $\hat{\boldsymbol\beta}$ is multivariate normal (p+1) dimensions, with
* $\text{E}(\hat{\boldsymbol\beta})={\bf C}\text{E}({\bf Y})=({\bf X}^T{\bf X})^{-1} {\bf X}^T{\bf X} {\boldsymbol\beta}={\boldsymbol\beta}$
* $\text{Cov}(\hat{\boldsymbol\beta})={\bf C}\text{Cov}({\bf Y}){\bf C}^T=({\bf X}^T{\bf X})^{-1} {\bf X}^T \sigma^2 {\bf I}(({\bf X}^T{\bf X})^{-1} {\bf X}^T)^T =({\bf X}^T{\bf X})^{-1} \sigma^2$.

So: $\hat{\beta}\sim N_{p+1}(\boldsymbol{\beta},\sigma^2({\bf X}^T{\bf X})^{-1})$.

<!-- [PropertiesBetahatMLR.pdf](https://www.math.ntnu.no/emner/TMA4268/2018v/PropertiesBetahatMLR.pdf) for a derivation (from TMA4267V2017) -->

---

$\hat{\beta}\sim N_{p+1}(\boldsymbol{\beta},\sigma^2({\bf X}^T{\bf X})^{-1})$.

\footnotesize
```{r}
coefficients(ozone.lm)
vcov(ozone.lm)
```
\normalsize

**Q:** Explain what all these numbers are!

---

## Estimator for $\sigma^2$ 

$$\begin{aligned} \text{RSS} &=\sum_{i=1}^n (y_i - \hat y_i)^2 = \sum_{i=1}^n (y_i - \hat \beta_0 - \hat \beta_1 x_{i1} - \hat \beta_2 x_{i2} -...-\hat \beta_p x_{ip} )^2 \\
&= \sum_{i=1}^n (y_i-{\bf x}_i^T \hat{\boldsymbol \beta})^2=(\bf{Y}-\bf{X}\hat{\boldsymbol{\beta}})^T(\bf{Y}-\bf{X}\hat{\boldsymbol{\beta}})\end{aligned}$$

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

---

## Distribution of regression parameters (contd.)

$\hat{\beta}\sim N_{p+1}(\boldsymbol{\beta},\sigma^2({\bf X}^T{\bf X})^{-1})$.

* unbiased
* covariance matrix dependent on the design (and $\sigma^2$)

Multicollinearity: columns in design matrix (that is, the covariates) are correlated, which may lead to difficulty in "identifying" the effect of each covariate on the response, and thus large variances (and covariances) for the elements of $\hat{\boldsymbol\beta}$.

The _variance inflation factor (VIF)_ is the ratio of the variance of $\hat{\beta}_j$ when fitting a model with the chosen covariates divided by the variance of $\hat{\beta}_j$ in a simple linear regression. 

* VIF=1: absence of collinearity
* VIF exceeding 5 or 10 might be problematic. 
* Solution: drop a covariate (that do not add much since it is correlated with other covariates).

---

\small
```{r,echo=c(5,6,7,8),out.width="45%"}
library(corrplot)
library(gamlss.data)
library(ElemStatLearn)
library(car)
oz1=as.data.frame(apply(ozone,2,scale,scale=FALSE))
fitoz=lm(ozone~temperature+wind+radiation, data=oz1)
vif(fitoz)
corrplot(cov2cor(vcov(fitoz)),cex.lab=0.7);
```
\normalsize

---

## Inference about $\beta_j$: confidence interval
$$\hat{\beta}\sim N_{p+1}(\boldsymbol{\beta},\sigma^2({\bf X}^T{\bf X})^{-1})$$

Statistic for inference about $\beta_j$, $c_{jj}$ is diagonal element corresponding to $\hat{\beta}_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-1}$$
$$ P(\hat{\beta}_j-t_{\alpha/2,n-p-1}\sqrt{c_{jj}}\hat{\sigma}
\le \beta_j \le \hat{\beta}_j+t_{\alpha/2,n-p-1}\sqrt{c_{jj}}\hat{\sigma})=1-\alpha$$

A $(1-\alpha)$% CI for $\beta_j$ is when we insert numerical values for the upper and lower limits: $[\hat{\beta}_j-t_{\alpha/2,n-p-1}\sqrt{c_{jj}}\hat{\sigma},\hat{\beta}_j+t_{\alpha/2,n-p-1}\sqrt{c_{jj}}\hat{\sigma}]$.

When we work with _large samples_ then $n-p-1$ 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-1}$. This is in done in our textbook.

---

\footnotesize
```{r}
fitoz=lm(ozone~temperature+wind+radiation, data=ozone)
confint(fitoz)
```
\normalsize

---

## Inference about $\beta_j$: hypothesis tests
$$H_0: \beta_j=0 \text{ vs } H_1: \beta_j\neq 0$$

Nothing new: using $t$-tests based on
$$ T_j=\frac{\hat{\beta}_j-\beta_j}{\sqrt{c_{jj}}\hat{\sigma}}\sim t_{n-p-1}$$

Again, $c_{jj}$ is diagonal element corresponding to $\hat{\beta}_j$ of $({\bf X}^T{\bf X})^{-1}$.

---

\footnotesize
```{r}
fitoz=lm(ozone~temperature+wind+radiation, data=ozone); summary(fitoz)
```
\normalsize

---

**Q:**

1. Where is hypothesis testing performed here, and which are the hypotheses rejected at level $0.01$?
2. Will the test statistics and $p$-values change if we change the regression model?
3. What is the relationship between performing an hypothesis test and constructing a CI interval?

** A:**
1. Column named `t value` gives numerical value for the t-statistic and column `*Pr(>|t|)` gives $p$-value for the two-sided test.
2. Yes. Unless our design matrix has orthogonal columns.
3. If a $(1-\alpha)$ 100% CI covers the hypothsized value (here 0) then this is a sensible value for the parameter it can be shown that then the $p$-value will be larger than $\alpha$. And, if a $(1-\alpha)$ 100% CI does not cover the hypothsized value (here 0) then this is not a sensible value for the parameter it can be shown that then the $p$-value will be smaller than $\alpha$.

---

## More complex hypotheses

Consider the hypotheses:
$$ H_0: \beta_1=\beta_2=\cdots= \beta_k =0 \text{ vs. } H_1: \text{at least one different from zero}.$$

This means we test a set of regression parameters is different from 0. 

This is used to

* let $k=p$ and thest if any parameters are different from 0 - this is called to _test if the regression is significant_.
* to compare two models - one where a the subset of the coefficients are omitted.

To do this $F$-tests are used.

---

### Is the regression significant?

$$ H_0: \beta_1=\beta_2=\cdots= \beta_p =0 \text{ vs. } H_1: \text{at least one different from zero}.$$

F-statistic

$$F=\frac{(\text{TSS-RSS})/p}{RSS/(n-p-1)} \sim F_{p,n-p-1}$$

When $H_0$ is false we expect that the numerator is larger than the denominator and thus $F$ is greater than 1. A $p$ value is calculated from the upper tail of the $F$-distribution

We find $\text{p-value}=P(F_{p,n-p-1}>f_0)$, where $f_0$ is the numerical value of $F$ inserted RSS, TSS from the data.

---

\footnotesize
```{r}
fitoz=lm(ozone~temperature+wind+radiation, data=ozone); summary(fitoz)
```
\normalsize

---

### Tested a subset by comparing two nested models

* Large model: RSS with $p+1$ regression parameters
* Small model: RSS$_0$ with $q+1$ regression parameters

$H_0:$ the additional $p-q$ coefficients in the large model are all zero vs. $H_1$: at least one different from zero.


$$F=\frac{(\text{RSS$_0$-RSS})/(p-q)}{RSS/(n-p-1)} \sim F_{p-q,n-p-1}$$

When $H_0$ is false we expect that the numerator is larger than the denominator and thus $F$ is greater than 1. A $p$ value is calculated from the upper tail of the $F$-distribution

We find $\text{p-value}=P(F_{p-q,n-p-1}>f_0)$, where $f_0$ is the numerical value of $F$ inserted RSS, TSS from the data.

In R we perform the test by fitting the two models `fit.large` and `fit.small` and use `anova(fit.small,fit.large).

---

\footnotesize
```{r}
fit.large=lm(ozone~temperature+wind+radiation, data=ozone)
fit.small=lm(ozone~temperature, data=ozone)
anova(fit.small,fit.large)
```
\normalsize

---

## Prediction intervals 

Once we have estimated the coeffients $\hat\beta_0$, $\hat\beta_1$,.., $\hat\beta_p$, we can make a prediction for a response value $Y_0$ for a new observation $\mathbf x_0 = (x_1, x_2, ..., x_p)$ as before:

$$\hat Y_0 = \hat \beta_0 + \hat \beta_1 x_1 + \hat \beta_2 x_2 + ... + \hat \beta_p x_p ={\bf x}_0^T \hat{\boldsymbol\beta}.$$
This is an intuitive point estimate.

Remember, one aim for regression was to "construct a model to predict the reponse from a set of (one or several) explanatory variables- more or less black box".

To assess the uncertainty in this prediction we present a prediction interval for the $Y_0$.

----

After some work (see "details on the derivation"):

$P({\bf x}_0^T \hat{\beta}-t_{\alpha/2,n-p-1}\hat{\sigma}\sqrt{1+{\bf x}_0^T({\bf X}^T{\bf X})^{-1}{\bf x}_0} \le Y_0 \le {\bf x}_0^T \hat{\beta}+t_{\alpha/2,n-p-1}\hat{\sigma}\sqrt{1+{\bf x}_0^T({\bf X}^T{\bf X})^{-1}{\bf x}_0})=1-\alpha$

A $(1-\alpha)$% PI for $Y_0$ is when we insert numerical values for the upper and lower limits: $[{\bf x}_0^T \hat{\beta}-t_{\alpha/2,n-p-1}\hat{\sigma}\sqrt{1+{\bf x}_0^T({\bf X}^T{\bf X})^{-1}{\bf x}_0}, {\bf x}_0^T \hat{\beta}+t_{\alpha/2,n-p-1}\hat{\sigma}\sqrt{1+{\bf x}_0^T({\bf X}^T{\bf X})^{-1}{\bf x}_0}]$.

---

PIs can be found in R using `predict` on an `lm` object, but make sure that `newdata` is a `data.frame` with the same names as the original data. 

**Example: Using the Munich rent index data**

We want to predict the rent - with PI - for an appartment with area 50, location 2 ("good"), nice bath and kitchen and with central heating.


---

\footnotesize
```{r}
require(gamlss.data)
fit=lm(rent~area+location+bath+kitchen+cheating,data=rent99)
newobs=rent99[1,]
newobs[1,]=c(NA,NA,50,NA,2,1,1,1,NA)
predict(fit,newdata=newobs,interval="prediction",type="response")
```
\normalsize

**Q** 

1. When is a prediction interval of interest? 

**A:** Always? Gives useful information in addition to a point prediction.

2. Explain the result from `predict` above.

**A:** Fit is point prediction, lwr is lower and upr is upper limit of the 95% PI (default with 95, argument if other.)

---

### Details in the derivation of the PI

We start to look at the difference between the unobserved response $Y_0$ (for a given covariate vector ${\bf x}_0$) and the point prediction $\hat{Y}_0$, $Y_0-\hat{Y}_0$.

First, we assume that the unobserved response at covariate ${\bf x}_0$ is independent of our previous observations and follows the same distibution, that is
$Y_0 \sim N({\bf x}_0^T \beta,\sigma^2)$. Further,

$$\hat{Y}_0={\bf x}_0^T \hat{\beta} \sim N({\bf x}_0^T \beta,\sigma^2 {\bf x}_0^T({\bf X}^T{\bf X})^{-1}{\bf x}_0).$$

Then, for $Y_0-{\bf x}_0^T \hat{\beta}$ we have

$\text{E}(Y_0-{\bf x}_0^T \hat{\beta})=0 \text{ and } \text{Var}(Y_0-{\bf x}_0^T \hat{\beta})=\text{Var}(Y_0)+\text{Var}({\bf x}_0^T \hat{\beta})=\sigma^2+\sigma^2{\bf x}_0^T({\bf X}^T{\bf X})^{-1}{\bf x}_0$

so that
$$Y_0-{\bf x}_0^T \hat{\beta}\sim N(0,\sigma^2 (1+{\bf x}_0^T({\bf X}^T{\bf X})^{-1}{\bf x}_0)) $$

---

Inserting our REML-estimate for $\sigma^2$ gives

$$T=\frac{Y_0-{\bf x}_0^T \hat{\beta}}{\hat{\sigma}\sqrt{1+{\bf x}_0^T({\bf X}^T{\bf X})^{-1}{\bf x}_0}}\sim t_{n-p-1}.$$

Then, we start with
$$ P(-t_{\alpha/2,n-p-1}\le \frac{Y_0-{\bf x}_0^T \hat{\beta}}{\hat{\sigma}\sqrt{1+{\bf x}_0^T({\bf X}^T{\bf X})^{-1}{\bf x}_0}} \le t_{\alpha/2,n-p-1})=1-\alpha$$
and solve so that $Y_0$ is in the middle, which gives

$P({\bf x}_0^T \hat{\beta}-t_{\alpha/2,n-p-1}\hat{\sigma}\sqrt{1+{\bf x}_0^T({\bf X}^T{\bf X})^{-1}{\bf x}_0} \le Y_0 \le {\bf x}_0^T \hat{\beta}+t_{\alpha/2,n-p-1}\hat{\sigma}\sqrt{1+{\bf x}_0^T({\bf X}^T{\bf X})^{-1}{\bf x}_0})=1-\alpha$

---

## Coefficient of determination, $R^2$

no change from simple linear regression.

$$R^2 = \frac{\text{TSS}-\text{RSS}}{TSS}= 1-\frac{\text{RSS}}{\text{TSS}}=1-\frac{\sum_{i=1}^n(y_i-\hat{y}_i)^2}{\sum_{i=1}^n(y_i-\bar{y}_i)^2}.$$ 

1. The interpretation of this coefficient is that the closer it is to 1 the better the fit to the data. If $R^2=1$ then all residuals are zero - that is, perfect fit to the data.

2. In a simple linear regression the $R^2$ equals the squared correlation coefficient between the reponse and the predictor. In multiple linear regression $R^2$ is the squared correlation coefficient between the observed and predicted response.

3. If we have two models M1 and M2, where model M2 is a submodel of model M1, then
$$ R^2_{M_1}\ge R^2_{M_2}.$$
This can be explained from the fact that $\text{RSS}_{M_1}\le \text{RSS}_{M_2}$.
(More in the Recommended exercises.)

## Model assessment and selection

### Quality measures

To assess the quality of the regression we can report the $R^2$ coefficient of determination. However, since adding covariates to the linear regression can not make the RSS larger, this means that adding covariates can not make the $R^2$ smaller. This means that RSS and $R^2$ are only useful measures for comparing models with the same number of regression parameters estimated.

If we consider two models with the same model complexity then RSS can be used to choose between (or compare) these models.

But, if we want to compare models with different model complexity we need to look at other measures of quality for the regression.

---

### $R^2$ adjusted 

$$R^2_{\text{adj}}=1-\frac{\frac{RSS}{n-p-1}}{\frac{TSS}{n-1}}=1-\frac{n-1}{n-p-1}(1-R^2)$$
Choose the model with the _largest_ $R^2_{\text{adj}}$.

---

### AIC Akaike information criterion

AIC is one of the most widely used criteria, and is designed for likelihood-based inference.
Let $l(\hat{\beta}_M,\tilde{\sigma}^2)$ be the maximum of the log-likelihood of the data inserted the maximum likelihood estimates for the regression parameters. Further, let $\lvert M \rvert$ be the number of estimated regression parameters in our model (and then we add 1 for the estimated variance).
$$\text{AIC} =-2 \cdot l(\hat{\beta}_M,\tilde{\sigma}^2)+2(\lvert M\rvert +1)$$

For a normal regression model:
$$\text{AIC} \propto n\ln(\tilde{\sigma}^2)+2(\lvert M\rvert +1)+C$$
where C is a function of $n$ (will be the same for two models for the same data set).
Remark that $\tilde{\sigma}^2=RSS/n$ - our ML estimator, so that the first term in the AIC is just a function of the RSS.

Choose the model with the minimum AIC.

Remark: for us here we have $2(\lvert M\rvert +1)=2(p+2)$

---

The likelihood - for the normal linear regression model:
$$l(\hat{\boldsymbol\beta},{\tilde{\sigma}^2})=\text{ln}(L(\hat{\boldsymbol\beta},{\sigma^2}))=-\frac{n}{2}\text{ln}(2\pi)-\frac{n}{2}\text{ln}\tilde{\sigma}^2-\frac{1}{2\tilde{\sigma}^2} ({\bf y}-{\bf X}\hat{\boldsymbol\beta})^T({\bf y}-
{\bf X}\hat{\boldsymbol\beta})$$

---

###BIC  Bayesian information criterion.

The BIC is also based on the likelihood (see notation above).
$$\text{BIC} =-2 \cdot l(\hat{\beta}_M,\tilde{\sigma}^2)+\ln(n)\cdot (\lvert M\rvert +1)$$

For a normal regression model:
$$ \text{BIC} \propto n\ln(\tilde{\sigma}^2)+\ln(n)(\lvert M\rvert +1)$$
Choose the model with the minimum BIC.

AIC and BIC are motivated in very different ways, but the final result for the normal regression model is very similar.
BIC has a larger penalty than AIC ($\log(n)$ vs. $2$), and will often give a smaller model (=more parsimonious models) than AIC.
In general we would not like a model that is too complex.

Remark: for us here we have $\ln(n)\cdot (\lvert M\rvert +1)=\ln(n)\cdot (p+2)$

---

## Subset selection procedures

Given a full multiple linear regression model the case is often that not all of the covariates are of equal importance for predicting the response. Some of the covariates could be removed from the model without affecting the fit. At the same time one would gain a more interpretable model with fewer parameters. We will here shortly discuss three subset selction procedures. These will be discussed in greater detail in Module 6.

**Forward selection**
The forward selection procedure stars with the null model (no covariates, only $\beta_0$). In a step-wise procedure, additional covariates are added one at a time. The inclusion of the covariate is based on a quality of fit measure (f.ex. $R^2_{\text{adj}}$ or AIC) and the covariate giving the greatest improve of the fit is added to the model at each step.

---

**Backward selection**
The backward selection proceudre has the opposite procedure: Here the starting point is the full multiple linear regression model (with all covariates included). At each step of the procedure, one covariate is removed from the model. The removal of this covariate is based on a quality of fit measure, where the covariate corresponding to the smallest decrease in the fit is removed at each step.

**Mixed selection**
Combine forward and backward.

**All subset selection**
All subset selection is a procedure where all possible combination of covariates are tested. This approach is efficient but computationally expensive.

Model selection is a major point in Module 6.

---

## Challenges - for model fit

1. Non-linearity of data
2. Correlation of error terms
3. Non-constant variance of error terms
4. Normality of error terms
5. Outliers
6. High leverage points
7. Collinearity

---

## Diagnostic plots

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

1. Plot the residuals 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, 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. As an additional aid a test for normality can be used, but must be interpreted with caution since for small sample sizes the test is not very powerful and for large sample sizes even very small deviances from normality will be labelled as significant.
4. Plot the residuals 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_.

---

### Different types of residuals

If can be shown that the vector of residuals, $\bf e=(e_1,e_2,\ldots,e_n)$ have a normal (singular) distribution with mean $\text{E}(\bf e)={\bf 0}$ and covariance matrix $\text{Cov}(\bf e)=\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.

**Q:** 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$ so they can not be used for this. However, we have made "predictions" of the errors - our residuals. And, we want to use our residuals to check the model assumptions. 

That is, we want to check that our errors are independent, homoscedastic (same variance for each observation), and not dependent on our covariates - and we want to use the residuals (observed) in place of the errors (unobserved). Then it would have been great if the residuals have these properties when the underlying errors have. To amend our problem we need to try to fix the residual so that they at least have equal variances. We do that by working with _standardized_ or _studentized residuals_.

---

**Standardized residuals:**

$$r_i=\frac{e_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{e_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.

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

---

### Diagnostic plots in `R`

More information on the plots here: <http://data.library.virginia.edu/diagnostic-plots/> and <http://ggplot2.tidyverse.org/reference/fortify.lm.html>

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

You can use the function `fortify.lm` in `ggplot2` to create a dataframe from an `lm`-object, which `ggplot` uses automatically when given a `lm`-object. This can be used to plot diagnostic plots.

For simplicity we use the Munch rent index with `rent` as response and only `area` as the only covariate. (You may change the model to a more complex one, and rerun the code chunks.)

\footnotesize
```{r tidy = TRUE, echo=showsolB}
fit <- lm(rent ~ area, data=rent99) # Run a regression analysis
format(head(fortify(fit)), digits = 4L)
# to show what fortify implicitly does in ggplot
# for more information ggplot2::fortify.lm
```

---

#### 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=showsolB}
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=showsolB}
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))
library(nortest)
ad.test(rstudent(fit))
```

<!-- (Not normal.) -->

----

#### 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}
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; even though data have extreme values, they might not be influential to determine the regression line (the results don't differ much if they are removed from the data set). These influential outliers can be seen as observations that does not get along with the trend in the majority of the observations. In `plot.lm`, dashed lines are used to indicate the Cook's distance, instead of using the size of the dots as is done here.

---

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=showsolB}
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.) -->

---

# Extensions and challenges in multiple regression

## Qualitative covariates
See [Rob Tibshirani explain - from ca 9 minutes](https://www.youtube.com/watch?v=3T6RXmIHbJ4&index=4&list=PL5-da3qGB5IBSSCPANhTgrw82ws7w_or9)

Qualitative predictors can be included in a linear regression model by introducing dummy variables

**Example**: consider our `rent` dataset with `rent` as reponse, and continuous covariate `area` and categorical covariate `location`. Let the `location` be a factor with levels `average, good, top`.

\footnotesize

```{r,results="hide"}
require(gamlss.data)
require(tidyverse)
require(GGally)

ds=rent99 %>%
  select(location, area,rent)
levels(ds$location)
# change to meaningful names
levels(ds$location)=c("average","good","top")
ggpairs(ds)
```
\normalsize

**Q**: comment on what you see in the `ggpairs` plot.

----

Categorical covariates may either be ordered or unordered. We will only consider unordered categories here.
In general, we could like to estimate regression coefficients for all levels for the categorical covariates. However, if we want to include an intercept in our model we can only include codings for one less variable than the number of levels we have - or else our design matrix will not have full rank.

**Q**: Assume you have a categorical variable with three levels. Check for yourself that making a design matrix with one intercept and three columns with dummy (0-1) variable coding will result in a matrix that is singular.

\footnotesize
```{r}
# make "wrong" dummy variable coding with 3 columns
n=length(ds$location)
X=cbind(rep(1,n),ds$area,rep(0,n),rep(0,n),rep(0,n))
X[ds$location=="average",3]=1
X[ds$location=="good",4]=1
X[ds$location=="top",5]=1
X[c(1,3,69),]
require(Matrix)
dim(X)
rankMatrix(X)
```
\normalsize

----

This is why we need to instead work with different ways of coding categorical variables. One solution is to not include an intercept in the model, but that is often not what we want. We will look at two other solutions - one where we decide on a reference category (that we not include in the coding, and therefore is kind of included in the intercept - this is called "treatment coding") and one where we require that the the sum of the coeffisients are zero (called "effect coding"). This mainly effects how we interpret parameter estimates and communicate our findings to the world. We will here restrict our discussion to "treatment coding".

---

If we fit a regression model with `lm` to the data with `rent` as response and `area` and `location` as covariates, a model matrix is made - and how to handle the categorical variable is either specified the call to `lm` in `contrasts=list(location="contr.treatment")` (or to model.matrix) or globally for all categorical variables with `options(contrasts=c("contr.treatment","contr.poly"))`- where first element give choice for unordered factor (then treatment contrast is default) and second for ordered (and then this polynomial contrast is default). We will only work with unordered factors now.

---

### Dummy variable coding

This is the default coding. The reference level is automatically chosen as the "lowest" level (sorted alphabetically). For our example this means that the reference category for location is "average".

$$x_{i \text{locationgood}} = \begin{cases} 1 \text{ if } i \text{ -th location}=\text{"good"} \\ 0 \text{ if } i \text{ -th location }\neq\text{ "good"} \end{cases}$$
$$x_{i \text{locationtop}} = \begin{cases} 1 \text{ if } i \text{ -th location}=\text{"top"} \\ 0 \text{ if } i \text{ -th location }\neq\text{ "top"} \end{cases}$$
$$\begin{aligned} y_i &= \beta_0 + \beta_1 x_{i \text{area}} + \beta_2 x_{i \text{locationgood}}+\beta_3 x_{i \text{locationtop}} + \varepsilon_i\\\\ &= \begin{cases} \beta_0 + \beta_1 x_{i \text{area}}+ \beta_2 + \varepsilon_i  &\text{ if } i \text{ -th location}=\text{"good"}  \\ \beta_0 + \beta_1 x_{i \text{area}}+ \beta_3 + \varepsilon_i &\text{ if } i \text{ -th location}=\text{"top"} \\ \beta_0 + \varepsilon_i &\text{ if } i \text{ -th location}=\text{"average"}\end{cases}\end{aligned}$$

If we instead wanted "good" to be reference category we could relevel the factor.

---

```{r}
X1=model.matrix(~area+location,data=ds)
X1[c(1,3,69),]
ds$locationRELEVEL=relevel(ds$location,ref="good")
X2=model.matrix(~area+locationRELEVEL,data=ds)
X2[c(1,3,69),]
```

So, what does this mean in practice? Model 1 has `average` as reference category and model 2 `good`.

```{r}
fit1=lm(rent~area+location,data=ds,contrasts = list(location="contr.treatment"))
summary(fit1)
fit2=lm(rent~area+locationRELEVEL,data=ds,contrasts = list(locationRELEVEL="contr.treatment"))
summary(fit2)
```

**Q**: Comment on the print-out. How do we interpret the intercept estimate?

---

## Interactions
See [Trevor Hastie explain](https://www.youtube.com/watch?v=IFzVxLv0TKQ&list=PL5-da3qGB5IBSSCPANhTgrw82ws7w_or9&index=5) 

To illustrate how interactions between covaraites 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)

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"}
library(ElemStatLearn)
library(knitr)
library(kableExtra)
data(ozne)
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:** 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 higly 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.

There reasons for this is that a model with interaction terms, but without the main effects is hard to interpret.

---

## 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}$$

---

\footnotesize
```{r,echo=showsolB}
temp.cat = ifelse(ozone$temperature < mean(ozone$temperature), "low", "high")
ozone2 = cbind(ozone, temp.cat)
kable(head(ozone2))
ozone.int2 = lm(ozone~wind + temp.cat+ temp.cat*wind, data=ozone2)
summary(ozone.int2)
```
\normalsize

---

```{r,echo=showsolB}
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")
```

---

## Nonlinearity

We may extend the linear model to handle non-linear relationships by using polynomial regression - as we did in Module 2 in our bias-variance trade-off example.

More on non-linearity in Module 7.

---

<!-- # Multiple linear regression vs. KNN-regression -->

<!-- read pages 104-109 in the textbook. -->

<!-- --- -->

# Optional: Projection matrices

First, we define predictions as $\hat{\bf Y}={\bf X}\hat{\beta}$, and inserted the ML (and LS) estimate we get $\hat{\bf Y}={\bf X}({\bf X}^T{\bf X})^{-1}{\bf X}^T{\bf Y}$.

We define the projection matrix
$$  {\bf H}={\bf X}({\bf X}^T{\bf X})^{-1} {\bf X}^T$$
called the _hat matrix_. This simplifies the notation for the predictions,
$$\hat{\bf Y}={\bf H}{\bf Y}$$ so the hat matrix is putting the hat on the response ${\bf Y}$.

In addition we define residuals as
\begin{align*}
\hat{\varepsilon}&={\bf Y}-\hat{\bf Y} \\
\hat{\varepsilon}&={\bf Y}-{\bf HY}=({\bf I-H}){\bf Y}\\
\end{align*}
so we have a second projection matrix
$$ {\bf I-H}={\bf I}-{\bf X}({\bf X}^T{\bf X})^{-1} {\bf X}^T $$

---

## Optional: Geometry of Least Squares (involving our two projection matrices)

* Mean response vector: $\text{E}({\bf Y})={\bf X}{\boldsymbol \beta}$
* As $\beta$ varies, ${\bf X}\beta$ spans the
    model plane of all linear combinations. I.e. the space
    spanned by the columns of ${\bf X}$: the column-space of
    ${\bf X}$.
* Due to random error (and unobserved covariates),
        ${\bf Y}$ is not exactly a linear combination of the
        columns of ${\bf X}$.
* LS-estimation chooses $\hat{\beta}$ such that
        ${\bf X}\hat{\beta}$ is  the point in
        the column-space of ${\bf X}$ that is closes to ${\bf Y}$.
* The residual vector
        $\hat{\varepsilon}={\bf Y}-\hat{{\bf Y}}=({\bf I-H}){\bf Y}$ is perpendicular
        to the column-space of ${\bf X}$.
* Multiplication by
        ${\bf H}={\bf X}({\bf X}^T{\bf X})^{-1}{\bf X}^T$ projects a
        vector onto the column-space of ${\bf X}$.
* Multiplication by
        ${\bf I-H}={\bf I}-{\bf X}({\bf X}^T{\bf X})^{-1}{\bf X}^T$ projects a
        vector onto the space perpendicular to the column-space of
        ${\bf X}$.

---

<!-- # Sadly not covered this year (but maybe next year - or you check for yourself) -->

<!-- * Big data - `big.lm` -->
<!-- * Fitting a linear model with with gradient descent - to motivate neural net estimation methods. -->

---

# Important results in multiple linear regression 

* Linear regression assumes a linear relationship between the response variable and the covariates.
* Simple linear regression has only one covariate and has the form $Y = \beta_0 + \beta_1 X + \varepsilon$.
* Muliple linear regression has $p$ covariates and has the form $Y = \beta_0 + \beta_1 X_1 + \beta_2 X_2 + ... + \beta_p X_p + \varepsilon$.
* Quantitative (categorical) covariates can be included using dummy variables.
* Correlations among the covariates can mask each others effects in the linear model.

* Parameter estimates can be obtained by minimizing the least squares (RSS) or maximum likelihood.
* We can calculate the standard errors of the parameter estimates, and use this to obtain confidence intervals.
* We can test the hypothesis of $H_0 : \beta_j = 0$ agains $H_1: \beta_j \neq 0$ by a t-test.
* We use the $F$-statistic to test if at least one of the covariates are useful.

---

* Not only additive effects: Interactions between covariates can be included in the model (also between qualitative and quantitative covariates).
* Transformations of the response variable or of a covariate can be useful if the relationship is not linear. A linear model can then be fit to the transformed variables.

* The overall accuracy of the model can be evaluated by calculating the $R^2$ statistic, AIC score and by using diagnostic plots.
* Model selection can not be based on RSS or $R^2$.
* Multiple linear regression might require subset selection if the number of covariates is high.

<!-- --- -->

<!-- # Summing up with Kahoot! -->

<!-- --- -->

# <a id="RecEx">Recommended exercises</a>

## Problem 1: Theoretical questions

###a)
A core finding is $\hat{\boldsymbol\beta}$.
\[ \hat{\boldsymbol\beta}=({\bf X}^T{\bf X})^{-1} {\bf X}^T {\bf Y}\]
with $\hat{\boldsymbol\beta}\sim N_{p}(\boldsymbol\beta,\sigma^2({\bf X}^T{\bf X})^{-1})$.

* Show that $\hat{\boldsymbol\beta}$ has this distribution with the given mean and covariance matrix. 
* What do you need to assume to get to this result? 
* What does this imply for the distribution of the $j$th element of $\hat{\beta}$? 
* In particular, how can we calculate the variance of $\hat{\beta}_j$?

###b) 
What is the interpretation of a 95% confidence interval? Hint: repeat experiment (on $Y$), on average how many CIs cover the true $\beta_j$?

###c)
What is the interpretation of a 95% prediction interval? Hint: repeat experiment (on $Y$) for a given ${\bf x}_0$.

###d)
Construct a 95% CI for ${\bf x}_0^T \beta$. Explain what is the connections between a CI for $\beta_j$, a CI for ${\bf x}_0^T \beta$ and a PI for $Y$ at ${\bf x}_0$.

###e)
Explain the difference between _error_ and _residual_.  What are the properties of the raw residuals? Why don't we want to use the raw residuals for model check? What is our solution to this?

###f)
Consider a multiple linear regression model $A$ and a submodel $B$ (all parameters in $B$ are in $A$ also). We say that $B$ is nested within $A$. Assume that regression parameters are estimated using least squares. Why is then the following true: RSS for model $A$ will always be smaller or equal to RSS for model $B$. And thus, $R^2$ for model $A$ can never be worse than $R^2$ for model B. (See also Problem 2d below.)

----

## Problem 2: Munich Rent index

###a)
Fit the regression model with first `rent` and then `rentsqm` as reponse and following covariates: `area`, `location` (dummy variable coding using location2 and location3, just write `as.factor(location)`), `bath`, `kitchen` and `cheating` (central heating).

Look at diagnostic plots for the two fits. Which response do you prefer?

Consentrate on the response-model you choose for the rest of the tasks.

###b) 
Explain what the parameter estimates mean in practice. In particular, what is the interpretation of the intercept?

###c) 
Go through the summary printout and explain all parts.

###d) 
Now we add random noise as a covariance, but simulating the IQ of the landlord of each appartment. Observe that $R^2$ increases (or stays unchanged) and RSS decreases (or stays the same) if we add IQ as covariate, but $R^2_{\text{adj}}$ decreases. What does this tell you about model selection and overfitting?

For the code - what is the connection between `sigma` and RSS?

```{r,results="hold"}
orgfit=lm(rent~area+as.factor(location)+bath+kitchen+cheating,data=rent99)
summary(orgfit)
set.seed(1) #to be able to reproduce results
n=dim(rent99)[1]
IQ=rnorm(n,100,16)
fitIQ=lm(rent~area+as.factor(location)+bath+kitchen+cheating+IQ,data=rent99)
summary(fitIQ)

summary(orgfit)$sigma
summary(fitIQ)$sigma

summary(orgfit)$r.squared
summary(fitIQ)$r.squared
summary(orgfit)$adj.r.squared
summary(fitIQ)$adj.r.squared
```


###e) 
We now want to use model selection to arrive at a good model. Start by defining which covariates you want to include and how to code them (use dummy variable coding of `location`). What about year of construction - is that a linear covariate? Maybe you want to make intervals in time instead? Linear or categorical for the time? What about the `district`? We leave that since we have not talked about how to use spatial covariates.

Hint: if you want to test out interval versions of year of construction the function `mutate` (from `dplyr`) is useful:
```{r, eval=FALSE}
rent99 <- rent99 %>% mutate(yearc.cat = cut(yearc, breaks = c(-Inf, seq(1920,2000,10)), labels = 10*1:9))
```
More on `dplyr`: Tutorial: <http://genomicsclass.github.io/book/pages/dplyr_tutorial.html>
and Cheat sheet (data wrangling): <https://www.rstudio.com/wp-content/uploads/2015/02/data-wrangling-cheatsheet.pdf> and dplyr in particular: <https://github.com/rstudio/cheatsheets/raw/master/source/pdfs/data-transformation-cheatsheet.pdf>

###f)
(More on this in Module 6.)

There are many ways to perform model selection for multiple linear regression. One possibility is best subsets, which can be done using the `regsubsets` function from library `leaps`. 
Assume your "full" model fitted by `lm` is called  `fit`. You may define `x` from `model.matrix(fit)[,-1]` (not including the intercept term), and then run `best=regsubsets(x=model.matrix(fit)[,-1],y=rent99$rent)` and look at `summary(best)`.
Explain the print-out (with all the stars). Using the Mallows Cp (named `cp` in the list from `summary(best)`) will give the same result at using AIC (which is not available in this function). What is your preferred model?
Hint: look at the R-code in Problem 2 (Figure 3) from the TMA4267V2017 exam: [pdf](https://www.math.ntnu.no/emner/TMA4267/2017v/Exam/eV2017Enew.pdf), and maybe the solutions for the interprtation [pdf](https://www.math.ntnu.no/emner/TMA4267/2017v/Exam/mergedLFV2017.pdf)

----

## Problem 3: Simulations in R

1. Make R code that shows the interpretation of a 95% CI for $\beta_j$. Hint: Theoretical question 1.

2. Make R code that shows the interpretation of a 95% PI for a new response at ${\bf x}_0$. Hint: Theoretical question 2.

c. For simple linear regression, simulate at data set with homoscedastic errors and with heteroscedastic errors. Here is a suggestion of one solution - not using `ggplot`. You use `ggplot`. Why this? To see how things looks when the model is correct and wrong.

```{r, eval=FALSE}
#Homoscedastic errore
n=1000
x=seq(-3,3,length=n)
beta0=-1
beta1=2
xbeta=beta0+beta1*x
sigma=1
e1=rnorm(n,mean=0,sd=sigma)
y1=xbeta+e1
ehat1=residuals(lm(y1~x))
plot(x,y1,pch=20)
abline(beta0,beta1,col=1)
plot(x,e1,pch=20)
abline(h=0,col=2)
#Heteroscedastic errors
sigma=(0.1+0.3*(x+3))^2
e2=rnorm(n,0,sd=sigma)
y2=xbeta+e2
ehat2=residuals(lm(y2~x))
plot(x,y2,pch=20)
abline(beta0,beta1,col=2)
plot(x,e2,pch=20)
abline(h=0,col=2)
```

4. All this fuss about raw, standardized and studentized residuals- does really matter in practice? Below is one example where the raw residuals are rather different from the standardized, but the standardized is identical to the studentized. Can you come up with a simulation model where the standardized and studentized are very different? Hint: what about at smaller sample size?

```{r, eval=FALSE}
n=1000
beta=matrix(c(0,1,1/2,1/3),ncol=1)
set.seed(123)
x1=rnorm(n,0,1); x2=rnorm(n,0,2); x3=rnorm(n,0,3)
X=cbind(rep(1,n),x1,x2,x3)
y=X%*%beta+rnorm(n,0,2)
fit=lm(y~x1+x2+x3)
yhat=predict(fit)
summary(fit)
ehat=residuals(fit); estand=rstandard(fit); estud=rstudent(fit)
plot(yhat,ehat,pch=20)
points(yhat,estand,pch=20,col=2)
```

----

# <a id="further"> Further reading </a>

* Need details on the simple linear regression: From TMA4240/TMA4245 Statistics we have the thematic page for [Simple linear regression (in Norwegian)](https://wiki.math.ntnu.no/tma4245/tema/begreper/regression).
* Need more advanced thory: Theoretical version (no simple linear regression) from TMA4315 Generalized linear models H2017: [TMA4315M2: Multiple linear regression](https://www.math.ntnu.no/emner/TMA4315/2017h/2MLR.html)
* Slightly different presentation (more focus on multivariate normal theory): [Slides and written material from TMA4267](https://www.math.ntnu.no/emner/TMA4267/2017v/TMA4267V2017Part2.pdf)
* And, same source, but now [Slides and written material from TMA4267 Linear Statistical Models in 2017, Part 3: Hypothesis testing and ANOVA](http://www.math.ntnu.no/emner/TMA4267/2017v/TMA4267V2017Part3.pdf)
* [Videoes on YouTube by the authors of ISL, Chapter 2](https://www.youtube.com/playlist?list=PL5-da3qGB5IBSSCPANhTgrw82ws7w_or9)

---

# <a id="Rpackages"> R packages used</a>

\footnotesize
```{r, eval=FALSE}
# packages to install before knitting this R Markdown file
# to knit the Rmd
install.packages("knitr")
install.packages("rmarkdown")

# nice tables in Rmd
install.packages("kableExtra")

# cool layout for the Rmd
install.packages("prettydoc") # alternative to github

#plotting
install.packages("ggplot2") # cool plotting
install.packages("ggpubr") # for many ggplots
install.packages("GGally") # for ggpairs
#datasets
install.packages("ElemStatLearn") # for ozone data set
install.packages("gamlss.data")#rent index data set here
#methods
install.packages("nortest")#test for normality - e.g. Anderson-Darling

install.packages("car") # vif
library(Matrix)
install.packages("reshape")
install.packages("corrplot")
install.packages("tidyverse")
```

