Simple Linear Regression - OLS.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "library(tidyverse)\n",
    "library(magrittr) ## so we can use more pipe operators, such as the pipe with assignment  %<>% \n",
    "library(olsrr) ## model diagnostics\n",
    "library(pwr) ## power analysis"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Table of Contents\n",
    "[Scatterplot](#scatter) <br>\n",
    "[Fitting a Model](#fit)<br>\n",
    "[Interpreting Model Output](#interpret)<br>\n",
    "[Post-hoc Checks](#posthoc)<br>\n",
    "    <ul>\n",
    "    <li>[Normality of Residuals](#residnorm)   </li>\n",
    "    <li>[Influential Outliers (Cook's D)](#cookd)       </li>\n",
    "    <li>[Homoscedasticity](#homo)        </li>\n",
    "    <li>[Independent Errors](#iid)         </li>\n",
    "    </ul>\n",
    "[Re-fit Model](#refit)<br>\n",
    "[Y-Intercept Issues](#yint)<br>\n",
    "    <ul>\n",
    "    <li>[Mean Centering](#meancent)     </li>\n",
    "    <li>[Standardized Coefficients](#stdized)      </li>\n",
    "    </ul>\n",
    "[Ordinal Variables](#ordinal)<br>\n",
    "[Categorical Variables](#categorical)<br>\n",
    "[Power](#power)<br>\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Simple Linear Regression\n",
    "\n",
    "As we discussed in class, we want to create a liner model that describes the relationship between a numerical IV (\"predictor\") and a numerical DV (\"outcome\").  The goal of this analysis is to describe the relationship between the two variables.  To describe this relationship we calculate the slope and intercept of the line that \"best fits\" the combined distribution of the data by minimizing the sum of squares of the residuals (the error remaining that is not explained by the model).  \n",
    "\n",
    "For our first example we’ll look at a small data set in which we’re interested in predicting the crime rate per million population based on socio-economic and demographic information at the state level.\n",
    "\n",
    "Let’s first import the data set and see what we’re working with."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## read table from an exernal website given the url.  This data is tab separated values (\"\\t\" denotes a tab) \n",
    "## and there is a header present (first row) that we want to use as our variable names\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\")\n",
    "glimpse(crime) ## inspect the data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Variable Descriptions\n",
    "\n",
    "I updated the included variable names to more descriptive names in the code above.  These are longer descriptions of the values of each of the variables.\n",
    "\n",
    "- crime.per.million: Crime rate: # of offenses reported to police per million population\n",
    "- young.males: The number of males of age 14-24 per 1000 population\n",
    "- is.south: Indicator variable for Southern states (0 = No, 1 = Yes)\n",
    "- average.ed: Mean # of years of schooling x 10 for persons of age 25 or older\n",
    "- exp.per.cap.1960: 1960 per capita expenditure on police by state and local government\n",
    "- exp.per.cap.1959: 1959 per capita expenditure on police by state and local government\n",
    "- labor.part: Labor force participation rate per 1000 civilian urban males age 14-24\n",
    "- male.per.fem: The number of males per 1000 females\n",
    "- population: State population size in hundred thousands\n",
    "- nonwhite: The number of non-whites per 1000 population\n",
    "- unemp.youth: Unemployment rate of urban males per 1000 of age 14-24\n",
    "- unemp.adult: Unemployment rate of urban males per 1000 of age 35-39\n",
    "- median.assets: Median value of transferable goods and assets or family income in tens of dollars\n",
    "- num.low.salary: The number of families per 1000 earning below 1/2 the median income\n",
    "\n",
    "<a id=\"scatter\"></a>\n",
    "## Scatterplot Matrix\n",
    "\n",
    "I first want to assess the basic relationship between the variables - especially between crime.per.million, our outcome variable, with each of the possible predictors.  For this I will use the `pairs()` function to generate a scatterplot matrix.  Because there are so many variables, I will print two matrices, including our predictor in both."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## print the matrix for the first 7 variables (columns)\n",
    "pairs(crime[1:7])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Because we know we are using the `crime.per.million` variable, we want to first focus on the relationship between that variable and the other variables in the dataset. For that we can focus on either the first row or first column of the matrix (they are identical).  While relationships between possible IVs will be an important factor in multiple regression, it is not important now.  \n",
    "\n",
    "In this set it looks like `exp.per.cap.1959` and `exp.per.cap.1960` have the strongest relationship with `crime.per.million` but let's review the rest of the variables"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## print the matrix for the crime.per.million variable (col 1) and the rest of the variables(columns) in the dataframe\n",
    "pairs(crime[c(1,8:14)])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here it appears that the variable with possibly the strongest relationship with `crime.per.million` is `median.assets`. \n",
    "\n",
    "We could look at the correlation between the variables too to assess possible relationships."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## we can use the same cor() function we did, except instead of providing two variables \n",
    "## we instead call the function on several variables. This provides a matrix of correlations.\n",
    "\n",
    "## Again, I'm going to run half the dataset so it's easier to read.\n",
    "\n",
    "cor(crime[1:7])\n",
    "cor(crime[c(1, 8:14)])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Both graphically and by correlation it looks like `exp.per.cap.1960`, the expenditure on state and local government in 1960, is the strongest predictor of crime.  This seems odd - it's a positive correlation, which means the more money that the state and local government has (and therefore the police jurisdictions have) the *higher* the number of crimes reported.  Can you think of a reason this might make sense?\n",
    "<a id=\"fit\"></a>\n",
    "## Fitting a model\n",
    "At this point we've decided to use the `exp.per.cap.1960` variable as our predictor, so we're all ready to run our model.  Remember the process of deciding the model we want to run is \"specifying the model.\"  We can report our model specification as:\n",
    "\n",
    "`crime.per.million` = $\\beta_0 + \\beta_1$`exp.per.cap.1960`\n",
    "\n",
    "By running the `lm()` function, we fit the linear model that tells us our $\\beta_0$ and $\\beta_1$ to fill into this formula, which defines the line that reflects the relationship between our IV and DV."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## run the regression using lm(DV ~ IV, data = df)\n",
    "\n",
    "mod1 <- lm(crime.per.million ~ exp.per.cap.1960, data=crime) ## here I've saved the resulting model to \"mod1\"\n",
    "summary(mod1) ## print the full summary of the model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id=\"interpret\"></a>\n",
    "## Interpreting model output\n",
    "We have a few things to look at here.  \n",
    "\n",
    "1. Our $\\beta_0$ coefficient. `(Intercept)` is self explanatory - it's the y-intercept of our line.  When `exp.per.cap.1960` equals 0 then `crime.per.million` is predicted to be 14.4.  The p-value (`Pr(>|t|)`) of the intercept is not important.  We'll return to the y-intercept later -- does it make sense for our IV to ever be zero?\n",
    "\n",
    "2. Our more important coefficient is the second one.  The row labeled `exp.per.cap.1960` gives us the information about our $\\beta_1$ or our \"coefficient with respect to expenses per capita in 1960.\"  This coefficient value (the `Estimate`) is the slope of our line.  We interpret this by saying __\"A one unit increase in expenditures per capita is associated with a 0.89 unit increase in crimes reported per million residents.\"__ I'm assuming that expenditures per capita is a dollar amount spent on state/local government per resident and crimes reported per million is the number crimes reported for every million residents.  So this means __\"For each dollar increase in expenditures per capita there is an increase of .89 crimes reported per million residents.\"__  The effect size of this variable is an unstandardized one, which means that it is in the units of the variable.  We have to decide ourselves if an increase of almost 1 crime reported (0.89) per million residents is substantively significant.\"\n",
    "\n",
    "3. The p-value of our second coefficient.  Unlike the intercept, we care about the statistical significance of our $\\beta_1$ / slope coefficient.  The null/alternative hypotheses associated with this test of significance of the coefficient are:\n",
    "\n",
    "$H_0: \\beta_1 = 0$ -- The coefficient is not significantly different from zero. <br>\n",
    "$H_A: \\beta_1 \\neq 0$ -- The coefficient is significantly different from zero.\n",
    "\n",
    "In this case the p-value (9.34 x 10^8) is much smaller than an alpha of 0.05, so we reject the null hypothesis.  This means that expenditures per capita is a significant predictor of crimes reported per million people.\n",
    "\n",
    "4. R-squared.  We will look at the Adjusted R-squared.  Recall that the adjustment is because the sample R-squared (in this output indicated as Multiple R-squared) is an inflated estimate of the population R-squared.  The adjustment allows us to make inference about the population without this inflation.  Our value here of 0.4611 indicates that `exp.per.cap.1960` accounts for 46% of the variance in `crimes.per.million`.  What is a \"good\" value of R-squared is varies by discipline, but this indicates a moderately strong relationship between our DV and IV.\n",
    "\n",
    "5. The F-statistic, or more specifically the p-value.  This test is our test of the value of our overall model.  The model we fit is compared to a \"null\" or \"empty\" model to determine if our model is a significant improvement over no model.  Our p-value here (9.34 x 10^8) is significantly lower than an alpha of 0.05, therefore our model is worth fitting.\n",
    "\n",
    "Overall, we can now define our model as:\n",
    "\n",
    "## `crime.per.million` = $14.4 + 0.89$`exp.per.cap.1960`\n",
    "\n",
    "We can theoretically use this relationship to predict values of `crime.per.million` for different states based on our knowing the value of `exp.per.cap.1960` of that state.  For example, if the value of `exp.per.cap.1960` of \"New State\" was 15.50, we could calculate the crimes reported per million residents with the linear equation:\n",
    "\n",
    "### `crime.per.million` of \"New State\" = $14.4 + 0.89(15.50) = 28.195$\n",
    "\n",
    "So \"New State\" would be predicted to have 28 crimes reported per million residents.\n",
    "\n",
    "Let's lay this line over a scatterplot of these two variables to visualize this linear relationship.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "options(repr.plot.width=5, repr.plot.height=4) ## plot size options for Jupyter notebook ONLY\n",
    "\n",
    "## use ggplot to create a scatterplot with a \"smoothed\" line defined by the lm.  We can show the se, or CI of our line as well.\n",
    "crime %>% ggplot(aes(x=exp.per.cap.1960, y=crime.per.million)) + ## indicate df, x and y variables.\n",
    "  geom_point()+ ## scatterplot\n",
    "  geom_smooth(method=lm, se=TRUE) + ## method is lm, show se.\n",
    "  expand_limits(x = 0, y = 0) + ## force ggplot to show us the y-intercept\n",
    "  geom_abline(slope = 0.89, intercept = 14.4) ## extend the line beyond where the datapoints are"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id=\"posthoc\"></a>\n",
    "## Post-hoc Checks\n",
    "We have a number of assumptions to check after running our model.  \n",
    "<a id=\"residnorm\"></a>\n",
    "### Normality of the Residuals\n",
    "The first is checking that the errors (or residuals) are normally distributed.  We do this by inspecting a QQ plot."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## plot() is used with our saved lm output.  the argument which allows us to print only the plot we specify, \n",
    "## the second plot is our QQ plot.\n",
    "plot(mod1, which = 2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can see that the residuals are normal, except for a moderate deviation from normality in the lower tail.\n",
    "<a id=\"cookd\"></a>\n",
    "### No influential outliers\n",
    "Next we need to determine if we have influential outliers.  We will look at the leverage vs. residuals plot to see if any of our observations exert undue influence on our model.  We will look for any observations fall outside of those dotted lines that indicate Cook's D."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot(mod1, which = 5) ## plot 5 is Residuals vs. Leverage"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It looks like our observation in row 29 is an extreme outlier and has inordinate leverage on our model.  We will probably want to exclude this observation and refit the model.  Observations 4 and 26 appear to also be outliers but are not past the red dotted line that indicates Cook's distance.  We can also calculate Cook's distance for our observations.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cooks.distance(mod1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The \"rule of thumb\" cutoff for influential observations is $4/n$, in this case $4/47$ or 0.0851, so we will only print those value that exceed this cutoff."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cd <- cooks.distance(mod1)\n",
    "lev <- cd > 0.0851\n",
    "cd[lev]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We have 5 potentially influential outliers, including the initial observation indicated - 29.  Note, this is the rowID of the observation, not the value.  An alternate way of visualizing potential outliers or observations with leverage is through a more detailed diagnostic plot available in the `olsrr` package."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ols_plot_resid_lev(mod1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Observations 4 and 26 have high leverage, but are not outliers.  Observations 2, 19, and 46 are potential outliers, but do not have leverage.  Only observation 29 is both an outlier and has high leverage.\n",
    "<a id=\"homo\"></a>\n",
    "### Homoscedasticity\n",
    "We have the assumption of homoscedasticity, or constant variance.  If our variance is non-constant (heteroscedastic) we will see a shape pattern in our plot of residuals vs. fitted (predicted) values.  Let's check the plot."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot(mod1, which = 1) ## residuals vs. fitted is plot #1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It does appear that there is a slight cone shape - the points are farther apart as x gets higher.  We can also conduct a Breusch Pagan Test of Heteroskedasticity (alternate spelling, also correct), using the `olsrr` package."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ols_test_breusch_pagan(mod1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here our p-value is below alpha, and therefore we have to reject the null hypothesis that our variance is constant, which means we have heteroscedasticity.  This means we have less precise estimates and biased significant tests.  There are ways of calculating \"robust\" standard errors that account for this, but is beyond the scope of this class.  So we will accept this as a limitation of our model fit.\n",
    "<a id=\"iid\"></a>\n",
    "### Errors are independent\n",
    "The final assumption we want to check is that the errors __*only*__ capture the natural variation not explained by our IV(s). Again we return to the plot of residuals vs. fitted."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# reprint resid vs. fitted.\n",
    "plot(mod1, which = 1) ## residuals vs. fitted is plot #1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There doesn't appear to be an apparent \"shape\" to the residuals that would indicate that our errors are not i.d.d. (Independent and identically distributed).\n",
    "<a id=\"refit\"></a>\n",
    "## Re-fit the Model\n",
    "Given everything we saw in our post-hoc review of assumptions, we have decided to remove our influential outlier, observation 29.  I will first remove that observation and then refit the model with the new dataframe."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dim(crime) ## what are the dimensions of crime df\n",
    "\n",
    "## what's the data in observation 29?\n",
    "crime[29,]\n",
    "## what's the distribution of exp.per.cap.1960\n",
    "summary(crime$exp.per.cap.1960)\n",
    "\n",
    "crime2 <- crime[-29,] ## save all of crime BUT observation 29 to a new df, crime2\n",
    "dim(crime2) ## what are the dimensions of crime2 df?  Should be one less row."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mod2 <- lm(crime.per.million ~ exp.per.cap.1960, data=crime2) ## same lm code, now with crime2 df\n",
    "summary(mod2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Our coefficient for `exp.per.cap.1960` changed slightly.  Now an increase of 1 dollar expenditure per capita is associated with an increase of 1.03 crimes reported per million people.  The coefficient is still significantly different from zero.  Our R-squared went up, now `exp.per.cap.1960` explains 52% of the variance in `crime.per.million`.  Keep in mind you will need to include in your report of the results that you removed an influential outlier and justify that removal.  If it wasn't a misreported value, you have potentially changed your model to fit your sample better, but potentially in the direction of not fitting your population as well.  If we had more observations with higher `exp.per.cap.1960` around 166, we could better fit the model and determine if this value is truly an influential outlier.  It all comes down to having a large enough sample size."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id=\"yint\"></a>\n",
    "## Issue: x = 0 is meaningless?\n",
    "It's perfectly fine to fit a model where the value of IV is not expected to ever realistically be 0, it just means our intercept is meaningless.  However, if we want our intercept to have a meaning, we could \"center\" our IV so that x = 0 has meaning.  \n",
    "<a id=\"meancent\"></a>\n",
    "### Mean-centering\n",
    "One way to center our IV is to adjust it so that x = 0 is the mean of the IV, and the values of x are the deviations from that mean.  This means that x will now have both positive and negative values.  The easiest way to mean center our IV is by using the scale() function.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## mean center IV using the scale() function with the argument scale = FALSE so that we're not adjusting the variance\n",
    "crime2 %<>% mutate(exp.1960.center = scale(exp.per.cap.1960, scale = FALSE))\n",
    "summary(crime2$exp.1960.center)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## refit model using centered IV\n",
    "mod3 <- lm(crime.per.million ~ exp.1960.center, data=crime2) ## same lm code, now with mean centered IV\n",
    "summary(mod3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice that the only thing that changed was our intercept, our slope, R-squared, and F-test are all identical.  Now our intercept reflects the mean of y (or the predicted value of y) when x is the mean of x.  So when the expenditure per capita is average, there are about 90 crimes reported per million people."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id=\"stdized\"></a>\n",
    "### Standardized Coefficients\n",
    "Another way to solve this problem is by calculating standardized coefficients.  This means both your x and y are standardized to have a mean of 0 and a variance of 1.  This complicates the interpretation of the model somewhat, but lets at least see how it's done."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## here I'm doing 2 things:\n",
    "## 1. using transmute so that crime3 only has the two variables I need calculated below\n",
    "## 2. not using scale = FALSE so that scale centers the mean AND makes the variance = 1\n",
    "\n",
    "crime3 <- crime2  %>% transmute(crimepm_std = scale(crime.per.million), exp1960_std = scale(exp.per.cap.1960))\n",
    "\n",
    "summary(crime3) ## look at dist of our variables"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## refit model with stdized data.\n",
    "\n",
    "options(scipen = 999) ## disable scientific notation\n",
    "\n",
    "mod4 <- lm(crimepm_std ~ exp1960_std, data=crime3) ## same lm code, now with crime3 df\n",
    "summary(mod4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So the only thing that has changed here is both of our coefficients.  They are now in terms of our standardized coefficients.  This complicates the interpretations a bit, since the units of our variables are now standard deviations, not the original units (crimes and dollars).  This is helpful if you're reporting the influence of one variable on another variable with uncommon units that are not understandable to the reader, however our interpretation of our coefficient for expenditures is now:\n",
    "\n",
    "__A one standard deviation increase in expenditure per capita is associated with a 0.73 standard deviation increase in crimes per million.__\n",
    "\n",
    "Our intercept is essentially zero, even though it doesn't appear as exactly zero in the output above.\n",
    "\n",
    "It's equivalent to talking about our variables in terms of z-scores.  This may be helpful in multiple regression when we want to talk about the relative importance of two different predictors that in an unstandardized format are measured in different units.\n",
    "<a id=\"ordinal\"></a>\n",
    "## Ordinal Variables\n",
    "We can use ordinal data as an IV (or maybe a DV) in our model, however only if it makes sense and the interpretation is meaningful. \n",
    "\n",
    "For this example we'll use the `mpg` example dataset built into R.  We'll look at two variables we've looked at before `cyl`: number of cylinders in the car, and `hwy`: highway mpg.  Cylinder has only four possible values - 4, 5, 6, or 8.  Technically it's ordinal."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "head(mpg)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "summary(lm(hwy ~ cyl, data = mpg))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's interpret this, shall we?\n",
    "\n",
    "1. Our intercept is 40.  This is the mpg when number of cylinders is 0.  Does that make sense?\n",
    "2. Our coefficient for `cyl` is -2.81.  The interpretation of that is __For every additional cylinder, there is a decrease of 2.8 miles per gallon.__ Are there such a thing as 7 cylinder cars?  \n",
    "3. Our R-squared is .57, which is pretty good.\n",
    "4. Our f-test indicates that our model fits the data well - it's better than a null model.  However, that doesn't mean that the model makes sense conceptually.\n",
    "\n",
    "Just because this model doesn't make sense doesn't mean all ordinal variables are bad and can't be used.  Just be careful to think through if it will be meaningful and appropriate given your variable."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id=\"categorical\"></a>\n",
    "## Categorical Variables\n",
    "It's probably more appropriate to treat `cyl` as a categorical variable.  Can we use categorical variables in a linear regression given our variables are supposed to be numerical?  Yes, but it requires some manipulation of the variable and some re-thinking of the interpretation of the model.\n",
    "\n",
    "We need our categorical variable to be a number.  How would we do this, even if it isn't ordinal?  We would create a series of \"dummy variables\" - 0/1 variables that indicate whether each observation is a member of that level.  So for cyl we would potentially have 4 dummy variables - cyl4, cyl5, cyl6, and cyl8 - with values of 1 if that car has that number of cylinders, and 0 otherwise."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### create dummies by hand\n",
    "\n",
    "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",
    "head(mpg2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "summary(lm(hwy ~ ., data = mpg2)) ## when we use . as the IV (after the ~) it fits the model using all of the columns\n",
    "                                  ## other than the DV specified before the ~\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now get an intercept and one row for each dummy variable.  Notice we get no coefficient for cyl8.  Why is that?  The pesky multicollinearity problem.  Basically all the information in cyl 8 can be described by the combination of the other three dummy variables.  So when we work with a categorical variable in a regression we need to \"leave one out\" or specify one level of the variable to be the reference value, to which each of the other levels are compared.  \n",
    " \n",
    "So let's interpret this output.\n",
    "\n",
    "1. The intercept is the mean of the reference group.  So here it's the mean of hwy among the cars with 8 cylinders."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mpg  %>% group_by(cyl)  %>%  summarize(mean_mpg = mean(hwy))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The other coefficients are the differences between that mean of the reference group (8 cylinders) and the mean of the group indicated.  So the coefficient for `cyl4` indicates that we would predict the `hwy` variable to equal 17.63 + 11.17 if the car has 4 cylinders.  This fitted / predicted value is the same thing as the mean of that group.\n",
    "\n",
    "2. All of the coefficients are significant.  This means that they are significantly different from zero, which in this case means that each of the group means is significantly different from the reference group - cars with 8 cylinders.\n",
    "\n",
    "3. The adjusted R-squared indicates that number of cylinders explains 58% of the variance in hwy mpg.  This means that `cyl` is a very good predictor of `hwy`.  Also it is not the same thing as the R-squared when we treated `cyl` as a number.\n",
    "\n",
    "4. The F-test indicates that our model is significantly better than a null model, which means again that `cyl` is a good predictor of `hwy`,"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Automatic Dummy Coding\n",
    "R is smart enough to automatically dummy code an IV in the model which is a factor.  By default the reference category will be the first factor."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mod5 <- lm(hwy ~ as.factor(cyl), data = mpg) ## we need to specifically convert cyl to a factor\n",
    "summary(mod5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Our model fit and R-squared is the same, because this is the exact same model, just with a different reference category.\n",
    "\n",
    "This time our intercept is the mean of the group of cars with 4 cylinders - 28.8.  Each of the other coefficients reflects the difference between the mean of that level and the reference group.  Here we see that the coefficient for 5 cylinders is not significantly different from 0 - there is not a significant difference in highway mpg between 4 and 5 cylinder cars.\n",
    "\n",
    "### Post-Hoc Plots\n",
    "\n",
    "Let's look at our plots for this model too."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "plot(mod5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We might be dealing with heteroscedasticity, but it's hard to tell given our x variable is categorical.  Our residuals are probably i.i.d. The QQ plot indicates that our residuals are fairly normal.  The residuals vs. leverage plot does not indicate that there are any influential outliers."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id=\"power\"></a>\n",
    "## Power Analysis\n",
    "Finally, we need to look at power.  For this we will use the `pwr.f2.test` function in the `pwr` package.  This function will work similarly as the others we've seen before, but our arguments will look slightly different.\n",
    "\n",
    "This is the power for the entire model - related to the omnibus F-test - and not the power of individual IV coefficients.\n",
    "\n",
    "For this function we have 5 values again:\n",
    "- u = the numerator degrees of freedom which is equal to the number of parameters in your model (not including the intercept)\n",
    "- v = n - u - 1 which is the total number of observations minus the u df minus one\n",
    "- f2 = f value computed from r-squared (see calculation below)\n",
    "- sig.level = chosen alpha\n",
    "- power = power\n",
    "\n",
    "Let's look at the first example and calculate the power of mod1."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## remind us of the summary of mod1\n",
    "summary(mod1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# u is equal to the number of parameters in model\n",
    "u = 1\n",
    "\n",
    "# v is equal to the total number of observations minus u minus 1\n",
    "v = dim(crime)[1] - u - 1\n",
    "\n",
    "rsq = summary(mod1)$adj.r.squared\n",
    "\n",
    "f = rsq/(1-rsq)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## what is the power of mod1?\n",
    "## pwr.f2.test(u, v, f2, sig.level, power)\n",
    "pwr.f2.test(u = u, v = v, f2 = f, sig.level = 0.05, power = NULL)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The power of my regression analysis is 0.999, which means I have a very small, near 0, probability of Type II error.\n",
    "\n",
    "Let's look at the power for `mod5`, my last categorical model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# remind us of summary of mod5\n",
    "summary(mod5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## how many parameters does mod 5 have?\n",
    "u5 = 3\n",
    "v5 = length(mpg$hwy) - u5 - 1\n",
    "rsq5 = summary(mod5)$adj.r.squared\n",
    "f5 = rsq5/(1-rsq5)\n",
    "## run power function\n",
    "pwr.f2.test(u = u5, v = v5, f2 = f5, sig.level = 0.05, power = NULL)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The power of `mod5` is 1, therefore there is no possibility of Type II error.\n",
    "\n",
    "What if we wanted to determine the n we needed to fit a model with a certain percentage of variance explained? Let's say we want to detect an effect of 30% of variance explained at a power of 0.8 and only one parameter/coefficient.  How do we determine n if we have no n argument?  We'll use v and work backwards."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "up = 1 # we specified one parameter\n",
    "rsqp = .3\n",
    "fp = rsqp/(1-rsqp)\n",
    "## run power function\n",
    "pwr.f2.test(u = up, v = NULL, f2 = fp, sig.level = 0.05, power = 0.8)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Rounding $v$ up to 19, we can calculate $n$ with this value, given $v = n - u - 1$ or $n = v + u + 1$.\n",
    "\n",
    "$n = 19 + 1 + 1 = 21$\n",
    "\n",
    "So we would need at least 21 observations to detect an underlying population relationship of 30% of variance in DV explained by our single IV with an alpha of 0.05 at 80% power."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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
}