Regression

class: center, middle, inverse, title-slide

.title[
# Regression
]
.subtitle[
## Part 3b
]
.author[
### Prof. Weldzius
]
.institute[
### Villanova University
]

---

# Agenda

1. Categorical Predictors

---

# Recap of Movie Analysis

``` r
require(tidyverse)
mv <- read_rds('../data/mv.Rds') %>%
  mutate(gross_log = log(gross),
         budget_log = log(budget))
```

- .blue[Theory]: the more a movie costs, the more it should make

- `\(X\)`: log(budget)

- `\(Y\)`: log(gross)

---

# Categorical Data

- Thus far, only using continuous variables

- But we can do regression with categorical data too!

- The Bechdel Test: 3 questions of a movie

1. Does it have two women in it?
  2. Who talk to each other?
  3. About something other than a man?

``` r
mv %>%
  count(bechdel_score)
```

```
## # A tibble: 5 × 2
##   bechdel_score     n
##           <dbl> <int>
## 1             0   251
## 2             1   895
## 3             2   366
## 4             3  1999
## 5            NA  4162
```

---

# Research Question

- .blue[Do movies that pass the Bechdel Test make more money?]

- .blue[Theory:] Women are ~50% of the population. Movies that pass the test are more appealing to women.
  
--

- .blue[Hypothesis:] Movies that pass the test make more money.
  
--

- .red[Wrangling:] Let's turn the `bechdel_score` variable into a binary

``` r
mv <- mv %>%
  mutate(bechdel_bin = ifelse(bechdel_score == 3,1,0)) %>%
  mutate(bechdel_factor=recode_factor(bechdel_bin,
                                      `1`="Pass",
                                      `0`="Fail",
                                      ))
```

---

# Regression

- We can add the binary factor to our regression

``` r
summary(lm(gross_log ~ bechdel_factor,mv))
```

```
## 
## Call:
## lm(formula = gross_log ~ bechdel_factor, data = mv)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -9.4532 -0.9132  0.3820  1.3392  4.2510 
## 
## Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)        17.73998    0.05278 336.105  < 2e-16 ***
## bechdel_factorFail  0.23349    0.08213   2.843  0.00451 ** 
## ---
## Signif. codes:  
## 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.98 on 2395 degrees of freedom
##   (5276 observations deleted due to missingness)
## Multiple R-squared:  0.003363,	Adjusted R-squared:  0.002947 
## F-statistic: 8.082 on 1 and 2395 DF,  p-value: 0.004507
```

---

# Regression

- Coefficient is positive

- What is the interpretation?

- Movies that fail make more money...
  
--

- ...than what?
  
--

- Movies that pass the Bechdel Test
  
--

- Categorical variables are **always interpreted in relation to the hold-out category**!

---

# Regression

- Movies that fail the test make more money!?

- **REMEMBER**: Correlation `\(\neq\)` causation

- What might explain this pattern?
  
--

- Budgets in a sexist Hollywood!
  
--

- Movies that fail the test get larger budgets
  
--

- Budgets are positively associated with gross
  
--

- So we want to "control" for budget by adding it to our regression

``` r
mBechCtrl <- lm(gross_log ~ budget_log + bechdel_factor,mv)
```

---

# Regression

``` r
summary(mBechCtrl)
```

```
## 
## Call:
## lm(formula = gross_log ~ budget_log + bechdel_factor, data = mv)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -8.6325 -0.5305  0.1287  0.6792  7.9370 
## 
## Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)         2.30814    0.34497   6.691 2.85e-11 ***
## budget_log          0.92089    0.01993  46.199  < 2e-16 ***
## bechdel_factorFail -0.18795    0.05254  -3.577 0.000356 ***
## ---
## Signif. codes:  
## 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.166 on 2055 degrees of freedom
##   (5615 observations deleted due to missingness)
## Multiple R-squared:  0.5106,	Adjusted R-squared:  0.5101 
## F-statistic:  1072 on 2 and 2055 DF,  p-value: < 2.2e-16
```

---

# Regression

- Our hypothesis is supported!

- What about non-binary categorical variables?

``` r
mv %>%
  count(rating)
```

