Linear Regression model with Python

Author: Matti Pastell
Date: 19.4.2013

Requirements

This en example of doing linear regression analysis using Python and statsmodels. We'll use the new formula API which makes fitting the models very familiar for R users. You'll also need Numpy, Pandas and matplolib.

The analysis can be published using Pweave 0.22 and later.

Import libraries

import pandas as pd
import numpy as np
import statsmodels.formula.api as sm
import statsmodels
import matplotlib.pyplot as plt

Statsmodels api seems to change often, check release version:

statsmodels.__version__
'0.8.0'

We'll use whiteside dataset from R package MASS. You can read the description of the dataset from the link, but in short it contains:

The weekly gas consumption and average external temperature at a house in south-east England for two heating seasons, one of 26 weeks before, and one of 30 weeks after cavity-wall insulation was installed.

Load dataset using Pandas:

url = 'https://raw.githubusercontent.com/mpastell/Rdatasets/master/csv/MASS/whiteside.csv'
whiteside = pd.read_csv(url)

Fitting the model

Let's see what the relationship between the gas consumption is before the insulation. See statsmodels documentation for more information about the syntax.

model = sm.ols(formula='Gas ~ Temp', data=whiteside, subset = whiteside['Insul']=="Before")
fitted = model.fit()
print(fitted.summary())
OLS Regression Results
==============================================================================
Dep. Variable:                    Gas   R-squared:
0.944
Model:                            OLS   Adj. R-squared:
0.941
Method:                 Least Squares   F-statistic:
403.1
Date:                Sun, 10 Sep 2017   Prob (F-statistic):
1.64e-16
Time:                        13:33:52   Log-Likelihood:
-2.8783
No. Observations:                  26   AIC:
9.757
Df Residuals:                      24   BIC:
12.27
Df Model:                           1
Covariance Type:            nonrobust
==============================================================================
                 coef    std err          t      P>|t|      [0.025
0.975]
------------------------------------------------------------------------------
Intercept      6.8538      0.118     57.876      0.000       6.609
7.098
Temp          -0.3932      0.020    -20.078      0.000      -0.434
-0.353
==============================================================================
Omnibus:                        0.296   Durbin-Watson:
2.420
Prob(Omnibus):                  0.862   Jarque-Bera (JB):
0.164
Skew:                          -0.177   Prob(JB):
0.921
Kurtosis:                       2.839   Cond. No.
13.3
==============================================================================

Warnings:
[1] Standard Errors assume that the covariance matrix of the errors is
correctly specified.

Plot the data and fit

Before = whiteside[whiteside["Insul"] == "Before"]
plt.plot(Before["Temp"], Before["Gas"], 'ro')
plt.plot(Before["Temp"], fitted.fittedvalues, 'b')
plt.legend(['Data', 'Fitted model'])
plt.ylim(0, 10)
plt.xlim(-2, 12)
plt.xlabel('Temperature')
plt.ylabel('Gas')
plt.title('Before Insulation')
<matplotlib.text.Text at 0x7fdff3a3a198>

Fit diagnostiscs

Statsmodels OLSresults objects contain the usual diagnostic information about the model and you can use the get_influence() method to get more diagnostic information (such as Cook's distance).

A look at the residuals

Histogram of normalized residuals

plt.hist(fitted.resid_pearson)
plt.ylabel('Count')
plt.xlabel('Normalized residuals')
<matplotlib.text.Text at 0x7fdff3aafa20>

Cooks distance

OLSInfluence objects contain more diagnostic information

influence = fitted.get_influence()
#c is the distance and p is p-value
(c, p) = influence.cooks_distance
plt.stem(np.arange(len(c)), c, markerfmt=",")
<Container object of 3 artists>

Statsmodels builtin plots

Statsmodels includes a some builtin function for plotting residuals against leverage:

from statsmodels.graphics.regressionplots import *
plot_leverage_resid2(fitted)
influence_plot(fitted)