---
title: "STAT 224 Chapter 9"
output:
  html_document:
    theme: cerulean
    css: styles.css
    code_folding: show
    df_print: kable
    number_sections: yes
    self_contained: yes
    toc: FALSE
    toc_depth: 4
    toc_float: FALSE
params:
  dotabs:  TRUE
---


```{r test-main, child = 'preamble.Rmd'}
```

# Chapter 9: Working with Collinear Data `r if(params$dotabs) paste('{.tabset}')`

### Example

#### Another Problem to Solve: (Multi)Collinearity

Collinearity (also called multicollinearity) refers to the situation where   
some or all predictors in the model are substantially correlated with each other

Linear relationships among predictors might not be only pairwise.   
Perhaps one of the predictor is linearly related to some subset of other predictors.

#### Example: Student achievement and school facilities

##### The scenario: Civil Rights Act of 1964

Separate (and allegedly equal) education still leads to differences in educational achievement.

Question: How is student achievement related to school resources?

A sample of 70 schools and their political/social context resulted in data on many variables.     
Data were combined into indices/measures of achievement and school characteristics

  + ACHV: An index of student achievement (the response variable)
  + SCHOOL: A measure of school resources
  + FAM: A measure of faculty credentials in each school
  + PEER: A measure of the influence of peer groups in each school
  
For achievement, high values are good.

The predictors are expected to be positively correlated with achievement.

##### The proposed model and estimate: A cautionary tale

$ACHV = \beta_0 + \beta_1 FAM + \beta_2 PEER + \beta_3 SCHOOL + \epsilon$

```{r}
achvData <- read.delim("http://statistics.uchicago.edu/~collins/data/RABE5/P236.txt")
glimpse(achvData)
head(achvData)
```


```{r}
lmfit.full <- lm(ACHV ~ FAM + PEER + SCHOOL, data = achvData)
tidy(lmfit.full)
glance(summary(lmfit.full))
```

**SEE ANYTHING ODD IN THE P-VALUES?**

<br> <br> <br> <br> <br> <br>

The overall F-test has a p-value of 0.0015.

$H_0: \beta_1 = \beta_2 = \beta_3 = 0$     
vs.      
$H_1:$ at least one $\beta_j \ne 0$

We reject $H_0$ and conclude that at least on predictor is significant.

But, **none** of the individual t-tests shows a significant predictor!!!


**SEE ANYTHING ODD ABOUT THE SIGNS OF THE REGRESSION COEFFICIENTS?**

```{r}
tidy(lmfit.full)
```

<br> <br> <br> <br> <br> <br>

The coefficient for SCHOOL is negative!

Hopefully, having more school resources is associated with higher achievement!

#### What's the problem then?

Usually, these types of oddities are the result of multicollinearity:    
Linear relationships among the predictor variables.

Are the predictors correlated?

```{r, fig.height=4}
ggpairs(achvData)
```

The predictors are strongly correlated with eachother.

### Issues

#### What's the big deal about multicollinearity?

  1. Since the predictors provide "overlapping" information, **results may be ambiguous**
     + A portion of the information from each predictor is redundant, but not all
     + Predictor variables might each serve as a proxy for the others in the    
       regression equation without affecting the total explanatory power.
     + Do we need all the predictors in the model?
     + If not, which ones should we drop? 
     + How does this choice effect inference, prediction, or forecasting?    
      (for prediction/forecasting, we are less concerned)
    
  2. **We lose the simple interpretation of a regression coefficient**   
    "while all other predictors are held constant"
  
  3. If we proceed in the presence of moderate multicollinearity
     + **Variance estimates are typically overstated**
       + So, standard errors typically overstated
       + Confidence intervals wider than needed
       + t-statistics smaller than they should be
     + These are missed opportunities to notice important predictors
    
  4. If we proceed in the presence of strong multicollinearity
     + **The estimates of $\beta$s can be very unstable**
     + $\widehat\beta$s may change substantially when other predictors are added/removed   
       or observations are added/deleted (even if not outliers)
    
##### Multicollinearity is not a modeling error

In fact, it cannot be detected by examination of residuals.

  + It is a condition of deficient data.
  + It arises from a lack of independent information (about the response) coming from predictors

Most any observed variables (data) have a little correlation with each other.

  + How much multicollinearity can we permit and still have a useful, valid model?
  + What do we do if multicollinearity is so serious that we cannot ignore it?
  
