statsmodelsLet’s explore linear regression using a familiar example dataset of student grades. Our goal will be to train a model to predict a student’s grade given the number of hours they have studied.
In this implementation, we will use the statsmodels package to achieve this.
Loading the data:
from pandas import read_csv
repo_url = "https://raw.githubusercontent.com/prof-rossetti/python-for-finance"
request_url = f"{repo_url}/main/docs/data/grades.csv"
df = read_csv(request_url)
df.head()| Name | StudyHours | Grade | |
|---|---|---|---|
| 0 | Arun | 10.00 | 50.0 |
| 1 | Sofia | 11.50 | 50.0 |
| 2 | Hassan | 9.00 | 47.0 |
| 3 | Zara | 16.00 | 97.0 |
| 4 | Liam | 9.25 | 49.0 |
Dropping null values:
df.dropna(inplace=True)
df.tail()| Name | StudyHours | Grade | |
|---|---|---|---|
| 17 | Tariq | 6.0 | 35.0 |
| 18 | Lakshmi | 10.0 | 48.0 |
| 19 | Maya | 12.0 | 52.0 |
| 20 | Yusuf | 12.5 | 63.0 |
| 21 | Zainab | 12.0 | 64.0 |
Exploring relationship between variables:
import plotly.express as px
px.scatter(df, x="StudyHours", y="Grade", height=350,
title="Relationship between Study Hours and Grades",
trendline="ols", trendline_color_override="red",
)Identifying the dependent and independent variables:
x = df["StudyHours"]
print(x.shape)
y = df["Grade"]
print(y.shape)(22,)
(22,)When using statsmodels, the documentation instructs us to manually add a column of ones (to help the model perform calculations related to the y-intercept):
import statsmodels.api as sm
x = sm.add_constant(x) # adding in a column of constants, as per the OLS docs
x.head()| const | StudyHours | |
|---|---|---|
| 0 | 1.0 | 10.00 |
| 1 | 1.0 | 11.50 |
| 2 | 1.0 | 9.00 |
| 3 | 1.0 | 16.00 |
| 4 | 1.0 | 9.25 |
Now we split the training and test sets:
from sklearn.model_selection import train_test_split
x_train, x_test, y_train, y_test = train_test_split(x, y, random_state=99)
print("TRAIN:", x_train.shape, y_train.shape)
print("TEST:", x_test.shape, y_test.shape)TRAIN: (16, 2) (16,)
TEST: (6, 2) (6,)Selecting a linear regression (OLS) model, and training it on the training data to learn the ideal weights:
import statsmodels.api as sm
model = sm.OLS(y_train, x_train, missing="drop")
print(type(model))
results = model.fit()
print(type(results))<class 'statsmodels.regression.linear_model.OLS'>
<class 'statsmodels.regression.linear_model.RegressionResultsWrapper'>It is possible to access the training results, including some summary statistics:
print(results.summary()) OLS Regression Results
==============================================================================
Dep. Variable: Grade R-squared: 0.865
Model: OLS Adj. R-squared: 0.855
Method: Least Squares F-statistic: 89.55
Date: Tue, 19 Nov 2024 Prob (F-statistic): 1.84e-07
Time: 20:54:26 Log-Likelihood: -56.257
No. Observations: 16 AIC: 116.5
Df Residuals: 14 BIC: 118.1
Df Model: 1
Covariance Type: nonrobust
==============================================================================
coef std err t P>|t| [0.025 0.975]
------------------------------------------------------------------------------
const -17.9242 7.078 -2.532 0.024 -33.105 -2.743
StudyHours 6.3637 0.672 9.463 0.000 4.921 7.806
==============================================================================
Omnibus: 0.338 Durbin-Watson: 2.313
Prob(Omnibus): 0.845 Jarque-Bera (JB): 0.089
Skew: -0.162 Prob(JB): 0.957
Kurtosis: 2.833 Cond. No. 34.5
==============================================================================
Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified./opt/hostedtoolcache/Python/3.11.10/x64/lib/python3.11/site-packages/scipy/stats/_axis_nan_policy.py:418: UserWarning:
`kurtosistest` p-value may be inaccurate with fewer than 20 observations; only n=16 observations were given.
In a regression analysis, each coefficient (the number associated with a feature in the model) has a corresponding t-statistic that tests whether the coefficient is meaningfully different from zero.
Interpreting the results:
T-statistic: Measures how many standard deviations the coefficient is away from zero. A larger t-statistic suggests that the coefficient is far from zero and potentially significant.
P-value (P>|t|): This tells you the probability that the observed t-statistic would occur if the coefficient were actually zero (the null hypothesis). If this probability is very small (typically < 0.05), it means it’s unlikely that the coefficient is zero, suggesting it is statistically significant.
Interpreting p-values:
The training results contain an r-squared score, however this represents the error for the training data. To get the real results of how the model generalizes to the test data, we will calculate the r-squared score and other metrics on the test results later.
The part of the training results we care about are the the learned weights (i.e. coefficients), which we use to arrive at the line of best fit:
print(results.params)
print("-----------")
print(f"y =", f"{results.params['StudyHours'].round(3)}x ",
f"+ {results.params['const'].round(3)}"
)const -17.924187
StudyHours 6.363725
dtype: float64
-----------
y = 6.364x + -17.924The training results also contain the fittedvalues (predictions), as well as the resid (residuals or errors). We can compare each of the predicted values against the actual known values, to verify the residuals for illustration purposes:
from pandas import DataFrame
# get all rows from the original dataset that wound up in the training set:
training_set = df.loc[x_train.index].copy()
# create a dataset for the predictions and the residuals:
training_preds = DataFrame({
"Predictions": results.fittedvalues,
"Residuals": results.resid
})
# merge the training set with the results:
training_set = training_set.merge(training_preds, how="inner",
left_index=True, right_index=True
)
# calculate error for each datapoint:
training_set["My Error"] = training_set["Grade"] - training_set["Predictions"]
training_set| Name | StudyHours | Grade | Predictions | Residuals | My Error | |
|---|---|---|---|---|---|---|
| 15 | Anika | 8.00 | 27.0 | 32.985616 | -5.985616 | -5.985616 |
| 0 | Arun | 10.00 | 50.0 | 45.713067 | 4.286933 | 4.286933 |
| 18 | Lakshmi | 10.00 | 48.0 | 45.713067 | 2.286933 | 2.286933 |
| 13 | Mei | 8.00 | 15.0 | 32.985616 | -17.985616 | -17.985616 |
| 12 | Priya | 9.00 | 37.0 | 39.349341 | -2.349341 | -2.349341 |
| 20 | Yusuf | 12.50 | 63.0 | 61.622380 | 1.377620 | 1.377620 |
| 7 | Kwame | 9.00 | 42.0 | 39.349341 | 2.650659 | 2.650659 |
| 21 | Zainab | 12.00 | 64.0 | 58.440518 | 5.559482 | 5.559482 |
| 16 | Elif | 9.00 | 36.0 | 39.349341 | -3.349341 | -3.349341 |
| 5 | Xia | 1.00 | 3.0 | -11.560462 | 14.560462 | 14.560462 |
| 4 | Liam | 9.25 | 49.0 | 40.940273 | 8.059727 | 8.059727 |
| 19 | Maya | 12.00 | 52.0 | 58.440518 | -6.440518 | -6.440518 |
| 9 | Takumi | 14.50 | 74.0 | 74.349831 | -0.349831 | -0.349831 |
| 8 | Fatima | 8.50 | 26.0 | 36.167479 | -10.167479 | -10.167479 |
| 3 | Zara | 16.00 | 97.0 | 83.895419 | 13.104581 | 13.104581 |
| 1 | Sofia | 11.50 | 50.0 | 55.258655 | -5.258655 | -5.258655 |
It is possible to calculate the training metrics ourselves, to verify the regression results summary we saw above:
from sklearn.metrics import r2_score, mean_squared_error, mean_absolute_error
r2 = r2_score(training_set["Grade"], training_set["Predictions"])
print("R^2:", round(r2, 3))
mae = mean_absolute_error(training_set["Grade"], training_set["Predictions"])
print("MAE:", round(mae, 3))
mse = mean_squared_error(training_set["Grade"], training_set["Predictions"])
print("MSE:", round(mse,3))
rmse = mse ** .5
print("RMSE:", rmse.round(3))R^2: 0.865
MAE: 6.486
MSE: 66.301
RMSE: 8.143Remember, these metrics only tell us about the predictive accuracy on the training data, and we care more about evaluating metrics on the test set.
Alright, we trained the model, but how well does it do in making predictions?
We use the trained model to make predictions on the unseen (test) data:
preds = results.get_prediction(x_test)
print(type(preds))
preds_df = preds.summary_frame(alpha=0.05)
preds_df<class 'statsmodels.regression._prediction.PredictionResults'>| mean | mean_se | mean_ci_lower | mean_ci_upper | obs_ci_lower | obs_ci_upper | |
|---|---|---|---|---|---|---|
| 17 | 20.258165 | 3.468123 | 12.819780 | 27.696550 | 0.161099 | 40.355231 |
| 14 | 80.713556 | 4.282255 | 71.529034 | 89.898079 | 59.906876 | 101.520237 |
| 2 | 39.349341 | 2.280844 | 34.457418 | 44.241265 | 20.049253 | 58.649430 |
| 10 | 80.713556 | 4.282255 | 71.529034 | 89.898079 | 59.906876 | 101.520237 |
| 6 | 55.258655 | 2.394199 | 50.123609 | 60.393701 | 35.895515 | 74.621795 |
| 11 | 69.577037 | 3.322979 | 62.449955 | 76.704119 | 49.593099 | 89.560975 |
When we use the summary_frame method on prediction results, it returns a DataFrame with several columns. Here’s a breakdown of what they mean:
Prediction (mean): This is the predicted value for each observation based on the model. It’s essentially the point prediction (ŷ) for the corresponding input.
Standard Error (mean_se): This stands for the standard error of the predicted mean. It measures the uncertainty associated with the predicted value due to sampling variability. A smaller mean_se indicates higher confidence in the predicted mean.
Confidence Interval (mean_ci_lower and mean_ci_upper): Represents the range in which the true mean prediction is likely to lie. For predicting the average value (e.g. “the average apple weight is between 140 and 160 grams”).
Prediction Interval (obs_ci_lower and obs_ci_upper): Represents the range in which an individual new observation is likely to lie. For predicting the range where individual values could fall (e.g. “an individual apple might weigh between 120 and 180 grams”).
Merging the actual values in, so we can compare predicted values vs actual values:
# reconstruct test set:
test_set = df.loc[x_test.index].copy()
# merge in prediction results:
test_set = test_set.merge(preds_df, how="inner",
left_index=True, right_index=True
)
test_set.rename(columns={"mean":"Prediction", "mean_se": "Standard Error",
"mean_ci_upper": "CI Upper", "mean_ci_lower": "CI Lower",
"obs_ci_upper": "PI Upper", "obs_ci_lower": "PI Lower",
}, inplace=True
)
test_set["My Error"] = test_set["Grade"] - test_set["Prediction"]
test_set| Name | StudyHours | Grade | Prediction | Standard Error | CI Lower | CI Upper | PI Lower | PI Upper | My Error | |
|---|---|---|---|---|---|---|---|---|---|---|
| 17 | Tariq | 6.00 | 35.0 | 20.258165 | 3.468123 | 12.819780 | 27.696550 | 0.161099 | 40.355231 | 14.741835 |
| 14 | Alberto | 15.50 | 70.0 | 80.713556 | 4.282255 | 71.529034 | 89.898079 | 59.906876 | 101.520237 | -10.713556 |
| 2 | Hassan | 9.00 | 47.0 | 39.349341 | 2.280844 | 34.457418 | 44.241265 | 20.049253 | 58.649430 | 7.650659 |
| 10 | Leila | 15.50 | 82.0 | 80.713556 | 4.282255 | 71.529034 | 89.898079 | 59.906876 | 101.520237 | 1.286444 |
| 6 | Carlos | 11.50 | 53.0 | 55.258655 | 2.394199 | 50.123609 | 60.393701 | 35.895515 | 74.621795 | -2.258655 |
| 11 | Jin | 13.75 | 62.0 | 69.577037 | 3.322979 | 62.449955 | 76.704119 | 49.593099 | 89.560975 | -7.577037 |
Visualizing the predictions:
import plotly.express as px
chart_df = test_set.copy()
# Create the plot
fig = px.scatter(chart_df, x="StudyHours", y="Grade", height=350,
title="True vs Predicted Grade (with Pred Interval and Confidence Interval)"
)
fig.add_scatter(x=chart_df["StudyHours"], y=chart_df['Prediction'],
mode='lines+markers',
name='Prediction (with PI)',
marker=dict(color='mediumpurple'), #size=10, symbol="x"
error_y=dict(type='data', symmetric=False,
array=chart_df['PI Upper'],
arrayminus=chart_df['PI Lower']),
legendrank=1)
fig.add_scatter(x=chart_df["StudyHours"], y=chart_df['Prediction'],
mode='lines+markers',
name='Prediction (with CI)',
marker=dict(color='lightblue',# size=10, #symbol="x"
),
error_y=dict(type='data', symmetric=False,
array=chart_df['CI Upper'],
arrayminus=chart_df['CI Lower']),
legendrank=1)
# Now add the scatter for the true values again, placing this plot behind
# (hack to get the actual values in front)
fig.add_scatter(x=chart_df["StudyHours"], y=chart_df['Grade'],
mode='markers',
name='True Value',
marker=dict(color='blue', size=10, symbol="x"),
legendrank=2)
fig.show()Evaluating the model using our evaluation metrics from sklearn:
from sklearn.metrics import r2_score, mean_squared_error, mean_absolute_error
y_pred = test_set["Prediction"]
r2 = r2_score(y_test, y_pred)
print("R^2:", round(r2, 3))
mae = mean_absolute_error(y_test, y_pred)
print("MAE:", round(mae, 3))
mse = mean_squared_error(y_test, y_pred)
print("MSE:", round(mse,3))
rmse = mse ** .5
print("RMSE:", rmse.round(3))R^2: 0.678
MAE: 7.371
MSE: 75.8
RMSE: 8.706We see the same results on the test set that we did with the sklearn regression implementation.