01_Regression-6.Rmd

---
title: "Chapter 1: Linear regression"
subtitle: "Multicollinearity"
author: "Joris Vankerschaver"
header-includes:
  - \useinnertheme[shadow=true]{rounded}
  - \usecolortheme{rose}
  - \setbeamertemplate{footline}[frame number]
  - \usepackage{color}
  - \usepackage{graphicx}
  - \usepackage{amsmath}
  - \graphicspath{{./images/01-linear-regression}}
output: 
  beamer_presentation:
    theme: "default"
    keep_tex: true
    includes:
      in_header: columns.tex
---

```{r include = FALSE}
needles <- read.table("./datasets/01-linear-regression/needles.txt", header = T, sep = "\t", dec = ".")
attach(needles)
model_l5 <- lm(length ~ nitrogen * phosphor + potassium 
               + phosphor * residu)
model_l8 <- lm(length ~ nitrogen *  phosphor + potassium)

body.fat <- read.table("./datasets/01-linear-regression/bodyfatNKNW.txt", header = T, dec = ".")
attach(body.fat)
```

## Multicollinearity 

- There is **multicollinearity** when 2 or more predictors are correlated
- \alert{Can possibly cause problems}: if there is strong correlation between 2 predictors $X_1$ and $X_2$, it becomes difficult to discern effect of $X_1$ of effect of $X_2$
\begin{exampleblock}{Example} 
If $X_1=X_2$, then 
\[E(Y|X_1,X_2)=\beta_0+\beta_1X_1+\beta_2X_2=\beta_0+(\beta_1+\beta_2)X_1\]
\end{exampleblock}

\begin{exampleblock}{Consequences}
\begin{itemize}
\item Numerically instable estimates
\item Estimates with large standard errors
\item Difficult interpretation of coefficients
\end{itemize}
\end{exampleblock}

## Diagnosing multicollinearity

Multicollinearity can be recognized through:


- **Instability**:
  - Large changes in coefficients after adding a predictor
  - Very wide confidence intervals
  - Unexpected results 
- **Strong correlation** between predictors:
  - Example: usually strong correlation between $X_f$ and $X_fX_s$
  - Can sometimes be eliminated by **centering** (subtracting the mean): 
\[
  X \to X-\bar X.
\]


## Impact of centering

\centering
```{r echo=FALSE}
x <- rnorm(100, 5, 1)
cx <- x - mean(x)
y <- x^2
cy <- cx^2
r1 <- cor(x,y)
r2 <- cor(cx,cy)
par(mfrow=c(1,2))
plot(x, y, pch = 20, main = paste("Correlation = ", round(r1,2)), xlab = "x", ylab = expression(x^2))
abline(lm(y ~ x), col = "blue")
plot(cx, cy, pch = 20, main = paste("Correlation = ", round(r2,2)), xlab = "cx", ylab = expression(cx^2))
abline(lm(cy ~ cx), col = "blue")
``` 

## Scatterplot matrix - before centering

\centering
```{r echo=FALSE}
nitrogenphospor <- nitrogen * phosphor 
minerals <- cbind(nitrogen, phosphor, potassium, nitrogenphospor)
colnames(minerals)[4] <- "nitrogen*phosphor"
pairs(minerals)
```

## Scatterplot matrix - after centering

\centering
```{r echo=FALSE}
cnitrogen <- nitrogen - mean(nitrogen)
cphosphor <- phosphor - mean(phosphor)
cpotassium <- potassium - mean(potassium)
cresidu <- residu - mean(residu)
cnitrogenphospor <- cnitrogen * cphosphor 
minerals2 <- cbind(cnitrogen, cphosphor, cpotassium, cnitrogenphospor)
colnames(minerals2)[4] <- "cnitrogen*cphosphor"
pairs(minerals2)
```

## Diagnosing multicollinearity

Previous diagnostics are \alert{limited}:
\begin{exampleblock}{Example} 
\begin{itemize}
\item Even if pairwise correlations between predictors $X_1,X_2,X_3$ low, there can be strong multicollinearity.
\item E.g., when strong correlation between $X_1$ and a linear combination of $X_2$ and $X_3$.
\end{itemize}
\end{exampleblock}
\begin{block}{Variance inflation factor for $k^{th}$ coefficient}
\[\textrm{VIF}_k=\left(1-R_k^2\right)^{-1}\]
with $R_k^2$ the $R^2$ of linear regression of $k^{th}$ predictor on other predictors.
\end{block}

## Interpretation VIF

- $\textrm{VIF}_k \geq 1$; $\textrm{VIF}_k=1$ if $k^{th}$ predictor \alert{not} linearly associated with other predictors.
- Expresses how much larger variance on $k^{th}$ coefficient is than when all predictors were independent.
- Average quadratic distance between estimated and true coefficients is proportionate with average VIF.
- Critical multicollinearity: maximum VIF of at least 10.


## Variance inflation factors

\centering
```{r echo=FALSE, message=FALSE}
library(car)
model_l8b <- lm(length ~ cnitrogen *  cphosphor + cpotassium)
model_l5b <- lm(length ~ cnitrogen * cphosphor + cpotassium + cphosphor * cresidu)
v1 <- vif(model_l5)
v2 <- vif(model_l5b)
par(mfrow=c(1,2))
plot(v1, pch = 20, main = "Before centering", ylab = "VIF", xlab = "Parameter number")
abline(h=10, lty = 2, col = "blue")
z <- seq(1.4,2.8,0.2)
plot(v2, pch = 20, main = "After centering", ylab = "VIF", xlab = "Parameter number", yaxp = c(1.4, 2.8, 7))
```

## Simpler interpretation of coefficients

Coefficients (without centering)
\small
```{r echo=FALSE}
summary(model_l8)$coefficients
```
\normalsize
Coefficients (with centering)
\small
```{r echo=FALSE}
summary(model_l8b)$coefficients
```

## Example: Prediction body fat

\begin{exampleblock}{}
\begin{itemize}
\item Determining percentage body fat difficult and expensive
\item Study investigates association between
\begin{itemize}
\item $Y$: body fat
\item $X_1$: triceps skinfold thickness
\item $X_2$: thigh circumference
\item $X_3$: midarm circumference
\end{itemize}
\item 20 healthy women between 25 and 34 years old
\end{itemize}
\end{exampleblock}

## Analysis in \texttt{R}

\footnotesize
```{r echo=FALSE}
model_bf <- lm(bodyfat ~ triceps.skinfold.thickness 
               + thigh.circumference + midarm.circumference)
summary(model_bf)
```

## Scatterplot matrix

\centering
```{r echo=FALSE}
pairs(body.fat[,-4])
```

## Variance inflation factors

\small
```{r echo=FALSE, warning=FALSE}
library(car)
vif_bodyfat <- vif(model_bf)
as.data.frame(vif_bodyfat)
```
\normalsize

- VIF on average 460.
- Large VIF for midarm circumference, although weakly correlated with other predictors.
- **How to correct for multicollinearity?**
  - Centering variables only valid option when higher order terms are in play.
  - Combine predictors, e.g., through principal component regression.
  - Ridge regression: allow some bias in exchange for increased precision and lower risk of overfitting.

## Multicollinearity and confounding

- A lot of textbooks advise to remove predictors from model in case of multicollinearity
- However, multicollinearity can also indicate strong confounding!