Multiple Linear Regression

# Multiple Linear Regression
## Model Assessment
### Dr. Maria Tackett
### 09.30.19

---

### [Click for PDF of slides](10-model-assessment.pdf)

---

### Announcements

- Lab 05 **due Tuesday at 11:59p**

- HW 03 **due Wednesday at 11:59p**
    
- [Reading 06](https://www2.stat.duke.edu/courses/Fall19/sta210.001/reading/reading-06.html) for Wednesday 
    
---

## R packages

```r
library(tidyverse)
library(knitr)
library(broom)
library(Sleuth3) # ex0824 data
library(cowplot) # use plot_grid function
```

---

### Log Transformations

---

## Respiratory Rate vs. Age

- A high respiratory rate can potentially indicate a respiratory infection in children. In order to determine what indicates a "high" rate, we first want to understand the relationship between a child's age and their respiratory rate.

- The data contain the respiratory rate for 618 children ages 15 days to 3 years.

- **Variables**: 
    - <font class="vocab">`Age`</font>: age in months
    - <font class="vocab">`Rate`</font>: respiratory rate (breaths per minute)

---

### Log transformation on `$y$`

```r
log_model <- lm(log_rate ~ Age, data = respiratory) 
  kable(tidy(log_model, conf.int = TRUE), format = "markdown", digits = 3)
```

|term        | estimate| std.error| statistic| p.value| conf.low| conf.high|
|:-----------|--------:|---------:|---------:|-------:|--------:|---------:|
|(Intercept) |    3.845|     0.013|   304.500|       0|     3.82|     3.870|
|Age         |   -0.019|     0.001|   -25.839|       0|    -0.02|    -0.018|

- <font class = "vocab">Slope: </font> For every one month increase in Age, we expect the median respiratory rate to be multiplied by a factor of `$\exp\{-0.019\} = 1.019$` breaths per minute.

- <font class = "vocab">Intercept: </font> The expected respiratory rate for a child who is 0 months old (a newborn) is `$\exp\{3.845\} = 46.76$` breaths per minute.

---

### Confidence interval for `$\beta_j$`

- The confidence interval for the coefficient of `$x$` describing its relationship with `$\log(y)$` is

`$$\hat{\beta}_j \pm t^* SE(\hat{\beta_j})$$`

- The confidence interval for the coefficient of `$x$` describing its relationship with `$y$` is

---

### Coefficient of `Age`

`$$[\exp\{-0.02\}, \exp\{-0.018\}] = [0.981, 0.982]$$` 
]

<font class = "vocab">Interpretation: </font> We are 95% confident that for each additional month in age, we can expect the median respiratory rate to be multiplied by a factor of 0.981 to 0.982.

---

### Log Transformation on `$x$`

.pull-left[
<img src="10-model-assessment_files/figure-html/unnamed-chunk-5-1.png" style="display: block; margin: auto;" />
]

.pull-right[
<img src="10-model-assessment_files/figure-html/unnamed-chunk-6-1.png" style="display: block; margin: auto;" />
]

- Try a transformation on `$X$` if the scatterplot shows some curvature but the variance is constant for all values of `$X$`

---

### Model with Transformation on `$x$`

- <font class="vocab">Intercept: </font> When `$\log(x)=0$`, `$(x=1)$`, `$y$` is expected to be `$\beta_0$` (i.e. the mean of `$y$` is `$\beta_0$`)

- <font class="vocab">Slope: </font> When `$x$` is multiplied by a factor of `$\mathbf{C}$`, `$y$` is expected to change by `$\boldsymbol{\beta_1}\mathbf{\log(C)}$` units, i.e. the mean of `$y$` changes by `$\boldsymbol{\beta_1}\mathbf{\log(C)}$`
    - *Example*: when `$x$` is multiplied by a factor of 2, `$y$` is expected to change by `$\boldsymbol{\beta_1}\mathbf{\log(2)}$` units

---

### Rate vs. log(Age)

---

### Rate vs. Age

1. Write the equation for the model of `$y$` regressed on `$\log(x)$`.

2. Interpret the intercept in the context of the problem.

3. Interpret the slope in terms of how the mean respiratory rate changes when a child's age doubles.