In principle, multicollinearity can be completely avoided in carefully designed experiments

  + The factors we expect to effect the response can be collected at values   
  that make the predictor variables statistically independent ("orthogonal").
  
  + Check out STAT 22200 Linear Models and Experimental Design (offered each Spring)
  
In practice, we select several predictors believed to be related to the outcome,   
and thus... likely to be associated with each other as well.
  
### Example
  
#### Example Again: Student achievement and school facilities

##### The scenario: Civil Rights Act of 1964

Question: How is student achievement related to school resources?

A sample of 70 schools and their political/social context resulted in data on many variables.     
Data were combined into indices/measures of achievement and school characteristics

  + ACHV: An index of student achievement (the response variable)
  + SCHOOL: A measure of school resources
  + FAM: A measure of faculty credentials in each school
  + PEER: A measure of the influence of peer groups in each school
  
For achievement, high values are good.

The predictors are expected to be positively correlated with achievement.

##### The proposed model and estimate

$ACHV = \beta_0 + \beta_1 FAM + \beta_2 PEER + \beta_3 SCHOOL + \epsilon$

```{r}
tidy(lmfit.full)
```

**REMINDERS:**

  + If the predictors provide "overlapping" information, results may be ambiguous   
    A portion of the information from each predictor may be redundant, but not all
  + Predictor variables might each serve as a proxy for the others in the    
    regression equation without affecting the total explanatory power.

If we proceed in the presence of strong multicollinearity   

  + The estimates of $\beta$'s can be very unstable
  + They may change substantially when other predictors are added/removed   
  
##### $\beta$s change when predictors added/removed

Start with models with just one predictor.
```{r}
lmfit.only.SCHOOL <- lm(ACHV ~ SCHOOL, data = achvData)
lmfit.only.FAM <- lm(ACHV ~ FAM, data = achvData)
lmfit.only.PEER <- lm(ACHV ~ PEER, data = achvData)
# lmfit.not.SCHOOL <- lm(ACHV ~ FAM + PEER, data = achvData)
# lmfit.not.PEER <- lm(ACHV ~ FAM + SCHOOL, data = achvData)
```

Adding PEER to a model with only SCHOOL:
```{r}
lmfit.not.FAM <- lm(ACHV ~ PEER + SCHOOL, data = achvData)
tidy(lmfit.only.SCHOOL)
tidy(lmfit.not.FAM)
tidy(lmfit.only.PEER)
```

Adding PEER to the model changes the sign of SCHOOL coefficient.    
The size of the SCHOOL coefficient is about the same,   
but the SE is 5 times larger!   
At the same time, the PEER t-statistic is more than cut in half.

The overall fit of the models $\left(R^2\right)$ is essentially unchanged.

$ACHV = \beta_0 + \beta_1 FAM + \beta_2 PEER + \beta_3 SCHOOL + \epsilon$
```{r}
glance(summary(lmfit.full))
```

$ACHV = \beta_0 + \beta_3 SCHOOL + \epsilon$
```{r}
glance(summary(lmfit.only.SCHOOL))
```

$ACHV = \beta_0 + \beta_2 PEER + \beta_3 SCHOOL + \epsilon$
```{r}
glance(summary(lmfit.not.FAM))
```

$ACHV = \beta_0 + \beta_2 PEER + \epsilon$
```{r}
glance(summary(lmfit.only.PEER))
```


<br>

If we were to keep just one variable, which would it be?

```{r}
kable(cor(achvData))
```

<!--

$ACHV = \beta_0 + \beta_1 FAM + \epsilon$
```{r}
glance(summary(lmfit.only.FAM))
```

$ACHV = \beta_0 + \beta_2 PEER + \epsilon$
```{r}
glance(summary(lmfit.only.PEER))
```

$ACHV = \beta_0 + \beta_3 SCHOOL + \epsilon$
```{r}
glance(summary(lmfit.only.SCHOOL))
```

-->

Technically, the PEER predictor is best (highest correlation with response),    
but for the purposes of the study,    
I suspect researchers would prefer to keep either FAM or SCHOOL.

### Detection

#### Detecting multicollinearity

##### Detection: Oddities in model estimates/tests

Look for overall (model) F-statistic significant
($H_0: \beta_1 = \beta_2 = \cdots = \beta_p = 0$),    
but  **ALL** individual $t$-tests nonsignificant.

However, a significant F-statistic ($H_0: \beta_1 = \beta_2 = \cdots = \beta_p = 0$)
and **SOME** nonsignificant $t$-tests for individual $\beta$s    
does **NOT** necessarily imply multicollinearity. 

