Week05_Statsmodels

July 17, 2024

1 Statistical analysis in Statsmodels

1.1 Overview
Package Statsmodels is a powerful Python library for performing statistical data analysis –
[Link] It provides classes and functions for estimating
and testing different statistical models, performing hypothesis tests, and data exploration.

1.2 Importing Statsmodels

As usual, to use the Statsmodels library, you need to import it along with other necessary libraries.
Important note: there are two different ways to define models in Statsmodels, so I import two
versions of the main API at the same time. You will see below how it works. On practice, you will
stick to one method and you will use one import only.

[1]: import [Link] as sm

import [Link] as smf
import numpy as np
import pandas as pd
import [Link] as plt
import seaborn as sns

2 Overview of Statsmodels Modules

Statsmodels provides a variety of modules for different statistical analyses. Saying “analysis”,
I mean really advanced statistical analysis – way more complex than you have seen before. In
most cases, these are models for multivariate statistical analysis – analysis of the data with many
variables at the same time. This is opposite to the univariate analysis you had in the course
Statistics for Data Science. Here are some of the included modules: - [Link]: The main
API for accessing various functionalities. - [Link]: Provides a formula-based
interface similar to R. - [Link].linear_model: Linear regression models. -
[Link].generalized_linear_model: Generalized linear models, including logistic
regression. - [Link]: Time series analysis. - [Link]: Experimental and
less-stable code.
For this presentation, we will focus on simple linear model (OLS regression) and logistic regression.
Something you already know really well.

1
 3 Simple OLS Linear Model
3.1 Defining the Model by Data
3.1.1 Example Data
Let’s create a simple synthetic dataset for demonstration:

[2]: # Generating example data

[Link](0)
X = [Link](100)
y = 2 * X + [Link](0, 0.2, 100)

# Print out some data

print('first 10 values of X:', X[:10])
print('first 10 values of y:', y[:10])

# Plot the data

[Link](X,y);

first 10 values of X: [0.5488135 0.71518937 0.60276338 0.54488318 0.4236548

0.64589411
0.43758721 0.891773 0.96366276 0.38344152]
first 10 values of y: [0.86459704 1.61054403 1.29865924 0.78251763 1.14496004
1.67096606
1.11093034 1.74756103 1.713175 0.97777338]

2
 3.1.2 Defining and Fitting the Model
There are two ways to define the model: by data and by formula. For the first method, we need
to use module sm as we defined it before and prepare the data properly. For the linear model,
it means we need to add a constant for the intercept. You should remember that a simple linear
model with one input variable has two parameters: slope and intercept. The slope is defined by
the independent variable X, intercept should be added manually as a number 1.

[3]: # Adding a constant for the intercept

X = sm.add_constant(X)

# Print out X data

print('first 10 values of X:\n', X[:10])

first 10 values of X:
[[1. 0.5488135 ]
[1. 0.71518937]
[1. 0.60276338]
[1. 0.54488318]
[1. 0.4236548 ]
[1. 0.64589411]
[1. 0.43758721]
[1. 0.891773 ]
[1. 0.96366276]
[1. 0.38344152]]
Now we are ready to create and fit the model. Check the help file for [Link](): the first parameter is
the dependent variable or target, and the second parameter is the independent variable or predictor.
There might be multiple predictors but we focus here on a simple linear model with one predictor
only.

[4]: # Define the model

model = [Link](y, X)

# Fit the model

results = [Link]()

# Print the summary

print([Link]())

OLS Regression Results

==============================================================================
Dep. Variable: y R-squared: 0.892
Model: OLS Adj. R-squared: 0.891
Method: Least Squares F-statistic: 810.4
Date: Wed, 17 Jul 2024 Prob (F-statistic): 3.49e-49

3
 Time: 12:39:18 Log-Likelihood: 19.429
