Enter the General Linear Model

class: center, middle

# The General Linear Model - Mixing in Categories, Nonlinearities, and All That

![](Images/mmi/acl-bg.jpg)

---
# The General Linear Model

`$$\Large \boldsymbol{Y} \sim \mathcal{N}(\boldsymbol{Y}, \sigma^2\boldsymbol{I})$$`
$$\Large \boldsymbol{\hat{Y}} = \boldsymbol{\beta X}  $$

-   This is the Matrix formulation of `$\sum\beta_jx_j$`  
  
-   This equation is huge. X can be anything - categorical,
    continuous, squared, sine, etc.

-   There can be straight additivity, or interactions

-   So far, the only model we've used with >1 predictor is ANOVA

---

# So Many Ways to Build a Model
.large[

1. Mixing Categorical and Continuous Variables
  
2. Nonlinearities

3. Interaction Effects
  
]

---

## Combining Categorical Variables and Continuous Variables (Analysis of Covariance)

`$$\Large \boldsymbol{Y} = \boldsymbol{\beta X}$$`

translates to

`$$\widehat{y_i} = \beta_0 + \beta_{1}x_1  + \sum_{j=2}^{k}\beta_{ij}x_{ij}$$`  
`$$x_{ij} = 0,1$$`

- Here, we have many levels of a group in `$x_{ij}$`

-   Often used to correct for a gradient or some continuous variable affecting outcome

-   OR used to correct a regression due to additional groups that may throw off slope estimates  
      - e.g. Simpson's Paradox: A positive relationship between test scores and academic performance can be masked by gender differences

---
# Neanderthals and Mixing Categorical and Continuous Variables

![:scale 60%](./Images/mlr/neanlooking.jpeg)

Who had a bigger brain: Neanderthals or us?

---
# The Means Look the Same...

---

# But there appears to be a Relationship Between Body and Brain Mass

---

# And Mean Body Mass is Different

---
class: center, middle

![](Images/mlr/redo_analysis.jpg)

---
# Mixing Categories and Continuous Variables: Analysis of Covariance

`$$lnbrain_i \sim \mathcal{N}(\widehat{lnbrain_i}, \sigma_2)$$`
`$$\widehat{lnbrain_i} = \beta_0 + \beta_{1}lnmass_1  + \sum_{j=2}^{k}\beta_{ij}species_{ij}$$`  
`$$species_{ij} = 0\; or\; 1\; if\; neandertal$$`

Evaluate a categorical effect(s), controlling for a *covariate* (parallel lines)

Groups modify the *intercept*.

---
# Fit with Our Engine: Ordinary Least Squares

```r
neand_lm <- lm(lnbrain ~ species + lnmass, 
               data=neand)

neand_dat <- augment(neand_lm)
```

---
# Does Your Model Match Your Data?

---

# Are There Weird Patterns in Your Residuals?

---

# Did You Modelboss Too Close to the Correlated Suns?

---

# Did Our Predictors Matter?

<table class="table table-striped" style="margin-left: auto; margin-right: auto;">
 <thead>
  <tr>
   <th style="text-align:left;"> term </th>
   <th style="text-align:right;"> sumsq </th>
   <th style="text-align:right;"> df </th>
   <th style="text-align:right;"> statistic </th>
   <th style="text-align:right;"> p.value </th>
  </tr>
 </thead>
<tbody>
  <tr>
   <td style="text-align:left;"> species </td>
   <td style="text-align:right;"> 0.028 </td>
   <td style="text-align:right;"> 1 </td>
   <td style="text-align:right;"> 6.204 </td>
   <td style="text-align:right;"> 0.017 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> lnmass </td>
   <td style="text-align:right;"> 0.130 </td>
   <td style="text-align:right;"> 1 </td>
   <td style="text-align:right;"> 29.276 </td>
   <td style="text-align:right;"> 0.000 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Residuals </td>
   <td style="text-align:right;"> 0.160 </td>
   <td style="text-align:right;"> 36 </td>
   <td style="text-align:right;"> NA </td>
   <td style="text-align:right;"> NA </td>
  </tr>
