Lecture12_MLR_PartTWO.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "library(tidyverse)\n",
    "library(olsrr) ## \"pretty\" ols plots\n",
    "library(magrittr) ## more pipe operators like  %<>% \n",
    "library(car) ## for the vif() function\n",
    "bb <- readRDS(\"bbdata.rds\")\n",
    "crime <- read.table(\"http://www.andrew.cmu.edu/user/achoulde/94842/data/crime_simple.txt\", sep = \"\\t\", header = TRUE)\n",
    "# Assign more meaningful variable names\n",
    "colnames(crime) <- c(\"crime.per.million\", \"young.males\", \"is.south\", \"average.ed\",\n",
    "                     \"exp.per.cap.1960\", \"exp.per.cap.1959\", \"labor.part\",\n",
    "                     \"male.per.fem\", \"population\", \"nonwhite\",\n",
    "                     \"unemp.youth\", \"unemp.adult\", \"median.assets\", \"num.low.salary\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Multicollinearity\n",
    "You may recall that there was an additional issue/assumption about OLS we discussed, which was not relevant with simple linear regression where we only had one IV.  This issue is multicollinearity - when two of our IVs are *too* correlated with each other.  They \"overlap\" in their prediction of the DV.\n",
    "\n",
    "#### What's the big deal?\n",
    "- redundancy in predicted variance - both IVs are trying to explain the same variation in the DV.\n",
    "- Low correlation in IVs = little to no effect to the model\n",
    "- Medium correlation in IVs = affects regression coefficients\n",
    "- High correlation = near perfect redundancy - the model \"blows up\"\n",
    "\n",
    "I've loaded the datasets we've used previously for our simple and multiple regression notebooks.  Let's look at some examples of multicollinearity.\n",
    "\n",
    "### 1. Improper dummy coding\n",
    "We saw this previously in the simple linear regression notebook - we can't include dummy codes for all levels in our model - we need to \"leave one out\" to serve as the reference variable."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mpg2 <- mpg  %>% select(cyl, hwy)  %>%  \n",
    "    mutate(cyl4 = ifelse(cyl == 4, 1, 0))  %>% \n",
    "    mutate(cyl5 = ifelse(cyl == 5, 1, 0))  %>% \n",
    "    mutate(cyl6 = ifelse(cyl == 6, 1, 0))  %>% \n",
    "    mutate(cyl8 = ifelse(cyl == 8, 1, 0))  %>% \n",
    "    select (-cyl) ## remove cyl from the new df\n",
    "summary(lm(hwy ~ ., data = mpg2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What if I instead made a mistake in my dummy coding and had overlapping groups?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mpg2 <- mpg  %>% select(cyl, hwy)  %>%  \n",
    "    mutate(cyl4 = ifelse(cyl == 4, 1, 0))  %>% \n",
    "    mutate(cyl5 = ifelse(cyl <= 5, 1, 0))  %>% ## less than or equal to 5 overlaps with cyl == 4\n",
    "    mutate(cyl6 = ifelse(cyl == 6, 1, 0))  %>% \n",
    "    select (-cyl) ## remove cyl from the new df\n",
    "summary(lm(hwy ~ ., data = mpg2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The coefficient for cyl4 is no longer significant - because cyl4 and cyl5 overlap.  This is why data cleaning is important preparation for model fitting.\n",
    "\n",
    "### 2. Including X that is computed from other Xs\n",
    "Let's look at our Big Brother data again.  We have total_hoh, total_vetos, AND total_wins.  Total wins is total comp wins, which *includes* HOH wins and veto wins.  Let's look at a small scatterplot matrix and correlation matrix."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "colnames(bb)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pairs(bb[c(\"days_in_game\", \"total_hoh\", \"total_vetos\", \"total_wins\")])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cor(bb[c(\"days_in_game\", \"total_hoh\", \"total_vetos\", \"total_wins\")])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Look at that 0.75 correlation between total_hoh and total_wins.  Let's see what bad stuff happens?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mcol <- lm(days_in_game ~ total_hoh + total_wins + total_vetos, data=bb) \n",
    "summary(mcol)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mod_w <- lm(days_in_game ~ total_hoh + total_vetos, data=bb) \n",
    "summary(mod_w)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "anova(mod_w, mcol)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3. Using more than 1 X that estimates similar/same measure\n",
    "For this one we go back to the crime data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pairs(crime[c(1, 5:6)])\n",
    "cor(crime[c(1, 5:6)])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mod_exp <- lm(crime.per.million ~ exp.per.cap.1959 + exp.per.cap.1960, data=crime) \n",
    "summary(mod_exp)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mod_exp_small <- lm(crime.per.million ~ exp.per.cap.1960, data=crime) \n",
    "summary(mod_exp_small)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "anova(mod_exp_small, mod_exp)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Because exp.per.cap.1960 and exp.per.cap.1959 are the same measure from two different years (sequential) there is little difference in the measures and they don't belong in the same model - either one will work to explain the relationship between expenditure and crimes reported."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Consequences of Multicollinearity\n",
    "\n",
    "OLS estimators are still unbiased, but\n",
    "- They have large variances and covariances, which leads to\n",
    "- Standard error increase, which leads to\n",
    "- Smaller coefficient t-scores, and\n",
    "- Confidence intervals tend to be wider, which leads to\n",
    "- Larger p-values, which leads to\n",
    "- Not rejecting the “zero null hypothesis”, which leads to\n",
    "- Greater Type II errors\n",
    "\n",
    "How can we check for multicollinearity?\n",
    "\n",
    "#### VIF (Variance Inflation Factor)\n",
    "\n",
    "Estimate of the amount of inflated variance in coefficient due to multicollinearity in the model. The higher the VIF, the less reliable the regression model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## calculating VIF\n",
    "vif(mcol)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "vif(mod_exp)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Interpreting VIF\n",
    "- Scores start at 1 – all clear\n",
    "- Scores >1 to 4: Still safe\n",
    "- Score of 5: Threshold for VIF issues (prefer 4 & under)\n",
    "- Score of 10 & up: Need to do something about it, massive multicollinearity\n",
    "\n",
    "When we reach a VIF of 4 that means the standard errors, which build the confidence intervals, would be twice as big\n",
    "as they would be otherwise."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mod_exp_vars <- lm(crime.per.million ~ exp.per.cap.1959 + exp.per.cap.1960 + median.assets + average.ed, data=crime)\n",
    "summary(mod_exp_vars)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "vif(mod_exp_vars)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### When is multicollinarity ok?\n",
    "1. control variables\n",
    "2. powers and products\n",
    "3. dummy groups with small frequencies\n",
    "\n",
    "##### So if we can fix it, why not just fix it?\n",
    "1. Std. Errors & p-values are ok before fix\n",
    "2. ‘Good’ reason to keep the offending variables (theory)\n",
    "3. Interpretive ease"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### *“Remember that all models are wrong; the practical question is how wrong do they have to be to not be useful.”*  \n",
    "– George Box\n",
    "\n",
    "<img src=\"images/box1.jpg\" width=\"400\" height=\"400\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Other error-based model evaluation measures.\n",
    "Beyond r-squared and our F-test of our model, which tells us about model fit, we can also look at a variety of error measures to evaluate our model.\n",
    "\n",
    "RSE - Residual Standard Error\n",
    "\n",
    "\n",
    "MSE - Mean Squared Error\n",
    "\n",
    "RMSE - Root Mean Square Error\n",
    "\n",
    "MAE - Mean Absolute Error\n",
    "\n",
    "RSE is based on the residual sum of squares.  It is the default error measure printed in the lm() summary.\n",
    "\n",
    "RSS = $\\sum(residuals)$ = $\\sum x_i - \\bar{x}$\n",
    "\n",
    "RSE is the $\\sqrt{RSS/df}$\n",
    "\n",
    "By taking the square root of RSS, and dividing by degrees of freedom (which is n - # of parameters - 1) we return the error to the *units* of Y, instead of squared units."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## calculate RSS\n",
    "RSS <- sum(residuals(mod_exp_vars)^2) ## mod_exp_vars is the last crime model we just fit\n",
    "RSS"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## use RSS to calculate RSE\n",
    "RSE <- sqrt(RSS / mod_exp_vars$df.residual)\n",
    "RSE"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What does this mean? It essentially works out to the standard deviation of our residuals - where RSS is the variance. \n",
    "\n",
    "General interpretation - smaller error is better.  We can compare RSE between models to see which predicts our y variable the best (which has the smallest standard deviation of the residuals).\n",
    "\n",
    "Moving to MSE - Mean Squared Error - it is the average of the residuals squared.\n",
    "\n",
    "RMSE is Root Mean Squared Error - like RSE is the square root of RSS, RMSE is the square root of MSE."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## calculate MSE\n",
    "MSE <- mean(residuals(mod_exp_vars)^2)\n",
    "MSE"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "RMSE <- sqrt(MSE)\n",
    "RMSE"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "RMSE is another measure of the error - the average magnitude of the residuals.  Because the residuals are squared before they are averaged (and then we ultimately take the square root), it can put high emphasis on larger errors.  This can be helpful when large errors are particularly concerning.\n",
    "\n",
    "And finally, MAE (Mean Absolute Error).  This time instead of squaring our residuals so that the positives don't cancel out the negatives, we use absolute value instead.  Therefore we just directly take the mean of the absolute value of the residuals - no need to square and then take the square root."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## calculate MAE\n",
    "MAE <- mean(abs(residuals(mod_exp_vars)))\n",
    "MAE"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice that MAE is lower, this is because unlike RMSE, each residual has equal weight.\n",
    "\n",
    "**In general, when interpreting any of these error measures, smaller is better.  We typically use these to compare two models, therefore look at the values of our error measurement for each model and compare.**  They are in the units of the observation, so we can sort of \"eyeball\" it and use our expert opinion to see if the difference is \"large enough\" for us to care.\n",
    "\n",
    "These measures become important when we begin to talk about model building based on reducing error."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Stepwise Selection\n",
    "Every good model starts with theory and careful review of the relevant literature, right?  If, instead, we want to fit the *best* model according to the fit statistics, we can ask R to add variables one by one until we reach some stopping threshold where adding another variable will not improve the model in a significant way.\n",
    "\n",
    "Is it a good way to build a model for the purposes of describing a phenomenon and understanding important relationships between variables? no - use theory\n",
    "\n",
    "Is it the way to build the best model for predictive purposes? yes.  The majority of machine learning supervised modeling is essentially stepwise regression for the purposes of increasing predictive power and reducing error.\n",
    "\n",
    "Stepwise regression with OLS == MACHINE LEARNING!\n",
    "\n",
    "<img src=\"images/ml.jpg\" width=\"400\" height=\"400\">\n",
    "\n",
    "Stepwise regression means adding or removing variables based on their contribution to the model until the \"perfect\" model is achieved.  This can be done by starting with an empty model and adding variables (\"forwards\"), starting with a full model and removing variables (\"backwards\") or by sequentially adding and/or removing variables (\"stepwise\").\n",
    "\n",
    "The most basic way of doing this is adding/removing variables based on p-value of the coefficient of those variables.  `olsrr` lets us do this with their variable selection functions\n",
    "\n",
    "First, we'll look at all possible models and plot the various error measures.  We'll use the crime dataset."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "names(crime)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "model <- lm(crime.per.million ~ average.ed + is.south + median.assets + \n",
    "            num.low.salary + exp.per.cap.1960 + exp.per.cap.1959 , data = crime) \n",
    "all <- ols_step_all_possible(model)\n",
    "plot(all)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "That's pretty computationally intensive if we want to look at all possible combinations of predictors, instead lets use forward regression, adding the \"best\" predictor based on the lowest p-value at each round."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# stepwise forward regression\n",
    "## remember that ~ . means include all possible IVs in the dataset\n",
    "model <- lm(crime.per.million ~ ., data = crime)\n",
    "ols_step_forward_p(model)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we instead add everything to the model and work backward based on removing the variable with the highest p-value, do we get the same model?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# stepwise backward regression\n",
    "## same model\n",
    "ols_step_backward_p(model)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We end up with 8 variables in the model, instead of the 6 in our forward regression.  We end up with essentially the same RMSE - recalling parsimony - which model should we keep?\n",
    "\n",
    "Finally, we can use stepwise regression in which we can go either direction at each step of the model fitting process."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ols_step_both_p(model)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this case we ended up with the same model as forward selection, since no variables were removed after they were added."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Overfitting\n",
    "Overfitting our models is less of a concern if our goal is simply to investigate the relationship between some variables and see if one variable is a significant predictor of the other, the way a social scientist would use a model to see which characteristics of people are related to their rating of feminists on a scale from 0-100.  If the goal is instead to predict individuals ratings of feminists, based on the characteristics of previous individuals ratings, we then have a prediction problem.  When modeling for the purpose of prediction, we want to have the model with the best predictive accuracy, but we don't want to make the model so specific to the data used to fit the model so that when we apply it to other data the predictive accuracy is poor.  By forming models based purely on model performance measures, we run the risk of \"overfitting\" models."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "R",
   "language": "R",
   "name": "ir"
  },
  "language_info": {
   "codemirror_mode": "r",
   "file_extension": ".r",
   "mimetype": "text/x-r-source",
   "name": "R",
   "pygments_lexer": "r",
   "version": "3.6.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}