No. Observations: 100 AIC: -34.86
Df Residuals: 98 BIC: -29.65
Df Model: 1
Covariance Type: nonrobust
==============================================================================
coef std err t P>|t| [0.025 0.975]
------------------------------------------------------------------------------
const 0.0444 0.039 1.149 0.253 -0.032 0.121
x1 1.9874 0.070 28.468 0.000 1.849 2.126
==============================================================================
Omnibus: 11.746 Durbin-Watson: 2.083
Prob(Omnibus): 0.003 Jarque-Bera (JB): 4.097
Skew: 0.138 Prob(JB): 0.129
Kurtosis: 2.047 Cond. No. 4.30
==============================================================================

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

3.2 Defining the Model by Formula

3.2.1 Example Data
The second method defines the model by formula. This is a very powerful method, and you will
see it again in R. Again, we need to prepare the data correctly. For the formula method, the data
should be Pandas DataFrame with predictor and target as columns. You don’t need to worry about
a constant for the intercept – it will be done automatically.

[5]: # Make data generated before in DataFrame format

df = [Link]({'X': X[:,1], 'y': y})

# Print some data

print([Link]())

# Plot the data

[Link](df, x = "X", y = "y");

X y
0 0.548814 0.864597
1 0.715189 1.610544
2 0.602763 1.298659
3 0.544883 0.782518
4 0.423655 1.144960

4
 3.2.2 Defining and Fitting the Model
String parameter 'y ~ X' in the function [Link]() is the formula. It references column names
in the data frame df and tells that column y depends on column X.
If there will be more predictors, for example, there would be independent variables (predictors) X,
W, Z in the data frame df, then the formula could be 'y ~ X + W + Z' – y depends on the linear
combination of X and W and Z, each predictor would have its own slope coefficient. We will talk
more about it when we get to R, formulas in R are more powerful and flexible than in Python.
Now, get back to our simple linear model.

[6]: # Define the model using a formula

model = [Link]('y ~ X', data=df)

# Fit the model

results = [Link]()

# Print the summary

print([Link]())

OLS Regression Results

5
==============================================================================
Dep. Variable: y R-squared: 0.892
Model: OLS Adj. R-squared: 0.891
Method: Least Squares F-statistic: 810.4
Date: Wed, 17 Jul 2024 Prob (F-statistic): 3.49e-49
Time: 12:39:18 Log-Likelihood: 19.429
No. Observations: 100 AIC: -34.86
Df Residuals: 98 BIC: -29.65
Df Model: 1
Covariance Type: nonrobust
==============================================================================
coef std err t P>|t| [0.025 0.975]
------------------------------------------------------------------------------
Intercept 0.0444 0.039 1.149 0.253 -0.032 0.121
X 1.9874 0.070 28.468 0.000 1.849 2.126
==============================================================================
Omnibus: 11.746 Durbin-Watson: 2.083
Prob(Omnibus): 0.003 Jarque-Bera (JB): 4.097
Skew: 0.138 Prob(JB): 0.129
Kurtosis: 2.047 Cond. No. 4.30
==============================================================================

Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly
specified.
Obviously, the result is the same as before as we used the same data for X and y. You got a lot of
numbers. What do they mean? What numbers are the most important?
You got several p-values in this report. Remember, every p-value means that there is a hypothesis
testing behind it with its own null- and alternative hypothesis. What are they here? What is your
conclusion?
In real-life data analysis, we don’t formally define our hypothesis testing. We skip that part for the
report. But we should have a really good understanding of all involved hypothesis testings (plural)
to be able to make the right conclusions about the model.
Finally, you need to provide an interpretation of the model. What is it? If there is no interpretation
in your report, it means you have no model at all. Data analytics is not about programming code
– it is about correct interpretations for your results. Obviously, the results should be correct too.

3.3 Diagnostic data visualisations

After fitting a linear model using Statsmodels, it’s important to perform diagnostic checks to ensure
that the assumptions of the linear regression are met. These assumptions include linearity, inde-
pendence, homoscedasticity, and normality of residuals. Statsmodels provides several diagnostic
plots to help assess these assumptions.
For any model, the variable results we got as a result of fitting the model is a very complex custom
class object. As with any other object in Python, we can use the command dir(results) to see

6
 available functionality.

[7]: print(dir(results))