4. Suppose a doctor has a patient who is currently 3 years old. Will this model provide a reliable prediction of the child's respiratory rate when her age doubles? Why or why not?
]

---

See [Log Transformations in Linear Regression](https://github.com/sta210-fa19/supplemental-notes/blob/master/log-transformations.pdf) for more details about interpreting regression models with log-transformed variables.

---

### Model Assessment & Selection
---

## Restaurant tips

What affects the amount customers tip at a restaurant?

- **Response:**
    - <font class="vocab">`Tip`</font>: amount of the tip
    
- **Predictors:**
    - <font class="vocab">`Party`</font>: number of people in the party
    - <font class="vocab">`Meal`</font>:  time of day (Lunch, Dinner, Late Night) 
    - <font class="vocab">`Age`</font>: age category of person paying the bill (Yadult, Middle, SenCit)

---

### Examining the Response Variable

---

### Response vs. Predictors

---

### Restaurant tips: model

```r
model1 <- lm(Tip ~ Party + Meal + Age , data = tips)
kable(tidy(model1),format="html",digits=3)
```

<table>
 <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;"> 1.254 </td>
   <td style="text-align:right;"> 0.394 </td>
   <td style="text-align:right;"> 3.182 </td>
   <td style="text-align:right;"> 0.002 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Party </td>
   <td style="text-align:right;"> 1.808 </td>
   <td style="text-align:right;"> 0.121 </td>
   <td style="text-align:right;"> 14.909 </td>
   <td style="text-align:right;"> 0.000 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> MealLate Night </td>
   <td style="text-align:right;"> -1.632 </td>
   <td style="text-align:right;"> 0.407 </td>
   <td style="text-align:right;"> -4.013 </td>
   <td style="text-align:right;"> 0.000 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> MealLunch </td>
   <td style="text-align:right;"> -0.612 </td>
   <td style="text-align:right;"> 0.402 </td>
   <td style="text-align:right;"> -1.523 </td>
   <td style="text-align:right;"> 0.130 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> AgeSenCit </td>
   <td style="text-align:right;"> 0.390 </td>
   <td style="text-align:right;"> 0.394 </td>
   <td style="text-align:right;"> 0.990 </td>
   <td style="text-align:right;"> 0.324 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> AgeYadult </td>
   <td style="text-align:right;"> -0.505 </td>
   <td style="text-align:right;"> 0.412 </td>
   <td style="text-align:right;"> -1.227 </td>
   <td style="text-align:right;"> 0.222 </td>
  </tr>
</tbody>
</table>
<br><br>

<center>
**Is this the best model to explain variation in Tips?**
</font>

---

## ANOVA table for regression

We can use the Analysis of Variance (ANOVA) table to decompose the variability in our response variable

|  | Sum of Squares | DF | Mean Square | F-Stat| p-value |
|------------------|----------------|--------------------|-------------|-------------|--------------------|
| Regression (Model) | `$$\sum\limits_{i=1}^{n}(\hat{y}_i - \bar{y})^2$$` | `$$p$$` | `$$\frac{MSS}{p}$$` | `$$\frac{MMS}{RMS}$$` | `$$P(F > \text{F-Stat})$$` |
| Residual | `$$\sum\limits_{i=1}^{n}(y_i - \hat{y}_i)^2$$` | `$$n-p-1$$` | `$$\frac{RSS}{n-p-1}$$` |  |  |
| Total | `$$\sum\limits_{i=1}^{n}(y_i - \bar{y})^2$$` | `$$n-1$$` | `$$\frac{TSS}{n-1}$$` |  |  |

The estimate of the regression variance, `$\hat{\sigma}^2 = RMS$`

---

## `$R^2$`

- **Recall**: `$R^2$` is the proportion of the variation in the response variable explained by the regression model
<br>

- `$R^2$` will always increase as we add more variables to the model 
  + If we add enough variables, we can always achieve `$R^2=100\%$`
<br>

- If we only use `$R^2$` to choose a best fit model, we will be prone to choose the model with the most predictor variables

---

## Adjusted `$R^2$`

- <font class="vocab">Adjusted `$R^2$`</font>: a version of `$R^2$` that penalizes for unnecessary predictor variables
<br>

- Similar to `$R^2$`, it measures the proportion of variation in the response that is explained by the regression model 
<br>

- Differs from `$R^2$` by using the mean squares rather than sums of squares and therefore adjusting for the number of predictor variables

---

## `$R^2$` and Adjusted `$R^2$`

`$$R^2 = \frac{\text{Total Sum of Squares} - \text{Residual Sum of Squares}}{\text{Total Sum of Squares}}$$`
<br>

.alert[
`$$Adj. R^2 = \frac{\text{Total Mean Square} - \text{Residual Mean Square}}{\text{Total Mean Square}}$$`
]
<br>

- `$Adj. R^2$` can be used as a quick assessment to compare the fit of multiple models; however, it should not be the only assessment!

- Use `$R^2$` when describing the relationship between the response and predictor variables

---

## Restaurant tips: ANOVA

- <font class="vocab">R output</font>

```r
kable(anova(model1), format = "markdown", digits = 3)
```

|          |  Df|   Sum Sq|  Mean Sq| F value| Pr(>F)|
|:---------|---:|--------:|--------:|-------:|------:|
|Party     |   1| 1188.636| 1188.636| 311.002|  0.000|
|Meal      |   2|   88.460|   44.230|  11.573|  0.000|
|Age       |   2|   13.032|    6.516|   1.705|  0.185|
|Residuals | 163|  622.979|    3.822|      NA|     NA|

- <font class="vocab">ANOVA table</font>

---

### Calculating `$R^2$` and Adj `$R^2$`

```r
#r-squared
mss <- 1290.12829
rss <- 622.97932	
tss <- mss + rss
(r_sq <- (tss - rss)/tss)
```

```
## [1] 0.6743626
```

```r
#adj r-squared
rms <- 3.821959	
tms <- tss/(nrow(tips)-1)
(adj_r_sq <- (tms - rms)/tms)
```

```
## [1] 0.6643738
```

---

### Restaurant tips: `$R^2$` and Adj. `$R^2$`

```r
glance(model1)
```

```
## # A tibble: 1 x 11
##   r.squared adj.r.squared sigma statistic  p.value    df logLik   AIC   BIC
##       <dbl>         <dbl> <dbl>     <dbl>    <dbl> <int>  <dbl> <dbl> <dbl>
## 1     0.674         0.664  1.95      67.5 6.14e-38     6  -350.  714.  736.
## # … with 2 more variables: deviance <dbl>, df.residual <int>
```
<br>

- Close values of `$R^2$` and Adjusted `$R^2$` indicate that the variables in the model are significant in understanding variation in the response, i.e. that there aren't a lot of unnecessary variables in the model

---

## ANOVA F Test

- Using the ANOVA table, we can test whether any variable in the model is a significant predictor of the response. We conduct this test using the following hypotheses:

.alert[
`$$\begin{aligned}&H_0: \beta_{1} =  \beta_{2} = \dots = \beta_p = 0 \\ 
&H_a: \text{at least one }\beta_j \text{ is not equal to 0}\end{aligned}$$`
]

<br>

- The statistic for this test is the `$F$` test statistic in the ANOVA table

- We calculate the p-value using an `$F$` distribution with `$p$` and `$(n-p-1)$` degrees of freedom

---

## ANOVA F Test in R

```r
model0 <- lm(Tip ~ 1, data = tips)
```

```r
model1 <- lm(Tip ~ Party + Meal + Age , data = tips)
```

```r
kable(anova(model0, model1), format="markdown", digits = 3)
```

| Res.Df|      RSS| Df| Sum of Sq|      F| Pr(>F)|
|------:|--------:|--:|---------:|------:|------:|
|    168| 1913.108| NA|        NA|     NA|     NA|
|    163|  622.979|  5|  1290.128| 67.511|      0|

**At least one coefficient is non-zero, i.e. at least one predictor in the model is significant**

---

### Testing subset of coefficients

- Sometimes we want to test whether a subset of coefficients are all equal to 0

- This is often the case when we want test 
    - whether a categorical variable with `$k$` levels is a significant predictor of the response
    - whether the interaction between a categorical and quantitative variable is significant

- To do so, we will use the  <font class="vocab3">Nested F Test </font>

---

## Nested F Test

- Suppose we have a full and reduced model:

`$$\begin{aligned}&\text{Full}: y = \beta_0 + \beta_1 x_1 + \dots + \beta_q x_q + \beta_{q+1} x_{q+1} + \dots \beta_p x_p \\
&\text{Red}: y = \beta_0 + \beta_1 x_1 + \dots + \beta_q x_q\end{aligned}$$`
<br>

- We want to test whether any of the variables `$x_{q+1}, x_{q+2}, \ldots, x_p$` are significant predictors. To do so, we will test the hypothesis:

.alert[
`$$\begin{aligned}&H_0: \beta_{q+1} =  \beta_{q+2} = \dots = \beta_p = 0 \\ 
&H_a: \text{at least one }\beta_j \text{ is not equal to 0}\end{aligned}$$`
]

---

## Nested F Test

- The test statistic for this test is

`$$F = \frac{(RSS_{reduced} - RSS_{full})\big/(p_{full} - p_{reduced})}{RSS_{full}\big/(n-p_{full}-1)}$$`
<br>

- Calculate the p-value using the F distribution with `$(p_{full} - p_{reduced})$` and `$(n-p_{full}-1)$` degrees of freedom

---

### Is `Meal` a significant predictor of tips?

<table>
 <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;"> 1.254 </td>
   <td style="text-align:right;"> 0.394 </td>
   <td style="text-align:right;"> 3.182 </td>
   <td style="text-align:right;"> 0.002 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Party </td>
   <td style="text-align:right;"> 1.808 </td>
   <td style="text-align:right;"> 0.121 </td>
   <td style="text-align:right;"> 14.909 </td>
   <td style="text-align:right;"> 0.000 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> AgeSenCit </td>
   <td style="text-align:right;"> 0.390 </td>
   <td style="text-align:right;"> 0.394 </td>
   <td style="text-align:right;"> 0.990 </td>
   <td style="text-align:right;"> 0.324 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> AgeYadult </td>
   <td style="text-align:right;"> -0.505 </td>
   <td style="text-align:right;"> 0.412 </td>
   <td style="text-align:right;"> -1.227 </td>
   <td style="text-align:right;"> 0.222 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> MealLate Night </td>
   <td style="text-align:right;"> -1.632 </td>
   <td style="text-align:right;"> 0.407 </td>
   <td style="text-align:right;"> -4.013 </td>
   <td style="text-align:right;"> 0.000 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> MealLunch </td>
   <td style="text-align:right;"> -0.612 </td>
   <td style="text-align:right;"> 0.402 </td>
   <td style="text-align:right;"> -1.523 </td>
   <td style="text-align:right;"> 0.130 </td>
  </tr>
</tbody>
</table>

---

### Tips data:  Nested F Test

`$$\begin{aligned}&H_0: \beta_{late night} = \beta_{lunch} = 0\\
&H_a: \text{ at least one }\beta_j \text{ is not equal to 0}\end{aligned}$$`

```r
reduced <- lm(Tip ~ Party + Age, data = tips)
```

```r
full <- lm(Tip ~ Party + Age + Meal, data = tips)
```

```r
kable(anova(reduced, full), format="markdown", digits = 3)
```

| Res.Df|     RSS| Df| Sum of Sq|     F| Pr(>F)|
|------:|-------:|--:|---------:|-----:|------:|
|    165| 686.444| NA|        NA|    NA|     NA|
|    163| 622.979|  2|    63.465| 8.303|      0|

**At least one coefficient associated with `Meal` is not zero. Therefore, `Meal` is a significant predictor of `Tips`.**

---

.question[
Why is it not good practice to use the individual p-values to determine a categorical variable with `$k > 2$` levels) is significant?

*Hint*: What does it actually mean if none of the `$k-1$` p-values are significant?
]