```
## # A tibble: 13 × 2
##    rating        n
##    <chr>     <int>
##  1 Approved      1
##  2 G           153
##  3 NC-17        23
##  4 Not Rated   283
##  5 PG         1254
##  6 PG-13      2114
##  7 R          3698
##  8 TV-14         1
##  9 TV-MA         9
## 10 TV-PG         5
## 11 Unrated      52
## 12 X             3
## 13 <NA>         77
```

---

# Categorical

- Let's first remove rarely-occurring ratings

``` r
mvAnalysis <- mv %>%
  filter(!rating %in% c('Approved','TV-14','TV-MA','TV-PG','X'))
```

---

# Categorical

``` r
summary(lm(gross_log ~ rating,mvAnalysis))
```

```
## 
## Call:
## lm(formula = gross_log ~ rating, data = mvAnalysis)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -9.1487 -0.9872  0.2300  1.3192  5.5249 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      19.1299     0.2560  74.727  < 2e-16 ***
## ratingNC-17      -2.7709     0.6880  -4.028 5.74e-05 ***
## ratingNot Rated  -4.2052     0.2898 -14.511  < 2e-16 ***
## ratingPG         -0.5539     0.2708  -2.045 0.040924 *  
## ratingPG-13      -0.9188     0.2611  -3.519 0.000439 ***
## ratingR          -2.4344     0.2599  -9.367  < 2e-16 ***
## ratingUnrated    -5.8053     0.4608 -12.599  < 2e-16 ***
## ---
## Signif. codes:  
## 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.916 on 3954 degrees of freedom
##   (3693 observations deleted due to missingness)
## Multiple R-squared:  0.2274,	Adjusted R-squared:  0.2262 
## F-statistic: 193.9 on 6 and 3954 DF,  p-value: < 2.2e-16
```

---

# Categorical

- Everything makes less money than the hold-out category!

- "G"-rated movies are powered by children
  
--

- What if we wanted to compare to a different reference category?

``` r
mvAnalysis <- mvAnalysis %>%
  mutate(rating = factor(rating,
                         levels = c('R','PG-13','PG','G','Not Rated')))

mRating2 <- lm(gross_log ~ rating,mvAnalysis)
```

---

# Categorical

``` r
summary(mRating2)
```

```
## 
## Call:
## lm(formula = gross_log ~ rating, data = mvAnalysis)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -9.1487 -0.9801  0.2305  1.3185  5.5249 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     16.69544    0.04491 371.793   <2e-16 ***
## ratingPG-13      1.51563    0.06839  22.162   <2e-16 ***
## ratingPG         1.88057    0.09925  18.948   <2e-16 ***
## ratingG          2.43445    0.26005   9.362   <2e-16 ***
## ratingNot Rated -1.77073    0.14310 -12.374   <2e-16 ***
## ---
## Signif. codes:  
## 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.917 on 3922 degrees of freedom
##   (3727 observations deleted due to missingness)
## Multiple R-squared:  0.2108,	Adjusted R-squared:   0.21 
## F-statistic: 261.9 on 4 and 3922 DF,  p-value: < 2.2e-16
```

---

# Cross Validation

- This is why `sample_n()` is useful

``` r
set.seed(123)
rmseRes_rating <- NULL
for(i in 1:100) {
  train <- mvAnalysis %>%
    group_by(rating) %>%
    sample_n(size = round(n()*.8),replace = F)
  test <- mvAnalysis %>% anti_join(train)
  
  m <- lm(gross_log ~ rating,train)
  rmseRes_rating <- test %>%
    mutate(preds = predict(m,newdata = test)) %>%
    summarise(rmse = sqrt(mean((gross_log - preds)^2,na.rm=T))) %>%
    bind_rows(rmseRes_rating)
}
*rmseRes_rating %>% summarise(rmse = mean(rmse))
```

```
## # A tibble: 1 × 1
##    rmse
##   <dbl>
## 1  1.92
```

---

# Categorical Analysis Practice

- Change the hold-out category to `Not Rated`, run the regression, and cross-validate the results.

``` r
# INSERT CODE HERE
```