['HC0_se', 'HC1_se', 'HC2_se', 'HC3_se', '_HCCM', '__class__', '__delattr__',

'__dict__', '__dir__', '__doc__', '__eq__', '__format__', '__ge__',
'__getattribute__', '__getstate__', '__gt__', '__hash__', '__init__',
'__init_subclass__', '__le__', '__lt__', '__module__', '__ne__', '__new__',
'__reduce__', '__reduce_ex__', '__repr__', '__setattr__', '__sizeof__',
'__str__', '__subclasshook__', '__weakref__', '_abat_diagonal', '_cache',
'_data_attr', '_data_in_cache', '_get_robustcov_results', '_get_wald_nonlinear',
'_is_nested', '_transform_predict_exog', '_use_t', '_wexog_singular_values',
'aic', 'bic', 'bse', 'centered_tss', 'compare_f_test', 'compare_lm_test',
'compare_lr_test', 'condition_number', 'conf_int', 'conf_int_el', 'cov_HC0',
'cov_HC1', 'cov_HC2', 'cov_HC3', 'cov_kwds', 'cov_params', 'cov_type',
'df_model', 'df_resid', 'diagn', 'eigenvals', 'el_test', 'ess', 'f_pvalue',
'f_test', 'fittedvalues', 'fvalue', 'get_influence', 'get_prediction',
'get_robustcov_results', 'info_criteria', 'initialize', 'k_constant', 'llf',
'load', 'model', 'mse_model', 'mse_resid', 'mse_total', 'nobs',
'normalized_cov_params', 'outlier_test', 'params', 'predict', 'pvalues',
'remove_data', 'resid', 'resid_pearson', 'rsquared', 'rsquared_adj', 'save',
'scale', 'ssr', 'summary', 'summary2', 't_test', 't_test_pairwise', 'tvalues',
'uncentered_tss', 'use_t', 'wald_test', 'wald_test_terms', 'wresid']
We can use information from the variable results to make diagnostic plots

[8]: # Residual Plot

[Link](figsize=(8, 4))
[Link]([Link], [Link])
[Link](y=0, color='r', linestyle='--')
[Link]('Fitted values')
[Link]('Residuals')
[Link]('Residual Plot')
[Link]()

7
 Residual Plot checks for homoscedasticity and linearity. Points should be randomly scattered
around the horizontal axis (residuals = 0). Random scatter indicates homoscedasticity (constant
variance). Patterns indicate non-linearity or heteroscedasticity (non-constant variance).

[9]: # Q-Q Plot

[Link](figsize=(8, 4))
[Link]([Link], line='s')
[Link]('Q-Q Plot')
[Link]();

<Figure size 800x400 with 0 Axes>

8
 Q-Q Plot checks if residuals are normally distributed. Points should lie on or close to the reference
line. Points on the line indicate normally distributed residuals. Deviations from the line indicate
departures from normality. As an addition to this graph, you can plot a histogram of residuals
[Link].

[10]: # Scale-Location Plot

[Link](figsize=(8, 4))
[Link]([Link], [Link]([Link]([Link])))
[Link]('Fitted values')
[Link]('Sqrt(|Residuals|)')
[Link]('Scale-Location Plot')
[Link]()

9
 Scale-Location Plot checks for homoscedasticity. It looks and works similar to the Residual Plot
above, but it looks on the square root of residuals. Interpretation is similar too. Points should be
randomly scattered horizontally. Random scatter indicates homoscedasticity. Patterns or trends
indicate heteroscedasticity.

[11]: # Leverage Plot

[Link](figsize=(8, 4))
[Link].influence_plot(results, criterion="cooks")
[Link]('Leverage Plot')
[Link]()

<Figure size 800x400 with 0 Axes>

10
Leverage Plot identifies influential data points – potential outliers – with high leverage and/or
high residuals. Points outside the “Cook’s distance” lines are highly influential. Points with high
leverage may disproportionately affect the model. I did not get lines for Cook’s distance on this
graph as our example data was “perfect” – there were no outliers. You will see them on real data.

