3.1.6.4. Simple Regression

Fit a simple linear regression using ‘statsmodels’, compute corresponding p-values.

# Original author: Thomas Haslwanter

import numpy as np
import matplotlib.pyplot as plt
import pandas

# For statistics. Requires statsmodels 5.0 or more
from statsmodels.formula.api import ols

# Analysis of Variance (ANOVA) on linear models
from statsmodels.stats.anova import anova_lm

Generate and show the data

x = np.linspace(-5, 5, 20)

# To get reproducible values, provide a seed value
rng = np.random.default_rng(27446968)

y = -5 + 3 * x + 4 * np.random.normal(size=x.shape)

# Plot the data
plt.figure(figsize=(5, 4))
plt.plot(x, y, "o")

[<matplotlib.lines.Line2D object at 0x7f9519d9af90>]

Multilinear regression model, calculating fit, P-values, confidence intervals etc.

# Convert the data into a Pandas DataFrame to use the formulas framework
# in statsmodels
data = pandas.DataFrame({"x": x, "y": y})

# Fit the model
model = ols("y ~ x", data).fit()

# Print the summary
print(model.summary())

# Perform analysis of variance on fitted linear model
anova_results = anova_lm(model)

print("\nANOVA results")
print(anova_results)

                            OLS Regression Results
==============================================================================
Dep. Variable:                      y   R-squared:                       0.875
Model:                            OLS   Adj. R-squared:                  0.868
Method:                 Least Squares   F-statistic:                     125.7
Date:                Fri, 15 Sep 2023   Prob (F-statistic):           1.50e-09
Time:                        19:08:01   Log-Likelihood:                -55.922
No. Observations:                  20   AIC:                             115.8
Df Residuals:                      18   BIC:                             117.8
Df Model:                           1
Covariance Type:            nonrobust
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
Intercept     -5.7794      0.934     -6.186      0.000      -7.742      -3.817
x              3.4512      0.308     11.211      0.000       2.804       4.098
==============================================================================
Omnibus:                        0.164   Durbin-Watson:                   3.251
Prob(Omnibus):                  0.921   Jarque-Bera (JB):                0.110
Skew:                           0.130   Prob(JB):                        0.946
Kurtosis:                       2.747   Cond. No.                         3.03
==============================================================================

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

ANOVA results
            df       sum_sq      mean_sq           F        PR(>F)
x          1.0  2194.066963  2194.066963  125.693214  1.495554e-09
Residual  18.0   314.203162    17.455731         NaN           NaN

Plot the fitted model

# Retrieve the parameter estimates
offset, coef = model._results.params
plt.plot(x, x * coef + offset)
plt.xlabel("x")
plt.ylabel("y")

plt.show()