Regression

class: center, middle, inverse, title-slide

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

---

# Agenda

1. Recap of Movie Analysis

2. Multiple Regression

3. Next time: Categorical Predictors

---

# Recap of Movie Analysis

``` r
require(tidyverse)

mv <- read_rds('../data/mv.Rds')
```

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

- If not, Hollywood would go out of business!
  
--

- `$X$`: budget

- `$Y$`: gross

---
# Area 51 vs. EEAAO

<center><img src="figs/Area51.jpg" width = "50%"></center>
<center><img src="figs/eeaao.jpg" width = "50%"></center>
---

# Step 1: Look

``` r
summary(mv %>% select(gross,budget))
```

```
##      gross               budget         
##  Min.   :7.140e+02   Min.   :     5172  
##  1st Qu.:1.121e+07   1st Qu.: 16865322  
##  Median :5.178e+07   Median : 37212044  
##  Mean   :1.402e+08   Mean   : 57420173  
##  3rd Qu.:1.562e+08   3rd Qu.: 77844746  
##  Max.   :3.553e+09   Max.   :387367903  
##  NA's   :3668        NA's   :4482
```

---

# Step 1: Look

``` r
mv %>%
  mutate(missing = ifelse(is.na(gross) | is.na(budget),1,0)) %>%
  group_by(year) %>%
  summarise(propMissing = mean(missing)) %>%
* ggplot(aes(x = year,y = propMissing)) + geom_line()
```

---

# Some quick wrangling

``` r
mv <- mv %>%
  drop_na(gross,budget)

mv %>%
  select(gross,budget) %>%
  glimpse()
```

```
## Rows: 3,179
## Columns: 2
## $ gross  <dbl> 73677478, 53278578, 723586629, 11490339, 62…
## $ budget <dbl> 93289619, 10883789, 160147179, 6996721, 139…
```

---

# Step 2: Univariate Viz

``` r
mv %>%
  select(title,gross,budget) %>%
  gather(metric,value,-title) %>%
  ggplot(aes(x = value,color = metric)) + 
  geom_density()
```

---

# More Wrangling?

- Univariate visualization higlights significant **skew** in both measures

- Most movies don't cost a lot and don't make a lot
  
--

- But there are a few blockbusters that pull the density way out
  
--

- Let's **wrangle** two new variables that take the log of these skewed measures

- Logging transforms skewed measures to more "normal" measures

``` r
mv <- mv %>%
  mutate(gross_log = log(gross),
         budget_log = log(budget))
```

---

# Step 2: Univariate Viz

``` r
mv %>%
  select(title,gross_log,budget_log) %>%
  gather(metric,value,-title) %>%
  ggplot(aes(x = value,color = metric)) + 
  geom_density()
```

---

# Step 3: Multivariate Viz

``` r
pClean <- mv %>%
  ggplot(aes(x = budget,y = gross)) + 
  geom_point() + 
  scale_x_log10(labels = scales::dollar) + 
  scale_y_log10(labels = scales::dollar) + 
  labs(title = "Movie Costs and Returns",
       x = "Costs (logged budget)",
       y = "Returns (logged gross)")
```

---

# Step 3: Multivariate Viz

``` r
pClean + geom_smooth(method = 'lm')
```

---

# Step 4: Regression!

``` r
m <- lm(gross_log ~ budget_log,data = mv)
summary(m)
```

```
## 
## Call:
## lm(formula = gross_log ~ budget_log, data = mv)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -8.2672 -0.6354  0.1648  0.7899  8.5599 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.26107    0.30953   4.074 4.73e-05 ***
## budget_log   0.96386    0.01786  53.971  < 2e-16 ***
## ---
## Signif. codes:  
## 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.281 on 3177 degrees of freedom
## Multiple R-squared:  0.4783,	Adjusted R-squared:  0.4782 
## F-statistic:  2913 on 1 and 3177 DF,  p-value: < 2.2e-16
```

---

# Step 5.1: Univariate Viz of Errors

- Errors `$\varepsilon = Y - \hat{Y}$`

- In `R`, can also get them via `resid()` function

``` r
mv <- mv %>%
  mutate(errors_manual = gross_log - predict(m),
         errors_resid = resid(m))
```

---

# Step 5.1: Univariate Viz of Errors

``` r
mv %>%
  ggplot(aes(x = errors_resid)) + 
  geom_histogram()
```

---

# Step 5.2: Multivariate Viz of Errors

``` r
mv %>%
  ggplot(aes(x = budget_log,y = errors_resid)) + 
  geom_point() + 
  geom_smooth()
```

---

# Step 5.3: Cross Validated RMSE

``` r
set.seed(123)
rmseBudget <- NULL
for(i in 1:100) {
  inds <- sample(1:nrow(mv),size = round(nrow(mv)/2),replace = F)
  
  train <- mv %>% slice(inds)
  test <- mv %>% slice(-inds)
  
  mTrain <- lm(gross_log ~ budget_log,train)
  
  test$preds <- predict(mTrain,newdata = test)
  
  rmse <- sqrt(mean((test$gross_log - test$preds)^2,na.rm=T))
  rmseBudget <- c(rmseBudget,rmse)
}

mean(rmseBudget)
```

```
## [1] 1.279899
```