4 Logistic Regression
Now, let’s have a look at a different model. The linear model above made prediction about the
numerical variable – target y was numerical. Logistic regression makes prediction for a categorical
variable – target y should be binary, like 0 and 1 or True and False.

4.1 Defining the Model by Data

4.1.1 Example Data
Again, we start with simple synthetic data. Let’s create a simple dataset for logistic regression.
Now we will have two predictors, two independent variables X1 and X2.

11
[12]: # Generating example data
[Link](0)
X1 = [Link](100)
X2 = [Link](100)
y = (X1 + X2 + [Link](0, 0.2, 100) > 1).astype(int)

# Combine two independent variables in one object

X = np.column_stack((X1, X2))

# Print out some data

print('first 10 values of X:\n', X[:10])
print('first 10 values of y:\n', y[:10])

# Plot the data

# Code c = (y + 1) is a "fancy" way to set two different colours
[Link](X1, X2, c = (y + 1));

first 10 values of X:
[[0.5488135 0.67781654]
[0.71518937 0.27000797]
[0.60276338 0.73519402]
[0.54488318 0.96218855]
[0.4236548 0.24875314]
[0.64589411 0.57615733]
[0.43758721 0.59204193]
[0.891773 0.57225191]
[0.96366276 0.22308163]
[0.38344152 0.95274901]]
first 10 values of y:
[1 0 1 1 0 1 1 1 0 1]

12
 4.1.2 Defining and Fitting the Model
Again we need to add a constant and then we can fit the model and see results.

[13]: # Adding a constant for the intercept

X = sm.add_constant(X)

# Define the model

model = [Link](y, X)

# Fit the model

results = [Link]()

# Print the summary

print([Link]())

Optimization terminated successfully.

Current function value: 0.307807
Iterations 8
Logit Regression Results
==============================================================================
Dep. Variable: y No. Observations: 100

13
 Model: Logit Df Residuals: 97
Method: MLE Df Model: 2
Date: Wed, 17 Jul 2024 Pseudo R-squ.: 0.5527
Time: 12:39:19 Log-Likelihood: -30.781
converged: True LL-Null: -68.814
Covariance Type: nonrobust LLR p-value: 3.037e-17
==============================================================================
coef std err z P>|z| [0.025 0.975]
------------------------------------------------------------------------------
const -8.2509 1.719 -4.800 0.000 -11.620 -4.882
x1 7.7653 1.674 4.639 0.000 4.484 11.046
x2 9.0683 1.902 4.767 0.000 5.340 12.797
==============================================================================

4.2 Defining the Model by Formula

4.2.1 Example Data
To define the model by formula we need to have a data frame.

[14]: # Generating example data in DataFrame format

df = [Link]({
'X1': X1,
'X2': X2,
'y': y
})

# Print some data

print([Link]())

# Plot the data

[Link](df, x = "X1", y = "X2", hue = "y");

X1 X2 y
0 0.548814 0.677817 1
1 0.715189 0.270008 0
2 0.602763 0.735194 1
3 0.544883 0.962189 1
4 0.423655 0.248753 0

14
 4.2.2 Defining and Fitting the Model

[15]: # Define and fit the model using a formula

model = [Link]('y ~ X1 + X2', data=df)
results = [Link]()

# Print the summary

print([Link]())

Optimization terminated successfully.

Current function value: 0.307807
Iterations 8
Logit Regression Results
==============================================================================
Dep. Variable: y No. Observations: 100
Model: Logit Df Residuals: 97
Method: MLE Df Model: 2
Date: Wed, 17 Jul 2024 Pseudo R-squ.: 0.5527
Time: 12:39:19 Log-Likelihood: -30.781
converged: True LL-Null: -68.814

15
 Covariance Type: nonrobust LLR p-value: 3.037e-17
==============================================================================
coef std err z P>|z| [0.025 0.975]
------------------------------------------------------------------------------
Intercept -8.2509 1.719 -4.800 0.000 -11.620 -4.882
X1 7.7653 1.674 4.639 0.000 4.484 11.046
X2 9.0683 1.902 4.767 0.000 5.340 12.797
==============================================================================
The results from both approaches are the same. Now you need to provide interpretations for these
results.
Similar to the linear model, the variable results from logistic regression analysis is a custom class
object – however different to the results object from the linear model. Again, we can use the
command dir(results) to see available functionality.

