{"id":1449,"date":"2024-12-04T18:46:23","date_gmt":"2024-12-04T18:46:23","guid":{"rendered":"https:\/\/www.allendowney.com\/blog\/?p=1449"},"modified":"2024-12-04T18:46:23","modified_gmt":"2024-12-04T18:46:23","slug":"multiple-regression-with-statsmodels","status":"publish","type":"post","link":"https:\/\/www.allendowney.com\/blog\/2024\/12\/04\/multiple-regression-with-statsmodels\/","title":{"rendered":"Multiple Regression with StatsModels"},"content":{"rendered":"\n
\n

This is the third is a series of excerpts from Elements of Data Science<\/em> which available from Lulu.com<\/a> and online booksellers. It’s from Chapter 10, which is about multiple regression. You can read the complete chapter here<\/a>, or run the Jupyter notebook on Colab<\/a>.<\/p>\n<\/blockquote>\n\n\n\n

In the previous chapter we used simple linear regression to quantify the relationship between two variables. In this chapter we\u2019ll get farther into regression, including multiple regression and one of my all-time favorite tools, logistic regression. These tools will allow us to explore relationships among sets of variables. As an example, we will use data from the General Social Survey (GSS) to explore the relationship between education, sex, age, and income.<\/p>\n\n\n\n

The GSS dataset contains hundreds of columns. We\u2019ll work with an extract that contains just the columns we need, as we did in Chapter 8. Instructions for downloading the extract are in the notebook for this chapter.<\/p>\n\n\n\n

We can read the DataFrame<\/code> like this and display the first few rows.<\/p>\n\n\n\n

import pandas as pd\n\ngss = pd.read_hdf('gss_extract_2022.hdf', 'gss')\ngss.head()\n<\/pre>\n\n\n\n
<\/th>year<\/th>id<\/th>age<\/th>educ<\/th>degree<\/th>sex<\/th>gunlaw<\/th>grass<\/th>realinc<\/th><\/tr><\/thead>
0<\/th>1972<\/td>1<\/td>23.0<\/td>16.0<\/td>3.0<\/td>2.0<\/td>1.0<\/td>NaN<\/td>18951.0<\/td><\/tr>
1<\/th>1972<\/td>2<\/td>70.0<\/td>10.0<\/td>0.0<\/td>1.0<\/td>1.0<\/td>NaN<\/td>24366.0<\/td><\/tr>
2<\/th>1972<\/td>3<\/td>48.0<\/td>12.0<\/td>1.0<\/td>2.0<\/td>1.0<\/td>NaN<\/td>24366.0<\/td><\/tr>
3<\/th>1972<\/td>4<\/td>27.0<\/td>17.0<\/td>3.0<\/td>2.0<\/td>1.0<\/td>NaN<\/td>30458.0<\/td><\/tr>
4<\/th>1972<\/td>5<\/td>61.0<\/td>12.0<\/td>1.0<\/td>2.0<\/td>1.0<\/td>NaN<\/td>50763.0<\/td><\/tr><\/tbody><\/table><\/figure>\n\n\n\n
We\u2019ll start with a simple regression, estimating the parameters of real income as a function of years of education. First we\u2019ll select the subset of the data where both variables are valid.<\/p>\n\n\n\n
data = gss.dropna(subset=['realinc', 'educ'])\nxs = data['educ']\nys = data['realinc']\n<\/pre>\n\n\n\n
Now we can use linregress<\/code> to fit a line to the data.<\/p>\n\n\n\n
from scipy.stats import linregress\nres = linregress(xs, ys)\nres._asdict()\n<\/pre>\n\n\n\n
{'slope': 3631.0761003894995,\n 'intercept': -15007.453640508655,\n 'rvalue': 0.37169252259280877,\n 'pvalue': 0.0,\n 'stderr': 35.625290800764,\n 'intercept_stderr': 480.07467595184363}\n<\/pre>\n\n\n\n
The estimated slope is about 3450, which means that each additional year of education is associated with an additional $3450 of income.<\/p>\n\n\n\n
Regression with StatsModels<\/h2>\n\n\n\n
SciPy doesn\u2019t do multiple regression, so we\u2019ll to switch to a new library, StatsModels. Here\u2019s the import statement.<\/p>\n\n\n\n
import statsmodels.formula.api as smf\n<\/pre>\n\n\n\n
To fit a regression model, we\u2019ll use ols<\/code>, which stands for \u201cordinary least squares\u201d, another name for regression.<\/p>\n\n\n\n
results = smf.ols('realinc ~ educ', data=data).fit()\n<\/pre>\n\n\n\n
The first argument is a formula string<\/strong> that specifies that we want to regress income as a function of education. The second argument is the DataFrame<\/code> containing the subset of valid data. The names in the formula string correspond to columns in the DataFrame<\/code>.<\/p>\n\n\n\n
The result from ols<\/code> is an object that represents the model \u2013 it provides a function called fit<\/code> that does the actual computation.<\/p>\n\n\n\n
The result is a RegressionResultsWrapper<\/code>, which contains a Series<\/code> called params<\/code>, which contains the estimated intercept and the slope associated with educ<\/code>.<\/p>\n\n\n\n
results.params\n<\/pre>\n\n\n\n
Intercept   -15007.453641\neduc          3631.076100\ndtype: float64\n<\/pre>\n\n\n\n
The results from Statsmodels are the same as the results we got from SciPy, so that\u2019s good!<\/p>\n\n\n\n
Multiple Regression<\/h2>\n\n\n\n
In the previous section, we saw that income depends on education, and in the exercise we saw that it also depends on age<\/code>. Now let\u2019s put them together in a single model.<\/p>\n\n\n\n
results = smf.ols('realinc ~ educ + age', data=gss).fit()\nresults.params\n<\/pre>\n\n\n\n
Intercept   -17999.726908\neduc          3665.108238\nage             55.071802\ndtype: float64\n<\/pre>\n\n\n\n
In this model, realinc<\/code> is the variable we are trying to explain or predict, which is called the dependent variable<\/strong> because it depends on the the other variables \u2013 or at least we expect it to. The other variables, educ<\/code> and age<\/code>, are called independent variables<\/strong> or sometimes \u201cpredictors\u201d. The +<\/code> sign indicates that we expect the contributions of the independent variables to be additive.<\/p>\n\n\n\n
The result contains an intercept and two slopes, which estimate the average contribution of each predictor with the other predictor held constant.<\/p>\n\n\n\n