There simply could be non-important variables in the model.

Look for a $\beta$ coefficient for a variable have an opposite sign    
from what you (and everyone else) would expect.   
(This may also be a result of omitting confounders: November 29 lecture.)

##### Detection: pairwise linear relationships (correlation, scatterplots)

```{r fig.height=4}
ggpairs(achvData)
```

Indeed, each predictor is positively correlated with the response, as expected.

<br>

But, each predictor also strongly correlated with the other predictors.

We can't even answer questions like:   
If faculty credentials are high, but school resources are low,
what is effect of peer influence on achievement?

There is no such data.  Could we obtain it?    
Probably not, because it appears that in this sample of public schools,    
no such situations exist.

These variables, perhaps by their design and practicality, are inextricably linked.   
Schools with more money have more money for resources    
and more money to hire higer-credentialed teachers.

##### Detection more complex linear relationships (partial $R^2$)

Not every multicollinearity situation can be detected by pairwise relationships:    
correlation, scatterplots.

Other measures include the partial $R^2$ and its companion:     
the variance inflation factor (VIF)

###### Calculation and Interpretation:

$R_j^2$ is the multiple correlation coefficient that results from predicting
predictor $x_j$ from a linear model with the other predictor variables.

For example, suppose there are 3 predictors (as in the achievement data)    
If $j=2$, we let $x_j = x_2$ temporarily be a response variable and consider the model

$X_2 = \beta_0 + \beta_1 x_1 + \beta_2 x_3 + \epsilon$

Find the $R^2$ from this model.    
Call it $R_j^2 = R_2^2$: the partial regression coefficient for $x_2$.


```{r}
lmfit.PEER <- lm(PEER ~ FAM + SCHOOL, data = achvData)
tidy(lmfit.PEER)
glance(summary(lmfit.PEER))

R2.PEER <- summary(lmfit.PEER)$r.squared
b0 <- tidy(lmfit.PEER)$estimate[1]
b1 <- tidy(lmfit.PEER)$estimate[2]
b2 <- tidy(lmfit.PEER)$estimate[3]
```

The partial correlation coefficient is $R_2^2 = `r R2.PEER`$    
$`r 100*R2.PEER`$% of the variability in PEER is explained by regression on FAM and SCHOOL:

$PEER \approx `r b0` + `r b1` FAM + `r b2` SCHOOL$

So, the data doesn't contain much information about how PEER varies    
when FAM and SCHOOL are **not** changed.

This complicates interpretation of the regression coefficient for PEER in the proposed model

$ACHV = \beta_0 + \beta_1 FAM + \beta_2 PEER + \beta_3 SCHOOL + \epsilon$

**REMINDER:** 
If we proceed in the presence of multicollinearity...   
We lose the simple interpretation of a regression coefficient   
    "while all other predictors are held constant"
    
##### Detection: Variance Inflation Factor (VIF)

Let $\displaystyle{\VIF_j = \frac{1}{1 - R_j^2}}$.

We call $\VIF_j$ the **variance inflaction factor** due to $x_j$.

The variance inflation factor for PEER: 
$\;\displaystyle{\VIF_j = \frac{1}{1 - R_j^2} 
\;=\; \frac{1}{1 - `r R2.PEER`} \;=\; `r 1 / (1-R2.PEER)`}$

R will calculate the VIF for you.  I used the `vif` function in the `car` package.

`library(car)`     
`lmfit <- lm(ACHV ~ FAM + PEER + SCHOOL, data = achvData)`    
`vif(lmfit)`

```{r}
kable(vif(lmfit.full))
```

<br> 

When the partial correlation coefficient $R_j^2$ is large, VIF is large.

...and if $R_j^2 = 0$ (orthogonal predictors), VIF = 1.

**Rough guide:** VIF $> 10$ indicates moderate multicollinearity is present.

<br>
**REMINDER:**
If we proceed in the presence of even moderate multicollinearity...

     + **Variance estimates are typically overstated**
       + So, standard errors typically overstated
       + Confidence intervals wider than needed
       + t-statistics smaller than they should be
     + These are missed opportunities to notice important predictors
    
<br>
    
The VIF is attempting to put a measure on whether overstated errors is a problem.

More specifically, the VIFs can be combined to measure the expected squared distance
(squared error) of OLS estimators from their true values

$\displaystyle D^2 = \sigma^2 \sum_{j=1}^p \VIF_j$