</tbody>
</table>

---

# What are the Association Between Predictors and Response?
<table class="table table-striped" style="margin-left: auto; margin-right: auto;">
 <thead>
  <tr>
   <th style="text-align:left;"> term </th>
   <th style="text-align:right;"> sumsq </th>
   <th style="text-align:right;"> df </th>
   <th style="text-align:right;"> statistic </th>
   <th style="text-align:right;"> p.value </th>
  </tr>
 </thead>
<tbody>
  <tr>
   <td style="text-align:left;"> species </td>
   <td style="text-align:right;"> 0.028 </td>
   <td style="text-align:right;"> 1 </td>
   <td style="text-align:right;"> 6.204 </td>
   <td style="text-align:right;"> 0.017 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> lnmass </td>
   <td style="text-align:right;"> 0.130 </td>
   <td style="text-align:right;"> 1 </td>
   <td style="text-align:right;"> 29.276 </td>
   <td style="text-align:right;"> 0.000 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Residuals </td>
   <td style="text-align:right;"> 0.160 </td>
   <td style="text-align:right;"> 36 </td>
   <td style="text-align:right;"> NA </td>
   <td style="text-align:right;"> NA </td>
  </tr>
</tbody>
</table>
---

# What Does it Look Like to See the Unique Effect of One Driver?

<table class="table table-striped" style="margin-left: auto; margin-right: auto;">
 <thead>
  <tr>
   <th style="text-align:left;"> term </th>
   <th style="text-align:right;"> estimate </th>
   <th style="text-align:right;"> std.error </th>
   <th style="text-align:right;"> statistic </th>
   <th style="text-align:right;"> p.value </th>
  </tr>
 </thead>
<tbody>
  <tr>
   <td style="text-align:left;"> (Intercept) </td>
   <td style="text-align:right;"> 5.188 </td>
   <td style="text-align:right;"> 0.395 </td>
   <td style="text-align:right;"> 13.126 </td>
   <td style="text-align:right;"> 0.000 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> speciesrecent </td>
   <td style="text-align:right;"> 0.070 </td>
   <td style="text-align:right;"> 0.028 </td>
   <td style="text-align:right;"> 2.491 </td>
   <td style="text-align:right;"> 0.017 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> lnmass </td>
   <td style="text-align:right;"> 0.496 </td>
   <td style="text-align:right;"> 0.092 </td>
   <td style="text-align:right;"> 5.411 </td>
   <td style="text-align:right;"> 0.000 </td>
  </tr>
</tbody>
</table>
--

- Intercept is the brain mass of a neanderthal with 0 body mass

- speciesrecent is the change in brain mass for a recent human, still assuming 0 body mass

- lnmass is the association of log body mass with log brain mass

---

# What Does it Mean Visually?

---

# We can run comparisons at equal body mass

---

# So Many Ways to Build a Model
.large[

1. Mixing Categorical and Continuous Variables
  
2. .red[ Nonlinearities]

3. Interaction Effects
  
]

---
background-image: url("Images/mlr/fires.jpg")
background-size: cover

<div style="bottom:100%; text-align:left; background:goldenrod"><h4>Five year study of wildfires & recovery in Southern California shrublands in 1993. 90 plots (20 x 50m)</h4>
(data from Jon Keeley et al.)</div>

---

class:center, middle

# What is the effect of cover on richness?

---

# The Linear Effect of Cover on Richness

---

# But... What's Up With That Bend?

---

# How Do We Add Nonlinearities to our Linear Model?

`$$y_{i} \sim N(\widehat{y_{i}}, \sigma^{2} )$$`
`$$\widehat{y_{i}} = \beta_{0} + \sum \beta_{j}x_{ij}$$`

- What if `$x_1$` is a linear term and `$x_2 = x_1^2$`?

- Adding nonlinear terms is just adding more predictors to a multiple linear regression model

- Sometimes to reduce collinearity/make our model make more sense, we need to center x - i.e. use `$x-\bar{x}$` - to reduce SE in parameters

