{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Linear regression\n",
"\n",
"files needed = ('sleep75.dta', 'wage1.dta')\n",
"\n",
"This notebook introduces us to the statsmodels package [(docs)](https://devdocs.io/statsmodels/), which provides functions for formulating and estimating statistical models. Econometrics is a prerequisite for this course, so this notebook will not address the models, per se, but will focus on how to take what you learned in econometrics class and use it in python. \n",
"\n",
"Most of you used STATA in your econometrics course. STATA is a great package for econometrics. Python can do most of what STATA can do, but STATA will have more specialized routines available. As python's popularity is grows the kinds of models you can estimated in grows, too. \n",
"\n",
"If STATA is your thing, this [page](http://rlhick.people.wm.edu/posts/comparing-stata-and-ipython-commands-for-ols-models.html) on Rob Hicks' website is a nice STATA to python concordance. "
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"import pandas as pd # for data handling\n",
"import numpy as np # for numerical methods and data structures\n",
"import matplotlib.pyplot as plt # for plotting\n",
"import seaborn as sea # advanced plotting\n",
"\n",
"import patsy # provides a syntax for specifying models \n",
"import statsmodels.api as sm # provides statistical models like ols, gmm, anova, etc...\n",
"import statsmodels.formula.api as smf # provides a way to directly spec models from formulas"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": []
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Reading Stata data files\n",
"\n",
"Most of you used Wooldridge's textbook in econometrics. I figure you would like to relive those happy times, so we can work through some of the problems in the Wooldridge textbook. \n",
"\n",
"On the plus side, the data that correspond to the Wooldridge problems are available to download and they are **ALREADY CLEANED.** \\[I contemplated adding some junk to the files to make it more interesting...\\]\n",
"\n",
"On the minus side, the files are in STATA's .dta format. \n",
"\n",
"Lucky for use, pandas has a method that [reads stata files](http://pandas.pydata.org/pandas-docs/stable/generated/pandas.read_stata.html). It also has methods for SQL, SAS,json,..."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Use pandas read_stata method to get the stata formatted data file into a DataFrame.\n",
"sleep = pd.read_stata('sleep75.dta')\n",
"\n",
"# Take a look...so clean!\n",
"sleep.head()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"sleep.info()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Describing models with patsy\n",
"The patsy package provides us with a formulaic syntax for defining models that uses strings. The basic syntax is \n",
"```\n",
"y ~ x1 + x2\n",
"```\n",
"which describes the model \n",
"\n",
"$y = \\beta_0 + \\beta_1x_1 + \\beta_2x_2 + \\epsilon$.\n",
"\n",
"Notice that I did not specify the constant. Patsy takes care of that automatically. Let's start slow and build up the regression, then we will see how to do it all in one shot. "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We start by passing our model, specified in a patsy string to `patsy.dmatrices( )` to create the *design matrices*. Our model is \n",
"\n",
"$$ sleep = \\beta_0 + \\beta_1 totwrk + \\beta_2 educ + \\beta_3 age + \\epsilon $$.\n",
"\n",
"\\[This is in problem 3, chapter 3 or Wooldrigde.\\]"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Pass the model formula and the associated data to create design matrices\n",
"y, X = patsy.dmatrices('sleep ~ totwrk + educ + age', sleep)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# What do we have?\n",
"print('X and y are of type:' , type(X), type(y))"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Peak at X and y\n",
"X"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"y"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"So `X` holds the independent variables and `y` holds the dependent variable. Note the addition of the intercept (the column of ones) to the `X` matrix. "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Building and estimating the model in statsmodels\n",
"\n",
"With the design matrices in hand, we can build **ordinary least squares** model in statsmodels. "
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Pass design matrices to OLS to spec an ordinary least squares model\n",
"sleep_model = sm.OLS(y, X)\n",
"type(model)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We can now estimate (fit) the model using the `.fit( )` method of the statsmodel object."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"res = model.fit() # Estimate the model and store the results in res\n",
"type(res) # What do we have here?\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# To see the summary report\n",
"print(res.summary())"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The more you work, the less you sleep. I feel better already.\n",
"\n",
"We can retrieve individual results from the RegressionResultsWrapper object. Try res.`TAB` to see what's in there. "
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"print('The parameters are:', res.params, '\\n')\n",
"print('The confidence intervals are:', res.conf_int(), '\\n')\n",
"print('The r-sqared is:', res.rsquared)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Directly specifying and estimating models with the formula.api\n",
"\n",
"We have built our model from the ground up\n",
"1. Create the design matrices with patsy\n",
"2. Create the model with statsmodel\n",
"3. Fit the model and obtain results\n",
"\n",
"That was helpful to get a sense of what is going on, but we can do all those steps in one line of code. We just pass the patsy string and the data directly to statsmodels and call fit.\n",
"\n",
"To do this, we use the `statsmodels.formula.api` methods, which we imported as `smf`. The formula api interprets the patsy formula for us, and creates the design matrices. "
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"res = smf.ols('sleep ~ totwrk + educ + age', data=sleep).fit()\n",
"print(res.summary())"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Data transformations\n",
"Patsy can handle common (and less common) regression tasks. Here are a few\n",
"\n",
"#### Interacting variables\n",
"Use '\\*' to interact two variables. Patsy will also include the two variables individually, too. Below, we interact education an age\n",
"\n",
"$$ sleep = \\beta_0 + \\beta_1 totwrk + \\beta_2 educ + \\beta_3 age + \\beta_4 age\\times educ + \\epsilon $$."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"res = smf.ols('sleep ~ totwrk + educ*age', data=sleep).fit()\n",
"print(res.summary())"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### logs \n",
"\n",
"We use the `np.log( )` method directly in the patsy syntax. Note that we loaded the numpy package above as np. Below, we use the logarithm of age in the model\n",
"\n",
"$$ sleep = \\beta_0 + \\beta_1 totwrk + \\beta_2 educ + \\beta_3 \\log(age) + \\epsilon $$."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"res = smf.ols('sleep ~ totwrk + educ + np.log(age)', data=sleep).fit()\n",
"print(res.summary())"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Fixed effects\n",
"When I was a kid, we called these dummy variables. Gender is coded {0,1} in the variable male. "
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"sleep['male'].head()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"res = smf.ols('sleep ~ totwrk + educ + age + C(male)', data=sleep).fit()\n",
"print(res.summary())"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The 'T.1.0' notation is a bit confusing in this context. It means it is giving the value for the '1.0' variable. Since we have included a constant, one of the categories gets dropped. \n",
"\n",
"To see things more clearly, let's recode the male variable. "
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"sleep['gender'] = sleep['male'].replace({1.0:'male', 0.0:'female'})\n",
"\n",
"res = smf.ols('sleep ~ totwrk + educ + age + C(gender)', data=sleep).fit()\n",
"print(res.summary())"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### No intercept \n",
"Use `-1` to kill the automatic intercept. Try is with our gender data to explicitly recover the female fixed effect. Now you have to do the subtraction yourself to see that males sleep, on average, 87.9933 more hours."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"res = smf.ols('sleep ~ totwrk + educ + age + C(gender) -1', data=sleep).fit()\n",
"print(res.summary())"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Practice\n",
"\n",
"Take a few minutes and try the following. Feel free to chat with those around you if you get stuck. The TA and I are here, too. \n",
"\n",
"Wooldridge problem C2 in chapter 6. \n",
"\n",
"1. Load wage1.dta\n",
"2. Scatter plot log(wage) against educ"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"3. Estimate \n",
"$$ \\log(wage) = \\beta_0 + \\beta_1 educ + \\epsilon$$"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"4. Scatter plot the residuals against education. The residuals are in the results object from the fit method. "
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"5. Looks heteroskedastic. We can change the covariance matrix type (which will correct the standard error calculation) using the `cov_type` parameter [(docs)](https://www.statsmodels.org/dev/generated/statsmodels.regression.linear_model.OLS.fit.html#statsmodels.regression.linear_model.OLS.fit). The types of covariance matrices are described in the [docs](https://www.statsmodels.org/dev/generated/statsmodels.regression.linear_model.RegressionResults.html).\n",
"\n",
"Try 'HC3' for your covariance matrix type."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"6. Scatter plot the data and add the regression line.\n",
"\n",
"\n",
"To plot the regression line you will need to create some x data and then apply the parameters. I used something like this\n",
"```python\n",
"y = [p.Intercept + p.educ*i for i in x]\n",
"```\n",
"where `p` hold the parameters from my results and x holds a few x data points. "
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"7. Go back and add the text 'log(w) = a + b*educ' to the northwest corner of your plot. Replace a and b with the estimated parameter values. "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"8. Let's add some more regressors. Estimate\n",
"$$ \\log(wage) = \\beta_0 + \\beta_1 educ + \\beta_2 exper + \\beta_3 exper^2 + \\beta_4I_m + \\epsilon$$\n",
"\n",
"where $I_m$ is a variable equal to 1 if the worker is a married."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"celltoolbar": "Attachments",
"kernelspec": {
"display_name": "Python 3",
"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.6.5"
}
},
"nbformat": 4,
"nbformat_minor": 2
}