---

# Thinking like a .blue[scientist]

- Our previous model predicted `gross` as a function of `budget`

- .blue[Theoretically], is this sensible?

--
    
  1. Bigger budgets &rarr; famous actors &rarr; mass appeal &rarr; more tickets
  
--

2. Bigger budgets &rarr; advertising money &rarr; mass appeal &rarr; more tickets
  
--

- But what if the movie is just...not good?

---

# Alternative .blue[Theory]

- Good movies make more money

- .blue[Theory]: good movies &rarr; recommendations &rarr; more tickets
  
--

- Predict gross with .red[IMDB rating] (`score`)

``` r
pIMDB <- mv %>%
  ggplot(aes(x = score,y = gross)) + 
  geom_point() + 
  labs(title = "Movie gross as a function of public perception",
       x = "IMDB score",
       y = "Gross (logged)") + 
  scale_y_log10(label = scales::dollar) + 
  geom_smooth(method = 'lm',se = F)
```

---

# Alternative .red[Model]

``` r
pIMDB
```

---

# Evaluating the Model

- Let's go straight to RMSE

- We can have `R` calculate errors for us with `residuals()` command

``` r
m2 <- lm(gross_log ~ score,mv)
error <- residuals(m2)
(rmseScore <- sqrt(mean(error^2)))
```

```
## [1] 1.753146
```

- Even worse!

---

# Multivariate Regression

- Recall that we can **model** our outcome with multiple **predictors**

`$$Y = \alpha + \beta_1 X_1 + \beta_2 X_2 + \dots + \varepsilon$$`

- How much better can we predict `gross` with **BOTH** `budget` and `score`?

``` r
m3 <- lm(gross_log ~ budget_log + score,mv)
error <- residuals(m3)
(rmseBudgScore <- sqrt(mean(error^2)))
```

```
## [1] 1.248817
```

---

# Comparing Models

- Which model best predicts movie revenues?

``` r
p <- data.frame(budget = mean(rmseBudget),
           IMDB = rmseScore,
           combined = rmseBudgScore) %>%
  gather(model,rmse) %>%
  ggplot(aes(x = reorder(model,rmse),y = rmse)) + 
  geom_bar(stat = 'identity') + 
  labs(title = "Best Regression Model",
       subtitle = "Predicting a movie's revenue",
       y = "RMSE (lower is better)",
       x = "Model") + 
  coord_flip()
```

---

# Comparing Models

- Which model best predicts movie revenues?

``` r
p
```

---

# Why RMSE?

- Want to understand how good / bad our model is

- Can use it to compare models

- Do we improve our model with `score`?

---

# Why RMSE?

``` r
set.seed(123)
bsRes <- NULL
for(i in 1:100) {
  inds <- sample(1:nrow(mv),size = round(nrow(mv)/2),replace = F)
  train <- mv %>% slice(inds)
  test <- mv %>% slice(-inds)
  
  mB <- lm(gross_log ~ budget_log,train)
  mS <- lm(gross_log ~ score,train)
  mC <- lm(gross_log ~ budget_log + score,train)
  
  bsRes <- test %>%
    mutate(pB = predict(mB,newdata = test),
           pS = predict(mS,newdata = test),
           pC = predict(mC,newdata = test)) %>%
    summarise(Budget = sqrt(mean((gross_log - pB)^2,na.rm=T)),
              Score = sqrt(mean((gross_log - pS)^2,na.rm=T)),
              Combined = sqrt(mean((gross_log - pC)^2,na.rm=T))) %>%
    bind_rows(bsRes)
}
```

---

# Quick Aside: alternative `code`

- `sample_n()` and `anti_join()`

``` r
set.seed(123)
bsRes <- NULL
for(i in 1:100) {
  train <- mv %>%
    sample_n(size = round(nrow(.)*.8),replace = F)
  test <- mv %>%
*   anti_join(train)

mB <- lm(gross_log ~ budget_log,train)
  mS <- lm(gross_log ~ score,train)
  mC <- lm(gross_log ~ budget_log + score,train)
  
  bsRes <- test %>%
    mutate(pB = predict(mB,newdata = test),
           pS = predict(mS,newdata = test),
           pC = predict(mC,newdata = test)) %>%
    summarise(Budget = sqrt(mean((gross_log - pB)^2,na.rm=T)),
              Score = sqrt(mean((gross_log - pS)^2,na.rm=T)),
              Combined = sqrt(mean((gross_log - pC)^2,na.rm=T))) %>%
    bind_rows(bsRes) 
}
```

---

# Why RMSE?

``` r
bsRes %>%
  summarise_all(mean,na.rm=T)
```

```
## # A tibble: 1 × 3
##   Budget Score Combined
##    <dbl> <dbl>    <dbl>
## 1   1.29  1.76     1.26
```

---

# Visualizing

``` r
bsRes %>%
  gather(model,rmse) %>%
  ggplot(aes(x = rmse,y = reorder(model,rmse))) + 
  geom_violin()
```

---

# Visualizing Practice

- Create another plot of the RMSE for each model, but use `facet_grid()` so each RMSE score is on its own plot.

- Also, visualize the density of the RMSE using `geom_density()`.

``` r
# INSERT CODE HERE
```