class: center, middle, inverse, title-slide .title[ # Regression ] .subtitle[ ## Part 2b ] .author[ ### Prof. Weldzius ] .institute[ ### Villanova University ] --- <style type="text/css"> .small .remark-code { /*Change made here*/ font-size: 85% !important; } .tiny .remark-code { /*Change made here*/ font-size: 50% !important; } </style> # Learning Goals 1. Evaluating a regression: Univariate and multivariate visualization of errors 3. Root Mean Squared Error (RMSE) 4. Cross Validation --- # Remember where we left off last time ``` 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 ## (4494 observations deleted due to missingness) ## Multiple R-squared: 0.4783, Adjusted R-squared: 0.4782 ## F-statistic: 2913 on 1 and 3177 DF, p-value: < 2.2e-16 ``` --- # Evaluation -- - Every regression line makes mistakes - If they didn't, they wouldn't be good at **reducing complexity**! -- - How bad do ours look? - How should we begin to answer this question!? -- - Are there patterns to the mistakes? - We **overestimate** gross for movies that cost between $1m and $10m - These are the "indies" -- - We also **underestimate** gross to the "blockbusters" -- - Why? --- # Understanding Regression Lines - Regression lines choose `\(\alpha\)` and `\(\beta\)` to minimize mistakes - Mistakes (aka "errors" or "residuals") are captured in the `\(\varepsilon\)` term -- - We can apply the **process** to these! ``` r # Wrangle data to drop missingness! mv_analysis <- mv %>% drop_na(gross_log,budget_log) m <- lm(gross_log ~ budget_log,data = mv_analysis) mv_analysis$predictions <- predict(m) mv_analysis$errors <- mv_analysis$gross_log - mv_analysis$predictions summary(mv_analysis$errors) ``` ``` ## Min. 1st Qu. Median Mean 3rd Qu. Max. ## -8.2672 -0.6354 0.1648 0.0000 0.7899 8.5599 ``` --- # Univariate Viz of Errors ``` r mv_analysis %>% ggplot(aes(x = errors)) + geom_histogram() + labs(x = 'Errors: Gross - Predicted',y = 'Number of movies') ``` <img src="10b_RegressionPart2_files/figure-html/unnamed-chunk-6-1.png" style="display: block; margin: auto;" /> --- # Univariate Viz of Errors - Note that they are on average zero -- - Don't feel too proud! Mean 0 error is baked into the method - More concerned about **skew**...there is evidence of overestimating -- - Can we do more? **Conditional Analysis** - Conditional on the **predictor** (the `\(X\)` variable) --- # Multivariate Viz of Errors - Ideal is where errors are unrelated to predictor -- - This **should** appear as a rectangular cloud of points around zero ``` r pIdeal ``` <img src="10b_RegressionPart2_files/figure-html/unnamed-chunk-8-1.png" style="display: block; margin: auto;" /> --- # Multivariate Viz of Errors - This is not the case for us! ``` r mv_analysis %>% ggplot(aes(x = budget,y = errors)) + geom_smooth() + geom_point() + geom_hline(yintercept = 0,linetype = 'dashed') + scale_x_log10(label = scales::dollar) ``` <img src="10b_RegressionPart2_files/figure-html/unnamed-chunk-9-1.png" style="display: block; margin: auto;" /> --- # Multivariate Viz of Errors - Evidence of a U-shape &rarr; underpredict low and high budgets, overpredict middle budgets - Ergo, our model is **not great**! -- - Could add additional predictors `\(X_2\)`, `\(X_3\)`, etc. - Next lecture! --- # RMSE - Univariate / Multivariate visualization of errors is **important** - But we want to summarize model quality in a simpler way -- - **RMSE**: summarizes model performance with a *single number* - Useful for comparing multiple models to each other --- # RMSE - **E**rror ( `\(\varepsilon\)` ): actual outcome ( `\(Y_i\)` ) - predicted outcome ( `\(\hat{Y}_i\)` ) - The "distance" between the data and the model -- - **S**quared: `\(\varepsilon^2\)` 1. Makes all values positive 2. .blue[Exaggerates] the presence of larger errors -- - **M**ean: average these squared errors - **R**oot: take their square root (.blue[un-exaggerate]) `$$RMSE = \sqrt{\frac{1}{n}\sum_{i = 1}^n(Y_i - \hat{Y}_i)^2}$$` --- # RMSE - **.red[E]**.red[rror] ( `\({\color{red}\varepsilon}\)` ): actual outcome ( `\(Y_i\)` ) - predicted outcome ( `\(\hat{Y}_i\)` ) - The "distance" between the data and the model - **S**quared: `\(\varepsilon^2\)` 1. Makes all values positive 2. Exaggerates the presence of larger errors - **M**ean: average these squared errors - **R**oot: take their square root (un-exaggerate) `$$RMS{\color{red}E} = \sqrt{\frac{1}{n}\sum_{i = 1}^n(\underbrace{Y_i - \hat{Y}_i}_{\varepsilon})^2}$$` --- # RMSE - **E**rror ( `\(\varepsilon\)` ): actual outcome ( `\(Y_i\)` ) - predicted outcome ( `\(\hat{Y}_i\)` ) - The "distance" between the data and the model - **.red[S]**.red[quared:] `\(\varepsilon^2\)` 1. Makes all values positive 2. Exaggerates the presence of larger errors - **M**ean: average these squared errors - **R**oot: take their square root (un-exaggerate) `$$RM{\color{red}S}E = \sqrt{\frac{1}{n}\sum_{i = 1}^n(\underbrace{\varepsilon)^2}_{S}}$$` --- # RMSE - **E**rror ( `\(\varepsilon\)` ): actual outcome ( `\(Y_i\)` ) - predicted outcome ( `\(\hat{Y}_i\)` ) - The "distance" between the data and the model - **S**quared: `\(\varepsilon^2\)` 1. Makes all values positive 2. Exaggerates the presence of larger errors - **.red[M]**.red[ean: average these squared errors] - **R**oot: take their square root (un-exaggerate) `$$R{\color{red}M}SE = \sqrt{\underbrace{\frac{1}{n}\sum_{i = 1}^n}_{M}(SE)}$$` --- # RMSE - **E**rror ( `\(\varepsilon\)` ): actual outcome ( `\(Y_i\)` ) - predicted outcome ( `\(\hat{Y}_i\)` ) - The "distance" between the data and the model - **S**quared: `\(\varepsilon^2\)` 1. Makes all values positive 2. Exaggerates the presence of larger errors - **M**ean: average these squared errors - **.red[R]**.red[oot: take their square root (un-exaggerate)] `$${\color{red}R}MSE = \sqrt{(MSE)}$$` --- # RMSE - RMSE is a **single measure that summarizes model performance** -- ``` r e <- mv_analysis$gross_log - mv_analysis$predictions se <- e^2 mse <- mean(se) rmse <- sqrt(mse) # Or (rmseBudget <- sqrt(mean(mv_analysis$errors^2))) ``` ``` ## [1] 1.280835 ``` - Is this good? --- # Predicting with uncertainty - Say we're talking to investors about a new movie that costs $10m -- - How do we plug 10m into our model? ``` r summary(m)$coefficients ``` ``` ## Estimate Std. Error t value Pr(>|t|) ## (Intercept) 1.2610666 0.30952898 4.074147 4.73126e-05 ## budget_log 0.9638585 0.01785871 53.971323 0.00000e+00 ``` - `\(\hat{Y}_i = \alpha + \beta * X\)` -- - `\(\alpha = 1.26\)` and `\(\beta = 0.96\)` -- - where `\(\hat{Y}_i\)` is predicted gross (log) and `\(X\)` is $10m budget (log) -- ``` r pred_gross_log <- 1.26 + 0.96*log(1e7) ``` --- # Predicted Gross -- - Again, convert back out of logged values with `exp()` ``` r scales::dollar(exp(pred_gross_log)) ``` ``` ## [1] "$18,501,675" ``` -- - Cool! We'll make $8.5m! -- - But we know our model isn't perfect -- - Need to adjust for it's errors via **RMSE** --- # Incorporating RMSE - Simple idea: add and subtract RMSE from this prediction ``` r pred_gross_log_ub <- 1.26 + 0.96*log(1e7) + rmseBudget pred_gross_log_lb <- 1.26 + 0.96*log(1e7) - rmseBudget scales::dollar(exp(c(pred_gross_log_ub,pred_gross_log_lb))) ``` ``` ## [1] "$66,599,457" "$5,139,861" ``` -- - So we'll either make a $56m profit or we'll lose almost $5m? -- - **CONCLUSION PART 2**: maybe our model isn't very good? --- # Introducing Cross Validation - We ran a model on the full data and calculated the RMSE -- - But this approach risks .red["overfitting"] -- - .red[Overfitting] is when we get a model that happens to do well on our specific data, but isn't actually that useful for predicting elsewhere. -- - "Elsewhere": Other periods, other movies, other datasets -- - .blue[Theory:] Why care about **external validity**? -- - What is the point of measuring relationship if they don't generalize? --- # Introducing Cross Validation - In order to avoid .red[overfitting], we want to "train" our model on one part of the data, and then "test" it on a different part of the data. -- - Model "can't see" the test data &rarr; better way to evaluate performance -- - Cross Validation: randomly split our data into a train set and test set -- - *Similar to bootstrapping* --- # Introducing Cross Validation (CV) ``` r set.seed(1021) # Create list of row numbers at random inds <- sample(1:nrow(mv_analysis), size = round(nrow(mv_analysis)/2), replace = F) # Use slice(inds) to get training data train <- mv_analysis %>% slice(inds) # Use slice(-inds) to get test data test <- mv_analysis %>% slice(-inds) ``` -- - We now have two datasets of roughly the same number of observations! --- # CV to Calculate RMSE -- - We want to estimate a model based on the **test** data -- - And evaluate RMSE based on the **train** data ``` r m2 <- lm(gross_log ~ budget_log,train) # predict() function on a new dataset test$preds <- predict(m2,newdata = test) # Now calculate RMSE on the new dataset e <- test$gross_log - test$preds se <- e^2 mse <- mean(se,na.rm=T) rmse <- sqrt(mse) rmse ``` ``` ## [1] 1.28959 ``` --- # CV to Calculate RMSE - We did worse with CV! This is a *feature* -- - We are not being overconfident -- - We are avoiding "overfitting" -- - Want to do this many times (like bootstrapping) --- # CV to Calculate RMSE ``` r set.seed(123) bsRes <- NULL for(i in 1:100) { inds <- sample(1:nrow(mv_analysis), size = round(nrow(mv_analysis)/2), replace = F) train <- mv_analysis %>% slice(inds) test <- mv_analysis %>% 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)) bsRes <- c(bsRes,rmse) } mean(bsRes) ``` ``` ## [1] 1.279899 ``` --- # CV to Calculate RMSE ``` r data.frame(rmseBS = bsRes) %>% ggplot(aes(x = rmseBS)) + geom_density() + geom_vline(xintercept = mean(bsRes),linetype = 'dashed') ``` <img src="10b_RegressionPart2_files/figure-html/unnamed-chunk-19-1.png" style="display: block; margin: auto;" /> --- # Cross Validation - In this example, we used a 50-50 split -- - Often, data scientists prefer an 80-20 split -- - **Improves** the model (80% of the data is more to learn from)... - ...but still **protects** against overfitting ``` r inds <- sample(1:nrow(mv_analysis), * size = round(nrow(mv_analysis)*.8), replace = F) ```