Multiple_Regression_PART_THREE.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "library(tidyverse)\n",
    "library(magrittr) # for special pipe operators like  %<>% \n",
    "library(olsrr) # for post-hoc plots and functions\n",
    "library(interactions) # for interaction plots\n",
    "library(jtools) # dependency for interactions plus has summ() function for different format of printing lm output\n",
    "library(ggpubr) # for ggline plot"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Interactions\n",
    "We briefly touched on interactions when we looked at two-way ANOVA.  Interaction effects account for a third variable that affects the impact of an IV on the DV.  This could be something like the effect of a treatment differing between genders.  We would include both the type of treatment, the gender variable, and the interaction of the two in our model.\n",
    "\n",
    "## Numerical by Categorical\n",
    "The easiest type of interaction to understand and interpret is an interaction between a numerical IV and a categorical IV.  Conceptually, we are looking at how the relationship between the two numerical variables differs depending on the level of the categorical variable.\n",
    "\n",
    "We often talk about interactions as the X1 by X2 interaction (treatment by gender).  We can also talk about it as the interaction effect between X1 and X2.  The individual coefficients are the \"main effects\" and the coefficients that represent the interaction are the interaction effects.\n",
    "\n",
    "For the first example we'll use the built-in R dataset `iris`.  This data includes information about iris flowers including lengths of various parts of the flowers (petals, etc.) and the species of iris.  \n",
    "\n",
    "<img src=\"images/iris-species.png\" width=\"600\" height=\"400\">"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data(iris)\n",
    "glimpse(iris)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## first fit a model without interaction\n",
    "iris_nointeract <- lm(Petal.Length ~ Petal.Width + Species, data = iris)\n",
    "summ(iris_nointeract)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## fit a model to predict Petal.Length with Petal.Width and Species, including an interaction between width and species\n",
    "## note we use the * to indicate we want both the main effects and interaction between the two IVs\n",
    "fitiris <- lm(Petal.Length ~ Petal.Width * Species, data = iris)\n",
    "summ(fitiris)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "options(repr.plot.width=6, repr.plot.height=4) ## plot size options for Jupyter notebook ONLY\n",
    "\n",
    "## use interactions package to create interaction plot. \n",
    "\n",
    "## modx is used to indicate the variable used to determine how the lines are drawn\n",
    "## in this case since our moderating variable is a factor, we get one line for each level of the factor.\n",
    "\n",
    "interact_plot(fitiris, pred = Petal.Width, modx = Species)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice how the slopes of setosa and virginica appear to be very similar - with different y-intercepts.  Versicolor, however, has a different slope.  This means that the interaction is important to include in the model, because the relationship between petal.width and petal.length differs by species.  A tell-tale sign of an interaction effect is when the lines cross, however it is not a requirement.  Another option we can use to better visualize the relationship between the lines and the actual data we can add the observations as points on the plot."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "options(repr.plot.width=6, repr.plot.height=4) ## plot size options for Jupyter notebook ONLY\n",
    "\n",
    "## add plot.points = TRUE\n",
    "\n",
    "interact_plot(fitiris, pred = Petal.Width, modx = Species, plot.points = TRUE)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's return back to the model summary and talk about interpretation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "summ(fitiris)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The first three coefficients (after the intercept) are the main effects of Petal.Width, Speciesversicolor, and Speciesvirginica.  Speciessetosa ends up being the \"reference group\" for the species categorical variable.  \n",
    "\n",
    "We can no longer interpret our main effects in the same way.  Previously, we would say - \"Holding all else constant, an increase of one unit of Petal.Width is associated with an average increase of 0.55 units of Petal.Length.\"  We **CANNOT** say that any longer with an interaction effect in the model.  Why not?\n",
    "\n",
    "Because as Petal.Width increases, both the main effect and the interaction effect impact the prediction of Petal.Length.  Commonly, when we want to interpret these models, we look at the graph or use example data to demonstrate how the fitted values differ at different levels of the predictor variables.\n",
    "\n",
    "We should also review the significance of the various predictors.  In the first model, without the interaction effect, all the main effects were significant.  Once we add the interaction into the model, Petal.Width and Speciesversicolor are no longer significant.  This means that it is actually the interaction effect that drives most of the variance in Petal.Length, and not the main effects.  The main effect of speciesvirginica, however, actually remains significant.  If we look at the plot, this would appear to be due to the fact that compared to the reference group, setosa, the observations of virginica are much larger, although the slope of the line between those two groups remains similar (the slope influenced by Petal.Length).\n",
    "\n",
    "What happens when we add another variable, Sepal.Width?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fitiris_sepal <- lm(Petal.Length ~ Petal.Width * Species + Sepal.Length, data = iris)\n",
    "summ(fitiris_sepal)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "interact_plot(fitiris_sepal, pred = Petal.Width, modx = Species, plot.points = TRUE)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can still interpret the coefficient of Sepal.Length in the similar way - \n",
    "\n",
    "\"Holding all else constant, a one unit increase of Sepal.Length is associated with, on average, a 0.54 unit increase in Petal.Length.\"  \n",
    "\n",
    "Why can we still do this?  Because none of the other coefficients depend on Sepal.Length, therefore we can hold all else equal.\n",
    "\n",
    "## Continuous x Continuous Interaction\n",
    "We can also add interactions between two continuous variables, although they are much harder to interpret - at least the coefficients themselves.  We can say that:\n",
    "\n",
    "1. There is a significant interaction between Petal.Width and Sepal.Length.\n",
    "2. The relationship between Petal.Width and Petal.Length depends on Sepal.Length (or vice versa)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fitiris_numint <- lm(Petal.Length ~ Petal.Width * Sepal.Length, data = iris)\n",
    "summ(fitiris_numint)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "interact_plot(fitiris_numint, pred = Petal.Width, modx = Sepal.Length, plot.points = TRUE)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What this plot does is it takes three different values of Sepal.Length (the modx variable we specify) and plots the line of Petal.Width and Petal.Length.  These three values are the mean of Sepal.Length, and +/- one SD.  We can reverse it if we wish and look at the slopes of Sepal.Length at different levels of Petal.Width."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## we can use jitter and point.shape to differentiate the different observations\n",
    "\n",
    "interact_plot(fitiris_numint, modx = Petal.Width, pred = Sepal.Length, plot.points = TRUE)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Categorical x Categorical Interactions\n",
    "These are the types of interactions we looked at in Two-Way ANOVAs, but instead of using aov() we can fit the model with lm() instead.  However, it's essentially the same model, different format.\n",
    "\n",
    "For that example in the ANOVA notebook we used the warpbreaks data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data(warpbreaks)\n",
    "summary(warpbreaks)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## remind ourselves of the interaction on the mean plot\n",
    "\n",
    "warpbreaks %>% ggline(x = \"tension\", y = \"breaks\", color = \"wool\",\n",
    "       add = c(\"mean_se\", \"jitter\"),\n",
    "       palette = c(\"#00BF7D\", \"#FF61C9\"), size = 1.5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## the ANOVA we fit previously\n",
    "summary(aov(breaks ~ wool * tension, data=warpbreaks)) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## now the lm model\n",
    "wool_int <- lm(breaks ~ wool * tension, data=warpbreaks)\n",
    "summ(wool_int)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cat_plot(wool_int, modx = wool, pred = tension, plot.points = TRUE)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`interactions` will create a plot, but the ggline plot displays the interaction in a more readable way.\n",
    "\n",
    "### Final note:\n",
    "We don't have to include both of the main effects with the interaction if we don't want to have all of those parameters in the model.  Here we fit a model with wool and only the interaction of wool:tension without including the main effect of tension."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "wool_part <- lm(breaks ~ wool + wool:tension, data=warpbreaks)\n",
    "summ(wool_part)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Quadratic Terms\n",
    "Sometimes the relationship between a variable and the outcome differs by the value of that same IV.  The most commonly seen example of that is age - often the very oldest and the very youngest people will have different relationship with the outcome vs. those in the middle.  It's essentially an interaction between that variable and itself.  Often on the scatterplot we will see a parabola that indicates the quadratic relationship between the IV and DV.\n",
    "\n",
    "When we were looking at the Big Brother data, we saw a possible curvilinear relationship in the errors of the model.  Let's see if adding a quadratic term will improve the model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bb <- readRDS(\"bbdata.rds\") ## load the data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "options(repr.plot.width=4, repr.plot.height=3) ## plot size options for Jupyter notebook ONLY\n",
    "\n",
    "bb %>% ggplot(aes(x=total_hoh, y=tenure)) + ## indicate df, x and y variables.\n",
    "  geom_point()+ ## scatterplot\n",
    "  geom_smooth(method=lm, se=TRUE) ## method is lm, show CI"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## draw the lm line including the quadratic term\n",
    "\n",
    "bb %>% ggplot(aes(x=total_hoh, y=tenure)) + ## indicate df, x and y variables.\n",
    "  geom_point()+ ## scatterplot\n",
    "        stat_smooth(method = \"lm\", formula = y ~ x + I(x^2), size = 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Fit the original model\n",
    "\n",
    "mod1 <- lm(tenure ~ total_hoh, data=bb) ## here I've saved the resulting model to \"mod1\"\n",
    "summ(mod1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Fit the model with the quadratic term \n",
    "\n",
    "## We need to surround the quadratic term with I() or add the precalculated value to the dataset\n",
    "mod2 <- lm(tenure ~ total_hoh + I(total_hoh^2), data=bb)\n",
    "summ(mod2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What we see is promising, the quadratic term is significant, and our adjusted r-squared has increased.  But is this model significantly better than the first model?  We can test that."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "anova(mod1, mod2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The model with the quadratic term is a significant improvement over the model without, justifying the addition of that term into the model.  \n",
    "\n",
    "Remember the shape of the error/residuals we saw with the first model?  The benefit of adding this term is that we improve the model in a way that improves our adherence to the assumptions.  Let's compare plots!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## compare QQ plots\n",
    "options(repr.plot.width=3, repr.plot.height=3) ## plot size options for Jupyter notebook ONLY\n",
    "\n",
    "ols_plot_resid_qq(mod1)\n",
    "ols_plot_resid_qq(mod2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Adding the quadratic term improved the normality of our residuals - they are now very close to a normal distribution - the deviation in the tails has been reduced.\n",
    "\n",
    "Let's look also at the shape of our errors to see if we still see a curvilinear relationship."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## compare resid vs fit plots\n",
    "options(repr.plot.width=5, repr.plot.height=4) ## plot size options for Jupyter notebook ONLY\n",
    "\n",
    "plot(mod1, which = 1) \n",
    "plot(mod2, which = 1) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We still have heteroscedasticity, but there is no longer a curvilinear relationship between the residuals and fitted values.  We have accounted for that remaining linear relationship that caused our errors to not be independent."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Transformations\n",
    "\n",
    "Another way we can adjust our variables to improve the fit of our models and reduce the violation of model assumptions.\n",
    "\n",
    "Both the predictor (IV) and outcome (DV) variables can be transformed.  The unfortunate result of transforming the variables is the complexity of the interpretation of the model.\n",
    "\n",
    "## Transformations to the y variable\n",
    "The most common transformation to the y variable is a log transformation.  This is typically done to mitigate the skewness in the distribution of the y variable.  This skewness is commonly seen in income variables, as we've seen previously.\n",
    "\n",
    "For these transformations we'll look at the boston housing dataset you used for your homework."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "boston <- readRDS(\"boston.rds\") ## load the data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "summary(boston$medv)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## quick histogram of medv\n",
    "hist(boston$medv)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hist(log(boston$medv))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice how log transforming medv makes the distribution look more normal and removes the skewness.  We will create the log transformation of medv as a new variable, then use it in fitting our model.  We will compare the basic model where medv is predicted by lstat, to a model where lstat predicts log_medv."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "boston$log_medv <- log(boston$medv)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "normmod <- lm(medv ~ lstat, data = boston)\n",
    "logmod <- lm(log_medv ~ lstat, data = boston)\n",
    "summ(normmod)\n",
    "summ(logmod)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The first thing I notice is that lstat is still a significant predictor of medv, and it's still a negative relationship.  \n",
    "\n",
    "Comparing the r-squared values, lstat predicts about 54% of the variance in medv, but it predicts 65% of the variance in log_medv.  We cannot use an F-test to compare the fit of the models because they have two different y variables.\n",
    "\n",
    "But it becomes tricky when we get to the interpretation of the coefficient.  We cannot interpret the coefficient in the second model in direct regard to medv. \n",
    "\n",
    "We would have to say:  \n",
    "\n",
    "\"A one unit increase in lstat yields a 0.05 decrease in log(medv).\"  Given that log median value its hard to understand, we need to do something else.\n",
    "\n",
    "**Only the dependent/response variable is log-transformed:**\n",
    "Exponentiate the coefficient, subtract one from this number, and multiply by 100. This gives the percent increase (or decrease) in the response for every one-unit increase in the independent variable. Example: the coefficient is 0.198. (exp(0.198) – 1) * 100 = 21.9. For every one-unit increase in the independent variable, our dependent variable increases by about 22%.\n",
    "(https://data.library.virginia.edu/interpreting-log-transformations-in-a-linear-model/)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## interpret the coefficient of lstat\n",
    "\n",
    "(exp(-0.05) - 1) * 100"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So for every unit increase in lstat (proportion of lower status residents in the census tract) median home value in 1000s decreases by about 5%.\n",
    "\n",
    "Let's see if this transformation improves our residuals/error."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## compare QQ plots\n",
    "options(repr.plot.width=3, repr.plot.height=3) ## plot size options for Jupyter notebook ONLY\n",
    "\n",
    "ols_plot_resid_qq(normmod)\n",
    "ols_plot_resid_qq(logmod)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The log transformation improves the normality of our residuals, let's also check the plot of residuals vs. fitted."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## compare resid vs fit plots\n",
    "options(repr.plot.width=5, repr.plot.height=4) ## plot size options for Jupyter notebook ONLY\n",
    "\n",
    "plot(normmod, which = 1) \n",
    "plot(logmod, which = 1) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The \"hook\" relationship evident in the plot from the \"normal\" model is no longer as evident in the second plot, therefore indicating that the errors are more normally distributed.  \n",
    "\n",
    "## Transformations to the x variable\n",
    "Sometimes a log transformation is beneficial for an x variable where we see \"heaping\" at the lower values (skewness).  The log transformation smooths out the distribution and makes it more even.  This time we will log transform the crim variable, and use that to predict medv."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "options(repr.plot.width=4, repr.plot.height=3) ## plot size options for Jupyter notebook ONLY\n",
    "\n",
    "pairs(boston[c(1,14)])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "boston$log_crim <- log(boston$crim)\n",
    "pairs(boston[c(14,17)])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can see the improvement in the scatterplot there is more of a linear relationship between the variables vs. the \"heaping\" seen previously.  Let's fit two models, one with crim and one with log_crim as the predictor.  We will use the untransformed version of medv as the outcome variable"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "normmod_crim <- lm(medv ~ crim, data = boston)\n",
    "logmod_crim <- lm(medv ~ log_crim, data = boston)\n",
    "summ(normmod_crim)\n",
    "summ(logmod_crim)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Crime rate remains a significant predictor of medv, but again the r-squared value of the model is higher when using log transformed predictor.  Again we cannot compare the models using F-test, because they are not nested.\n",
    "\n",
    "**Only independent/predictor variable(s) is log-transformed.** Divide the coefficient by 100. This tells us that a 1% increase in the independent variable increases (or decreases) the dependent variable by (coefficient/100) units. Example: the coefficient is 0.198. 0.198/100 = 0.00198. For every 1% increase in the independent variable, our dependent variable increases by about 0.002. For x percent increase, multiply the coefficient by log(1.x). Example: For every 10% increase in the independent variable, our dependent variable increases by about 0.198 * log(1.10) = 0.02.\n",
    "\n",
    "So, in our case, for every 1% increase in crime rate, there is a -0.0193 unit (1000s of dollars) increase in medv.  Therefore medv decreases by about $19.\n",
    "\n",
    "Again, let's look at our errors\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## compare QQ plots\n",
    "options(repr.plot.width=3, repr.plot.height=3) ## plot size options for Jupyter notebook ONLY\n",
    "ols_plot_resid_qq(normmod_crim)\n",
    "ols_plot_resid_qq(logmod_crim)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The log transformation of the predictor doesn't seem to have improved the normality of the residuals, However..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## compare resid vs fit plots\n",
    "options(repr.plot.width=5, repr.plot.height=4) ## plot size options for Jupyter notebook ONLY\n",
    "\n",
    "plot(normmod_crim, which = 1) \n",
    "plot(logmod_crim, which = 1) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There is no longer a linear relationship evident in the plot of residuals vs. fitted, therefore our errors are now independent.  In addition, it appears that the variance is fairly constant in this model.  \n",
    "\n",
    "## Transforming both x and y\n",
    "We can log transform both x and y.  Let's try that using medv and crim."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pairs(boston[c(15,17)])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "normmod_both <- lm(medv ~ crim, data = boston)\n",
    "logmod_both <- lm(log_medv ~ log_crim, data = boston)\n",
    "summ(normmod_both)\n",
    "summ(logmod_both)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The r-squared doubles between the first and second model.  This model is a considerable improvement.\n",
    "\n",
    "**Both dependent/response variable and independent/predictor variable(s) are log-transformed.** Interpret the coefficient as the percent increase in the dependent variable for every 1% increase in the independent variable. Example: the coefficient is 0.198. For every 1% increase in the independent variable, our dependent variable increases by about 0.20%. For x percent increase, calculate 1.x to the power of the coefficient, subtract from 1, and multiply by 100. Example: For every 20% increase in the independent variable, our dependent variable increases by about (1.20 0.198 – 1) * 100 = 3.7 percent.\n",
    "\n",
    "So - for every 1% increase in crime rate, medv decreases by about 0.11%\n",
    "\n",
    "Again, we'll check our post-hoc plots."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## compare QQ plots\n",
    "options(repr.plot.width=3, repr.plot.height=3) ## plot size options for Jupyter notebook ONLY\n",
    "ols_plot_resid_qq(normmod_both)\n",
    "ols_plot_resid_qq(logmod_both)\n",
    "\n",
    "## compare resid vs fit plots\n",
    "options(repr.plot.width=5, repr.plot.height=4) ## plot size options for Jupyter notebook ONLY\n",
    "\n",
    "plot(normmod_both, which = 1) \n",
    "plot(logmod_both, which = 1) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "These look fairly similar to the previous set of models where only crim was log transformed, however, there may be a small improvement to the normality of the residuals.  \n",
    "\n",
    "These examples highlight the importance of working with a model to improve the fit.  When you have a theory about the relationship between variables and the predictors you believe to influence your outcome, if your model doesn't fit that well and has violation of assumptions, instead of giving up on the model or changing variables (p-hacking) we instead work to improve the model fit to refine the model.  Log-transformation is not the only transformation available.\n",
    "\n",
    "For IVs/predictors:\n",
    "\n",
    "<img src=\"images/IVtrans.jpg\" width=\"600\" height=\"400\">\n",
    "\n",
    "See the following for more information about possible transformations:\n",
    "- https://newonlinecourses.science.psu.edu/stat462/node/155/\n",
    "- Brief tutorial on Box-Cox Transformations: https://rpubs.com/bskc/288328\n",
    "- Tukey's Ladder of Powers: http://onlinestatbook.com/2/transformations/tukey.html"
   ]
  },
  {
   "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
}