Linear Regression Diagnostic in Python with StatsModels
import numpy as np import statsmodels import seaborn as sns from matplotlib import pyplot as plt %matplotlib inline
While linear regression is a pretty simple task, there are several assumptions for the model that we may want to validate. I follow the regression diagnostic here, trying to justify four principal assumptions, namely LINE in Python:
- Independence (This is probably more serious for time series. I’ll pass it for now)
- Equal variance (or homoscedasticity)
I learnt this abbreviation of linear regression assumptions when I was taking a course on correlation and regression taught by Walter Vispoel at UIowa. Really helped me to remember these four little things!
statsmodels itself contains useful modules for regression diagnostics. In addition to those, I want to go with somewhat manual yet very simple ways for more flexible visualizations.
Let’s go with the depression data. More toy datasets can be found here. For simplicity, I randomly picked 3 columns.
> df = statsmodels.datasets.get_rdataset("Ginzberg", "car").data > df = df[['adjdep', 'adjfatal', 'adjsimp']] > df.head(2)
Linear regression is simple, with
statsmodels. We are able to use R style regression formula.
> import statsmodels.formula.api as smf > reg = smf.ols('adjdep ~ adjfatal + adjsimp', data=df).fit() > reg.summary()
|Date:||Tue, 01 May 2018||Prob (F-statistic):||1.82e-10|
Now let’s try to validate the four assumptions one by one
Linearity & Equal variance
Both can be tested by plotting residuals vs. predictions, where residuals are prediction errors.
> pred_val = reg.fittedvalues.copy() > true_val = df['adjdep'].values.copy() > residual = true_val - pred_val
> fig, ax = plt.subplots(figsize=(6,2.5)) > _ = ax.scatter(residual, pred_val)
It seems like the corresponding residual plot is reasonably random. To confirm that, let’s go with a hypothesis test, Harvey-Collier multiplier test, for linearity
> import statsmodels.stats.api as sms > sms.linear_harvey_collier(reg) Ttest_1sampResult(statistic=4.990214882983107, pvalue=3.5816973971922974e-06)
Several tests exist for equal variance, with different alternative hypotheses. Let’s go with Breusch-Pagan test as an example. More can be found here. Small p-value (
pval below) shows that there is violation of homoscedasticity.
> _, pval, __, f_pval = statsmodels.stats.diagnostic.het_breuschpagan(residual, df[['adjfatal', 'adjsimp']]) > pval, f_pval (6.448482473013975e-08, 2.21307383960346e-08)
Usually assumption violations are not independent of each other. Having one violations may lead to another. In this case, we see that both linearity and homoscedasticity are not met. Possible data transformation such as log, Box-Cox power transformation, and other fixes may be needed to get a better regression outcome.
We can apply normal probability plot to assess how the data (error) depart from normality visually:
> import scipy as sp > fig, ax = plt.subplots(figsize=(6,2.5)) > _, (__, ___, r) = sp.stats.probplot(residual, plot=ax, fit=True) > r**2 0.9523990893322951
The good fit indicates that normality is a reasonable approximation.