---

# A Nonlinear Effect of Cover?

```r
mod_sq <- lm(rich ~ cover + I(cover^2), data = keeley)
```

---

# Which is Better? Maybe A Cubic? Cross-Validate with LOO!

<table class="table table-striped" style="margin-left: auto; margin-right: auto;">
 <thead>
  <tr>
   <th style="text-align:right;"> order </th>
   <th style="text-align:right;"> mse </th>
  </tr>
 </thead>
<tbody>
  <tr>
   <td style="text-align:right;"> 1 </td>
   <td style="text-align:right;"> 210.27 </td>
  </tr>
  <tr>
   <td style="text-align:right;"> 2 </td>
   <td style="text-align:right;"> 200.16 </td>
  </tr>
  <tr>
   <td style="text-align:right;"> 3 </td>
   <td style="text-align:right;"> 208.85 </td>
  </tr>
  <tr>
   <td style="text-align:right;"> 4 </td>
   <td style="text-align:right;"> 196.36 </td>
  </tr>
  <tr>
   <td style="text-align:right;"> 5 </td>
   <td style="text-align:right;"> 211.24 </td>
  </tr>
  <tr>
   <td style="text-align:right;"> 6 </td>
   <td style="text-align:right;"> 1323.11 </td>
  </tr>
</tbody>
</table>

---
# 2 or 4? Or Other Nonlinear Terms?

---

class: center, middle

# Cool story. Other than nonlinear transformations of y or x, are there other forms of nonlinearity we should watch for?

---

# So Many Ways to Build a Model
.large[

1. Mixing Categorical and Continuous Variables
  
2. Nonlinearities

3. .red[Interaction Effects]
  
]
---

# What if Fire Severity Only Affects Species Richness in Old Stands?

---
# Interaction Effects

- The effect of one predictor cannot be known without knowing the level of another

- Common in nature!

- Some are inevitable (e.g., the effects of disturbance depend on something being there to disturb)

- Some are.... quite tricky, but the spice of biological life!

- Exercise: Think of an interaction effect you have encountered in biology or elsewhere!

---

# Interaction Effects in Models

`$$y_{i} \sim N(\widehat{y_{i}}, \sigma^{2} )$$`
`$$\widehat{y_{i}} = \beta_{0} +  \beta_{1}x_{1i} +  \beta_{2}x_{2i} +   \beta_{3}x_{1i}x_{2i}$$` 
`$$x_{2i} = 0,1$$`
--

If we consider `$x_1 * x_2$` a variable in it's own right - i.e. `$x_3$` then..... this is just.......

`$$\widehat{y_{i}} = \beta_{0} +  \sum\beta_{j}x_{ij}$$`

IT'S ALL THE SAME THING!

---

# Interaction Effects in Models - and F-Tests!

```r
mod_int <- lm(rich ~ firesev * age_break, data = keeley)
```

<table class="table table-striped" style="margin-left: auto; margin-right: auto;">
 <thead>
  <tr>
   <th style="text-align:left;"> term </th>
   <th style="text-align:right;"> sumsq </th>
   <th style="text-align:right;"> df </th>
   <th style="text-align:right;"> statistic </th>
   <th style="text-align:right;"> p.value </th>
  </tr>
 </thead>
<tbody>
  <tr>
   <td style="text-align:left;"> firesev </td>
   <td style="text-align:right;"> 1918.314 </td>
   <td style="text-align:right;"> 1 </td>
   <td style="text-align:right;"> 10.122 </td>
   <td style="text-align:right;"> 0.002 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> age_break </td>
   <td style="text-align:right;"> 202.176 </td>
   <td style="text-align:right;"> 1 </td>
   <td style="text-align:right;"> 1.067 </td>
   <td style="text-align:right;"> 0.305 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> firesev:age_break </td>
   <td style="text-align:right;"> 1060.005 </td>
   <td style="text-align:right;"> 1 </td>
   <td style="text-align:right;"> 5.593 </td>
   <td style="text-align:right;"> 0.020 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Residuals </td>
   <td style="text-align:right;"> 16298.526 </td>
   <td style="text-align:right;"> 86 </td>
   <td style="text-align:right;"> NA </td>
   <td style="text-align:right;"> NA </td>
  </tr>