[16]: print(dir(results))

['__class__', '__delattr__', '__dict__', '__dir__', '__doc__', '__eq__',

'__format__', '__ge__', '__getattribute__', '__getstate__', '__gt__',
'__hash__', '__init__', '__init_subclass__', '__le__', '__lt__', '__module__',
'__ne__', '__new__', '__reduce__', '__reduce_ex__', '__repr__', '__setattr__',
'__sizeof__', '__str__', '__subclasshook__', '__weakref__', '_cache',
'_data_attr', '_data_in_cache', '_get_endog_name', '_get_robustcov_results',
'_get_wald_nonlinear', '_transform_predict_exog', '_use_t', 'aic', 'bic', 'bse',
'conf_int', 'converged', 'cov_kwds', 'cov_params', 'cov_type', 'df_model',
'df_resid', 'f_test', 'fittedvalues', 'get_distribution', 'get_influence',
'get_margeff', 'get_prediction', 'im_ratio', 'info_criteria', 'initialize',
'k_constant', 'llf', 'llnull', 'llr', 'llr_pvalue', 'load', 'method',
'mle_retvals', 'mle_settings', 'model', 'nobs', 'normalized_cov_params',
'params', 'pred_table', 'predict', 'prsquared', 'pvalues', 'remove_data',
'resid_dev', 'resid_generalized', 'resid_pearson', 'resid_response', 'save',
'scale', 'score_test', 'set_null_options', 'summary', 'summary2', 't_test',
't_test_pairwise', 'tvalues', 'use_t', 'wald_test', 'wald_test_terms']
And then we can use this information for diagnostic and interpretation. For example,

[17]: # Parameters of the model

print("Model parameters:\n", [Link], "\n", sep = "")

# They are very useful for better interpretations of the model as for the␣
↪logistic

# regression we prefer to discuss odds ratio rather than raw coefficients

print("Odds ratio:\n", [Link]([Link]), sep = "")

Model parameters:
Intercept -8.250864
X1 7.765251
X2 9.068314
dtype: float64

16
 Odds ratio:
Intercept 0.000261
X1 2357.249566
X2 8675.979560
dtype: float64

[18]: # Extract predicted values and compare them to true values of the target

# Make predictions
pred = ([Link]() > 0.5).astype(int)

# Compute the confusion matrix using [Link]

conf_matrix = [Link](y, pred, rownames=['Actual'], colnames=['Predicted'],␣
↪margins=False)

print("Confusion matrix:\n", conf_matrix, "\n", sep = "")

# Calculate accuracy
print("Accuracy: ", [Link](conf_matrix).sum() / conf_matrix.sum().sum())

Confusion matrix:
Predicted 0 1
Actual
0 37 8
1 6 49

Accuracy: 0.86
There are other tools available to calculate confusion matrix, accuracy and many different things.
You don’t need to do that “manually”. You will see more advanced tools next week.

5 Datasets in Statsmodels
There are a number of datasets provided in the Statsmodels package –
[Link] You
can use them for testing and exercises. For example, load Longley dataset –
[Link] – and try to see if there are
any relationships between population, employment and size of armed forces.

[19]: # Load the dataset object

data = [Link].load_pandas()

# Extract the actual data from the dataset object

df = [Link]

# Have a look on the top of the data frame

[Link]()

17
[19]: TOTEMP GNPDEFL GNP UNEMP ARMED POP YEAR
0 60323.0 83.0 234289.0 2356.0 1590.0 107608.0 1947.0
1 61122.0 88.5 259426.0 2325.0 1456.0 108632.0 1948.0
2 60171.0 88.2 258054.0 3682.0 1616.0 109773.0 1949.0
3 61187.0 89.5 284599.0 3351.0 1650.0 110929.0 1950.0
4 63221.0 96.2 328975.0 2099.0 3099.0 112075.0 1951.0

[ ]:

18