econmt-diagnostics.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(regression-diagnostics)=\n",
    "# Regression diagnostics and visualisations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Introduction\n",
    "\n",
    "In this chapter, you'll see some diagnostics and visualisations that are closely tied to regressions.\n",
    "\n",
    "Most of this chapter will rely on the amazing [**statsmodels**](https://www.statsmodels.org/stable/index.html) package, including for some of the examples. Some of the material in this chapter follows [Grant McDermott](https://grantmcdermott.com/)'s excellent notes and the [Library of Statistical Translation](https://lost-stats.github.io/)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Imports\n",
    "\n",
    "Let's import some of the packages we'll be using:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import scipy.stats as st\n",
    "import statsmodels.api as sm\n",
    "import statsmodels.formula.api as smf\n",
    "from lets_plot import *\n",
    "\n",
    "LetsPlot.setup_html()\n",
    "\n",
    "# Set seed for random numbers\n",
    "seed_for_prng = 78557\n",
    "prng = np.random.default_rng(\n",
    "    seed_for_prng\n",
    ")  # prng=probabilistic random number generator"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "51a55374",
   "metadata": {
    "tags": [
     "remove-cell"
    ]
   },
   "outputs": [],
   "source": [
    "import matplotlib_inline.backend_inline\n",
    "\n",
    "# Plot settings\n",
    "plt.style.use(\n",
    "    \"https://github.com/aeturrell/coding-for-economists/raw/main/plot_style.txt\"\n",
    ")\n",
    "matplotlib_inline.backend_inline.set_matplotlib_formats(\"svg\")\n",
    "\n",
    "# Set max rows displayed for readability\n",
    "pd.set_option(\"display.max_rows\", 6)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Regression plots\n",
    "\n",
    "**statsmodels** has a number of built-in plotting methods to help you understand how well your regression is capturing the relationships you're looking for. Let's see a few examples of these using **statsmodels** built-in Statewide Crime data set:\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "crime_data = sm.datasets.statecrime.load_pandas()\n",
    "print(sm.datasets.statecrime.NOTE)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, let's look at a Q-Q plot to get a sense of how the variables are distributed. This uses **scipy**'s stats module. The default distribution is normal but you can use any that **scipy** supports."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "st.probplot(crime_data.data[\"murder\"], dist=\"norm\", plot=plt);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Clearly, this is not quite normal and there are some serious outliers in the tails.\n",
    "\n",
    "Let's run take a look at the unconditional relationship we're interested in: how murder depends on high school graduation. We'll use [**lets-plot**](https://lets-plot.org/)'s `geom_smooth()` to do this but bear in mind it will only run a linear model of `'murder ~ hs_grad'` and ignore the other covariates."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "p = (\n",
    "    ggplot(crime_data.data, aes(y=\"murder\", x=\"hs_grad\"))\n",
    "    + geom_point()\n",
    "    + geom_smooth(method=\"lm\")\n",
    ")\n",
    "p.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can take into account those other factors by using a partial regression plot that asks what does $\\mathbb{E}(y|X)$ look like as a function of $\\mathbb{E}(x_i|X)$? (Use `obs_labels=False` to remove data point labels.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with plt.rc_context({\"font.size\": 5}):\n",
    "    sm.graphics.plot_partregress(\n",
    "        endog=\"murder\",\n",
    "        exog_i=\"hs_grad\",\n",
    "        exog_others=[\"urban\", \"poverty\", \"single\"],\n",
    "        data=crime_data.data,\n",
    "        obs_labels=True,\n",
    "    )\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "At this point, the results of the regression are useful context."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "results_crime = smf.ols(\n",
    "    \"murder ~ hs_grad + urban + poverty + single\", data=crime_data.data\n",
    ").fit()\n",
    "print(results_crime.summary())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Putting the multicollinearity problems to one side, we see that the relationship shown in the partial regression plot is also implied by the coefficient on `hs_grad` in the regression table."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also look at an in-depth summary of one exogenous regressor and its relationship to the outcome variable. Each of these types of regression diagnostic are available individually, or for all regressors at once, too. The first panel is the chart we did with **lets-plot** rendered differently (and, one could argue, more informatively). Most of the plots below are self-explanatory except for the third one, the CCPR (Component-Component plus Residual) plot. This provides a way to judge the effect of one regressor on the response variable by taking into account the effects of the other independent variables. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(8, 6), dpi=150)\n",
    "\n",
    "sm.graphics.plot_regress_exog(results_crime, \"hs_grad\", fig=fig)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**statsmodels** can also produce influence plots of the 'externally studentised' residuals vs. the leverage of each observation as measured by the so-called hat matrix $X(X^{\\;\\prime}X)^{-1}X^{\\;\\prime}$ (because it puts the 'hat' on $y$). Externally studentised residuals are residuals that are scaled by their standard deviation. High leverage points could exert an undue influence over the regression line, but only if the predicted $y$ values of a regression that was fit with them excluded was quite different. In the example below, DC is having a big influence."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with plt.rc_context({\"font.size\": 6}):\n",
    "    sm.graphics.influence_plot(results_crime)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, it's nice to be able to see plots of our coefficients along with their standard errors. There isn't a built-in **statsmodels** option for this, but happily it's easy to extract the results of regressions in a sensible format. Using the `results` object from earlier, and excluding the intercept, we can get the coefficients from `results.params[1:]` and the associated errors from `results.bse[1:]`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Put the results into a dataframe with Name, Coefficient, Error\n",
    "res_df = (\n",
    "    pd.concat([results_crime.params[1:], results_crime.bse[1:]], axis=1)\n",
    "    .reset_index()\n",
    "    .rename(columns={\"index\": \"Name\", 0: \"Coefficient\", 1: \"Error\"})\n",
    ")\n",
    "res_df[\"ymax\"] = res_df[\"Coefficient\"] + res_df[\"Error\"]\n",
    "res_df[\"ymin\"] = res_df[\"Coefficient\"] - res_df[\"Error\"]\n",
    "# Plot the coefficient values and their errors\n",
    "(\n",
    "    ggplot(res_df)\n",
    "    + geom_point(aes(\"Name\", \"Coefficient\"))\n",
    "    + geom_errorbar(aes(x=\"Name\", ymin=\"ymin\", ymax=\"ymax\"))\n",
    ").show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Binscatter\n",
    "\n",
    "We're going to follow the excellent example in this [blog post](https://towardsdatascience.com/goodbye-scatterplot-welcome-binned-scatterplot-a928f67413e4) by [Matteo Courthoud](https://www.linkedin.com/in/matteo-courthoud/).\n",
    "\n",
    "First we'll generate some fake data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "def dgp_marketplace(N=10_000):\n",
    "    # Does the firm sells only online?\n",
    "    online = prng.binomial(1, 0.5, N)\n",
    "    # How many products does the firm have\n",
    "    products = 1 + prng.poisson(1, N)\n",
    "    # What is the age of the firm\n",
    "    t = prng.exponential(0.5 * products, N)\n",
    "    # Sales\n",
    "    sales = 1e3 * prng.exponential(\n",
    "        products\n",
    "        + np.maximum(\n",
    "            (1 + 0.3 * products + 4 * online) * t - 0.5 * (1 + 6 * online) * t**2, 0\n",
    "        ),\n",
    "        N,\n",
    "    )\n",
    "    # Generate the dataframe\n",
    "    df = pd.DataFrame(\n",
    "        {\"age\": t, \"sales\": sales, \"online\": online, \"products\": products}\n",
    "    )\n",
    "    return df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sales = dgp_marketplace(N=10_000)\n",
    "sales.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We have pulled down information on 10,000 firms. For each firm we know:\n",
    "\n",
    " - `age`: the age of the firm\n",
    " - `sales`: the monthly sales from last month\n",
    " - `online`: whether the firm is only active online\n",
    " - `products`: the number of products that the firm offers\n",
    "\n",
    "Suppose we are interested in understanding the relationship between age and sales. What is the life-cycle of sales?\n",
    "\n",
    "Let’s start with a simple scatterplot of sales vs age."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "(ggplot(sales, aes(x=\"age\", y=\"sales\")) + geom_point(color=\"blue\", alpha=0.4)).show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This plot is far too crowded to get a real sense of the relationship here. \n",
    "\n",
    "A linear regression can help us unpick what the relationship really is"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "smf.ols(\"np.log(sales) ~ np.log(age)\", sales).fit().summary().tables[1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There is a pretty strong positive relationship here that wasn't evident in the scatter plot. BUT it's possible that this relationships depends on whether the firms are online-only or not. We can condition on the other variables and do another regression:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "smf.ols(\n",
    "    \"np.log(sales) ~ np.log(age) + C(online) + products\", sales\n",
    ").fit().summary().tables[1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The coefficient now looks very different having conditioned on the other confounders. We can see that sales still increase with age, but we don't know if this relationship is linear or not. We could add extra terms (like age squared) to explore this.\n",
    "\n",
    "Or, we could use a *binned scatterplot* to analyse the relationship non-parametrically. The binned scatterplot divides the conditioning variable, `age` in our example, into equally sized bins or quantiles, and then plots the conditional mean of the dependent variable, `sales` in our example, within each bin.\n",
    "\n",
    "We'll be using the [**binsreg**](https://nppackages.github.io/binsreg/) package for this {cite:t}`cattaneo2019binscatter`.\n",
    "\n",
    "Binned scatterplots do not just provide conditional means for optimally chosen intervals; they can also provide inference for these means. In particular, they will give confidence intervals around each data point. In the **binsreg** package, the `ci` option adds confidence intervals to the estimation results.\n",
    "\n",
    "As the **binsreg** package is not written in the most Pythonic way, we'll adorn it a bit to make it more easily used with dataframes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import binsreg\n",
    "\n",
    "\n",
    "def binscatter(**kwargs):\n",
    "    # Estimate binsreg\n",
    "    est = binsreg.binsreg(**kwargs, noplot=True)\n",
    "    # Retrieve estimates\n",
    "    df_est = pd.concat([d.dots for d in est.data_plot])\n",
    "    df_est = df_est.rename(columns={\"x\": kwargs.get(\"x\"), \"fit\": kwargs.get(\"y\")})\n",
    "    # Add confidence intervals\n",
    "    if \"ci\" in kwargs:\n",
    "        df_est = pd.merge(df_est, pd.concat([d.ci for d in est.data_plot]))\n",
    "        df_est = df_est.drop(columns=[\"x\"])\n",
    "        df_est[\"ci\"] = df_est[\"ci_r\"] - df_est[\"ci_l\"]\n",
    "    # Rename groups\n",
    "    if \"by\" in kwargs:\n",
    "        df_est[\"group\"] = df_est[\"group\"].astype(\n",
    "            kwargs.get(\"data\")[kwargs.get(\"by\")].dtype\n",
    "        )\n",
    "        df_est = df_est.rename(columns={\"group\": kwargs.get(\"by\")})\n",
    "\n",
    "    return df_est"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Estimate binsreg\n",
    "br_est = binscatter(x=\"age\", y=\"sales\", data=sales, ci=(3, 3))\n",
    "br_est.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now plot the results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "(\n",
    "    ggplot(br_est, aes(x=\"age\", y=\"sales\", ymin=\"ci_l\", ymax=\"ci_r\"))\n",
    "    + geom_point()\n",
    "    + geom_errorbar()\n",
    "    + ggtitle(\"Binscatter: sales as a function of age\")\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now see that the relationship is non-linear with a sharp increase in `sales` at the beginning of the lifetime of a firm, followed by a plateau.\n",
    "\n",
    "As noted, this relationship may be modified by other variables though. So let's now condition on those. \n",
    "\n",
    "To condition on `products`, we pass `w=[\"products\"]` to the `binscatter()` function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Estimate binsreg\n",
    "br_est = binscatter(x=\"age\", y=\"sales\", w=[\"products\"], data=sales, ci=(3, 3))\n",
    "br_est.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's look at the plot of this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "(\n",
    "    ggplot(br_est, aes(x=\"age\", y=\"sales\", ymin=\"ci_l\", ymax=\"ci_r\"))\n",
    "    + geom_point()\n",
    "    + geom_errorbar()\n",
    "    + ggtitle(\"Binscatter: sales as a function of age conditioned on products\")\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Conditional on the number of `products`, the shape of the `sales` life-cycle changes further. Now, after an initial increase in `sales`, we observe a gradual decrease over time.\n",
    "\n",
    "Do online-only firms have different `sales` life-cycles with respect to mixed online-offline firms? We can produce different binned scatterplots by group using the option `by`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Estimate binsreg\n",
    "br_est = binscatter(\n",
    "    x=\"age\", y=\"sales\", w=[\"products\"], by=\"online\", data=sales, ci=(3, 3)\n",
    ")\n",
    "br_est[\"online\"] = br_est[\"online\"].astype(\"boolean\")\n",
    "br_est.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "(\n",
    "    ggplot(br_est, aes(x=\"age\", y=\"sales\", color=\"online\", ymin=\"ci_l\", ymax=\"ci_r\"))\n",
    "    + geom_point()\n",
    "    + geom_errorbar()\n",
    "    + ggtitle(\"Binscatter: sales as a function of age conditioned on products\")\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From the binned scatterplot, we can see that online products have on average shorter lifetimes, with a higher initial peak in sales, followed by a sharper decline.\n",
    "\n",
    "Hopefully this shows you how you can use **binsreg** to get a better understanding of causal relationships in your data.\n",
    "\n",
    "You can find out more about bin scatters in this [video](https://www.youtube.com/watch?v=fg9T2gPZCIs&ab_channel=PaulGoldsmith-Pinkham) by Paul Goldsmith-Pinkham."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Specification curve analysis\n",
    "\n",
    "When specifying a model, modellers have many options. These can be informed by field intelligence, priors, and even misguided attempts to find a significant result. Even with the best of intentions, research teams can reach entirely different conclusions using the same, or similar, data because of different choices made in preparing data or in modelling it.\n",
    "\n",
    "There’s formal evidence that researchers really do make different decisions; one study {cite:ps}`silberzahn2018many` gave the same research question - whether soccer referees are more likely to give red cards to dark-skin-toned players than to light-skin-toned players - to 29 different teams. From the abstract of that paper:\n",
    "\n",
    "> Analytic approaches varied widely across the teams, and the estimated effect sizes ranged from 0.89 to 2.93 (Mdn = 1.31) in odds-ratio units. Twenty teams (69%) found a statistically significant positive effect, and 9 teams (31%) did not observe a significant relationship. Overall, the 29 different analyses used 21 unique combinations of covariates. Neither analysts’ prior beliefs about the effect of interest nor their level of expertise readily explained the variation in the outcomes of the analyses. Peer ratings of the quality of the analyses also did not account for the variability.\n",
    "\n",
    "So not only were different decisions made, there seems to be no clearly identifiable reason for them. There is usually scope for reasonable alternative model specifications when estimating coefficients, and those coefficients will vary with those specifications. \n",
    "\n",
    "Specification curve analysis {cite:ps}`simonsohn2020specification` looks for a more exhaustive way of trying out alternative specifications. The three steps of specification curve analysis are:\n",
    "\n",
    "1. identifying the set of theoretically justified, statistically valid, and non-redundant analytic specifications;\n",
    "\n",
    "2. displaying alternative results graphically, allowing the identification of decisions producing different results; and\n",
    "\n",
    "3. conducting statistical tests to determine whether as a whole results are inconsistent with the null hypothesis.\n",
    "\n",
    "For a good example of specification curve analysis in action, see this recent Nature Human Behaviour paper {cite:ps}`orben2019association` on the association between adolescent well-being and the use of digital technology.\n",
    "\n",
    "We'll use the [**specification curve analysis**](https://aeturrell.github.io/specification_curve) package to do the first two, which you can install with `pip install specification_curve`. To demonstrate the full functionality, we'll use the \"mpg\" dataset and create a second, alternative 'hp' that is a transformed version of the original."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mpg = pd.read_csv(\n",
    "    \"https://raw.githubusercontent.com/LOST-STATS/lost-stats.github.io/source/Data/mtcars.csv\",\n",
    "    dtype={\"model\": str, \"mpg\": float, \"hp\": float, \"disp\": float, \"cyl\": \"category\"},\n",
    ")\n",
    "mpg.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mpg[\"hp_boxcox\"], _ = st.boxcox(mpg[\"hp\"])\n",
    "mpg[\"lnhp\"] = np.log(mpg[\"hp\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mpg.info()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's create a specification curve. We need to specify the data, the different outcome variables we'd like to try, `y_endog`; the different possible versions of the main regressor of interest, `x_exog`; the possible controls, `controls`; any controls that should always be included, `always_include`; and any categorical variables to include class-by-class, `cat_expand`. Some of these accept lists of variables as well as single reggressors. The point estimates that have confidence intervals which include zero are coloured in grey, instead of blue. There is also an `exclu_grps` option to exclude certain combinations of regressors, and you can pass alternative estimators to fit, for example `fit(estimator=sm.Logit)`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import specification_curve as specy\n",
    "\n",
    "sc = specy.SpecificationCurve(\n",
    "    mpg,\n",
    "    y_endog=\"mpg\",\n",
    "    x_exog=[\"lnhp\", \"hp_boxcox\"],\n",
    "    controls=[\"drat\", \"qsec\", \"cyl\", \"gear\"],\n",
    "    always_include=[\"gear\"],\n",
    "    cat_expand=\"cyl\",\n",
    ")\n",
    "sc.fit()\n",
    "sc.plot()"
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Tags",
  "kernelspec": {
   "display_name": "codeforecon",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.16"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}