</tbody>
</table>

---

# Interaction Effects What does it mean?

<table class="table table-striped" style="margin-left: auto; margin-right: auto;">
 <thead>
  <tr>
   <th style="text-align:left;"> term </th>
   <th style="text-align:right;"> estimate </th>
   <th style="text-align:right;"> std.error </th>
   <th style="text-align:right;"> statistic </th>
   <th style="text-align:right;"> p.value </th>
  </tr>
 </thead>
<tbody>
  <tr>
   <td style="text-align:left;"> (Intercept) </td>
   <td style="text-align:right;"> 53.464 </td>
   <td style="text-align:right;"> 6.581 </td>
   <td style="text-align:right;"> 8.123 </td>
   <td style="text-align:right;"> 0.000 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> firesev </td>
   <td style="text-align:right;"> 0.117 </td>
   <td style="text-align:right;"> 1.626 </td>
   <td style="text-align:right;"> 0.072 </td>
   <td style="text-align:right;"> 0.943 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> age_breakold stand </td>
   <td style="text-align:right;"> 16.348 </td>
   <td style="text-align:right;"> 8.961 </td>
   <td style="text-align:right;"> 1.824 </td>
   <td style="text-align:right;"> 0.072 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> firesev:age_breakold stand </td>
   <td style="text-align:right;"> -4.730 </td>
   <td style="text-align:right;"> 2.000 </td>
   <td style="text-align:right;"> -2.365 </td>
   <td style="text-align:right;"> 0.020 </td>
  </tr>
</tbody>
</table>

- The intercept is the number of species when fire severity is 0 and stand age is young

- The firesev effect is the effect of fire for young stands

- The age effect is the increase in species in an old stand - but only if fire severity is 0

- The interaction is the change in the fire severity slope when the stand is old

---
# This is Why I Hate Staring at Coefficients

---
class: center, middle

# Maybe if Age Was Continuous....

---

# Models with Continuous Interactions - The Same Model!
`$$y_{i} \sim \mathcal{N}(\widehat{y_{i}}, \sigma^2)$$`
`$$\widehat{y_{i}} = \beta_{0} +  \beta_{1}x_{1i} +  \beta_{2}x_{2i} +  \beta_{3}x_{1i}x_{2i}$$`

---

# Fit and Go!

.large[

```r
mod_int_c <- lm(rich ~ firesev * age, data = keeley)
```
]

<table class="table table-striped" style="margin-left: auto; margin-right: auto;">
 <thead>
  <tr>
   <th style="text-align:left;"> term </th>
   <th style="text-align:right;"> sumsq </th>
   <th style="text-align:right;"> df </th>
   <th style="text-align:right;"> statistic </th>
   <th style="text-align:right;"> p.value </th>
  </tr>
 </thead>
<tbody>
  <tr>
   <td style="text-align:left;"> firesev </td>
   <td style="text-align:right;"> 1361.509 </td>
   <td style="text-align:right;"> 1 </td>
   <td style="text-align:right;"> 7.877 </td>
   <td style="text-align:right;"> 0.006 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> age </td>
   <td style="text-align:right;"> 466.926 </td>
   <td style="text-align:right;"> 1 </td>
   <td style="text-align:right;"> 2.701 </td>
   <td style="text-align:right;"> 0.104 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> firesev:age </td>
   <td style="text-align:right;"> 2228.687 </td>
   <td style="text-align:right;"> 1 </td>
   <td style="text-align:right;"> 12.894 </td>
   <td style="text-align:right;"> 0.001 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Residuals </td>
   <td style="text-align:right;"> 14865.094 </td>
   <td style="text-align:right;"> 86 </td>
   <td style="text-align:right;"> NA </td>
   <td style="text-align:right;"> NA </td>
  </tr>