---

## Practice with Interactions

<table>
 <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;"> 1.2764989 </td>
   <td style="text-align:right;"> 0.4910882 </td>
   <td style="text-align:right;"> 2.5993270 </td>
   <td style="text-align:right;"> 0.0102086 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Party </td>
   <td style="text-align:right;"> 1.7947980 </td>
   <td style="text-align:right;"> 0.1715003 </td>
   <td style="text-align:right;"> 10.4652753 </td>
   <td style="text-align:right;"> 0.0000000 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> AgeSenCit </td>
   <td style="text-align:right;"> 0.4007889 </td>
   <td style="text-align:right;"> 0.3969295 </td>
   <td style="text-align:right;"> 1.0097230 </td>
   <td style="text-align:right;"> 0.3141431 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> AgeYadult </td>
   <td style="text-align:right;"> -0.4701634 </td>
   <td style="text-align:right;"> 0.4197146 </td>
   <td style="text-align:right;"> -1.1201978 </td>
   <td style="text-align:right;"> 0.2642977 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> MealLate Night </td>
   <td style="text-align:right;"> -1.8454674 </td>
   <td style="text-align:right;"> 0.7089728 </td>
   <td style="text-align:right;"> -2.6030159 </td>
   <td style="text-align:right;"> 0.0101039 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> MealLunch </td>
   <td style="text-align:right;"> -0.4608832 </td>
   <td style="text-align:right;"> 0.8651044 </td>
   <td style="text-align:right;"> -0.5327487 </td>
   <td style="text-align:right;"> 0.5949421 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Party:MealLate Night </td>
   <td style="text-align:right;"> 0.1108600 </td>
   <td style="text-align:right;"> 0.2846584 </td>
   <td style="text-align:right;"> 0.3894491 </td>
   <td style="text-align:right;"> 0.6974586 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Party:MealLunch </td>
   <td style="text-align:right;"> -0.0500822 </td>
   <td style="text-align:right;"> 0.2825586 </td>
   <td style="text-align:right;"> -0.1772455 </td>
   <td style="text-align:right;"> 0.8595384 </td>
  </tr>