That's nice, but we can't calculate this from the data    
since multicollinearity causes $\widehat{\sigma}^2$ to overstate the variance.   
It's a circular reasoning issue.  
We can't estimate how much the variance is inflated by using an inflated variance estimate.

We can measure what we would expect if the predictors were orthogonal (a null hypothesis)   
and see how far away we are from that hypothesis using the VIFs.

<br>

Suppose the predictor variables are all orthogonal to each other.   

Then, $\VIF_j = 1$ for $j = 1, 2, \ldots, p$    

and $D^2 = \sigma^2 \sum_{j=1}^p \VIF_j = \sigma^2 \sum_{j=1}^p 1 =  \sigma^2 p$

Then, the ratio
$\displaystyle \frac{D^2}{D^2 \textrm{ if predictors orthogonal}}
\;=\; \frac{\sigma^2 \sum_{j=1}^p \VIF_j}{\sigma^2 p}
\;=\; \frac{1}{p} \sum_{j=1}^p \VIF_j
\;=\; \overline{\VIF}$

So, the average of the VIFs measures the squared error of the OLS estimators from their true values    
relative to the size of that error if the predictors were orthogonal.

For example, when the predictors are orthogonal and the regression model valid,   
then software can calculate $\se(\widehat{\beta}_j)$ and we would consider it to be accurate.

If multicollinearity is present,     
we estimate the true $\se(\widehat{\beta}_j)$ is overestimated by a factor of $\sqrt{\textrm{average } \VIF}$.

```{r}
VIFs <- vif(lmfit.full)
p <- length(VIFs)
p
avgVIF = sum(VIFs) / p
```

For the achievement data, $\overline{\VIF} = `r avgVIF`$.

So, the squared error in the OLS estimates is $`r round(avgVIF,0)`$ times larger    
than it would be if the predictor variables were orthogonal.

...and the reported standard errors are $\sqrt{`r round(avgVIF,0)`} = `r round(sqrt(avgVIF),1)`$
times larger than they should be.

Interpretation is qualitative.     

We don't know the distribution of $\overline{\VIF}$ (even if the predictors are orthogonal).    
It's an open problem in statistics, but not one we are clamoring to solve.    
The analyst's judgement should suffice.

### Remedies

#### What to do about multicollinearity?

The bad news: We cannot "fix" or eliminate multicollinearity for a given model.

For example, no transformation will help.

Some ideas...

  1. Get more and better data
  2. Omit redundant variables
  3. Centering
  4. Try "constrained regression" (Section 10.4)   
    This is the only section in Chapter 10 covered in this course.
  5. Principal components (not covered in this course, see Sections 10.1-10.3)

##### Get more and better data

Get more and/or better data with predictors that lie in the areas about which the current data are not informative. 

This is not always feasible. 

Stratified sampling to assure representativeness may aid in maximizing information in the data.

###### For the achievement data

We can't even answer questions like:   
If faculty credentials are high, but school resources are low,
what is effect of peer influence on achievement?

There is no such data.  Could we obtain it?    
Probably not, because it appears that in this sample of public schools,    
no such situations exist.

##### Omit redundant variables

Highly collinear predictors do not lend additional information    
while using degrees of freedom and creating other problems. 

There may be extra-statistical reasons (such as cost, reliability, etc) to just omit some predictors. 

Consequences can be evaluated.

##### Centering (Chapter 3)

Some degree of collinearity can be introduced when quadratic or interaction terms are added to a model. 

These can be often reduced by centering.    
Center predictors before constructing powers (square, etc) and interaction terms involving them. 

##### Try "constrained regression" (Section 10.4)

If predictors are functionally related, we can construct new variables that are   
functions of the predictors and include the new variables in the model instead.

This is equivalent to defining constraints on the coefficients.  

Omitting $x_j$ completely is a special case equivalent to fitting the model with the constraint that $\beta_j=0$.  

##### Principle components (not in this course)

Principal Components is an analysis technique performed on the $X$ matrix to discern possible structural effects. 

The method involves finding functions of $x$'s that are orthogonal and "best" capture the variability in $Y$.    
(Somewhat outside of scope of this course, see Chapter 10).


#### Prediction vs. explanation when multicollinearity

Explanation/interpretation of coefficients is pretty much impossible in the presence of strong multicollinearity, but predictions (within the scope of the data) are probably just fine.

Of course, the standard errors of those predictions will be inflated.

#### Other examples

  + Textbook example on French import and domestic consumption.
  + Other examples in Prof. Chen's lecture notes
  