</tbody>
</table>

---

# Did You Make Sure your Train Wasn't On Fire?

---

# But Was This a Good Idea? AIC Edition

<table class="table table-striped" style="margin-left: auto; margin-right: auto;">
 <thead>
  <tr>
   <th style="text-align:left;"> Modnames </th>
   <th style="text-align:right;"> K </th>
   <th style="text-align:right;"> AIC </th>
   <th style="text-align:right;"> Delta_AIC </th>
   <th style="text-align:right;"> AICWt </th>
  </tr>
 </thead>
<tbody>
  <tr>
   <td style="text-align:left;"> fire*age </td>
   <td style="text-align:right;"> 5 </td>
   <td style="text-align:right;"> 725.035 </td>
   <td style="text-align:right;"> 0.000 </td>
   <td style="text-align:right;"> 0.991 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> fire + age </td>
   <td style="text-align:right;"> 4 </td>
   <td style="text-align:right;"> 735.608 </td>
   <td style="text-align:right;"> 10.573 </td>
   <td style="text-align:right;"> 0.005 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> age </td>
   <td style="text-align:right;"> 3 </td>
   <td style="text-align:right;"> 736.034 </td>
   <td style="text-align:right;"> 10.998 </td>
   <td style="text-align:right;"> 0.004 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> fire </td>
   <td style="text-align:right;"> 3 </td>
   <td style="text-align:right;"> 740.506 </td>
   <td style="text-align:right;"> 15.470 </td>
   <td style="text-align:right;"> 0.000 </td>
  </tr>
</tbody>
</table>

---

# Coefficients are Clearer

<table class="table table-striped" style="margin-left: auto; margin-right: auto;">
 <thead>
  <tr>
   <th style="text-align:left;"> term </th>
   <th style="text-align:right;"> estimate </th>
   <th style="text-align:right;"> std.error </th>
   <th style="text-align:right;"> statistic </th>
   <th style="text-align:right;"> p.value </th>
  </tr>
 </thead>
<tbody>
  <tr>
   <td style="text-align:left;"> (Intercept) </td>
   <td style="text-align:right;"> 42.935 </td>
   <td style="text-align:right;"> 7.851 </td>
   <td style="text-align:right;"> 5.469 </td>
   <td style="text-align:right;"> 0.000 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> firesev </td>
   <td style="text-align:right;"> 3.106 </td>
   <td style="text-align:right;"> 1.863 </td>
   <td style="text-align:right;"> 1.667 </td>
   <td style="text-align:right;"> 0.099 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> age </td>
   <td style="text-align:right;"> 0.827 </td>
   <td style="text-align:right;"> 0.313 </td>
   <td style="text-align:right;"> 2.641 </td>
   <td style="text-align:right;"> 0.010 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> firesev:age </td>
   <td style="text-align:right;"> -0.230 </td>
   <td style="text-align:right;"> 0.064 </td>
   <td style="text-align:right;"> -3.591 </td>
   <td style="text-align:right;"> 0.001 </td>
  </tr>
</tbody>
</table>

- Intercept is when all predictors are 0  
  <br>
- Additive slopes are when the OTHER predictor is 0  
<br>
- Interaction is how they modify each other

---

# What if a Predictor Value of 0 Makes No Sense: Centering

.large[

`$$X_i - \bar{X}$$`

]

- Additive coefficients are now the effect of a predictor at the mean value of the other predictors

- Intercepts are at the mean value of all predictors

- Visualization will keep you from getting confused!

---

# Model Visualization: The Effect of One Predictor at Different Levels of the Other

---

# Final Thoughts

- Whew! That's a lot of ways to build a model!

- Note, we are still thinking in terms of normal error - and an additive model

- There are plenty of extensions 
     - nonlinear non-additive models  
     - Generalized Linear Models with Non-Normal Error  
     
--

- But, additive Gaussian models are a safe starting point

- AND extensions.... really all just have linear models at their core

- It's the BIOLOGY and your intuition of what model to build that should **always** be the core of your modeling process