</tbody>
</table>

.question[
1. What is the baseline level for `Meal`?
2. How do we expect the mean tips to change when `Meal == "Late Night"`, holding Age and Party constant? 
4. How does the slope of `Party` change when `Meal == "Late Night"`, holding Age and Party constant?
]
---

### Nested F test for interactions

**Are there any significant interaction effects with Party in the model?**

```r
reduced <- lm(Tip ~ Party + Age + Meal, data = tips)
```

```r
full <- lm(Tip ~ Party + Age + Meal + Age* Party + Meal * Party, 
           data = tips)
```

```r
kable(anova(reduced, full ), format="markdown", digits = 3)
```

| Res.Df|     RSS| Df| Sum of Sq|     F| Pr(>F)|
|------:|-------:|--:|---------:|-----:|------:|
|    163| 622.979| NA|        NA|    NA|     NA|
|    159| 615.380|  4|       7.6| 0.491|  0.742|

---

### Final model for now

We conclude that there are no significant interactions with `Party` in the model. Therefore, we will use the original model that only included main effects.

<table>
 <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;"> 1.254 </td>
   <td style="text-align:right;"> 0.394 </td>
   <td style="text-align:right;"> 3.182 </td>
   <td style="text-align:right;"> 0.002 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Party </td>
   <td style="text-align:right;"> 1.808 </td>
   <td style="text-align:right;"> 0.121 </td>
   <td style="text-align:right;"> 14.909 </td>
   <td style="text-align:right;"> 0.000 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> AgeSenCit </td>
   <td style="text-align:right;"> 0.390 </td>
   <td style="text-align:right;"> 0.394 </td>
   <td style="text-align:right;"> 0.990 </td>
   <td style="text-align:right;"> 0.324 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> AgeYadult </td>
   <td style="text-align:right;"> -0.505 </td>
   <td style="text-align:right;"> 0.412 </td>
   <td style="text-align:right;"> -1.227 </td>
   <td style="text-align:right;"> 0.222 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> MealLate Night </td>
   <td style="text-align:right;"> -1.632 </td>
   <td style="text-align:right;"> 0.407 </td>
   <td style="text-align:right;"> -4.013 </td>
   <td style="text-align:right;"> 0.000 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> MealLunch </td>
   <td style="text-align:right;"> -0.612 </td>
   <td style="text-align:right;"> 0.402 </td>
   <td style="text-align:right;"> -1.523 </td>
   <td style="text-align:right;"> 0.130 </td>
  </tr>
</tbody>
</table>