Data Exploration

Linear Regression with `statsmodels`

Let’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.

Data Loading

Loading the data:

Code

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",
)

Data Splitting

X/Y Split

Identifying the dependent and independent variables:

x = df["StudyHours"]
print(x.shape)

y = df["Grade"]
print(y.shape)

(22,)
(22,)

Adding Constants

Note

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

Train Test Split

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,)

Model Selection and Training

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:                Mon, 04 Nov 2024   Prob (F-statistic):           1.84e-07
Time:                        01:37:43   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.

Interpreting P-values

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:

A low p-value (typically less than 0.05) suggests that you can reject the null hypothesis, meaning the coefficient is statistically significant and likely has an impact on the dependent variable.
A high p-value (greater than 0.05) indicates that the coefficient is not statistically significant, implying that the feature may not contribute meaningfully to the model.

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.924

The 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.143

Remember, these metrics only tell us about the predictive accuracy on the training data, and we care more about evaluating metrics on the test set.

Model Predictions and Evaluation

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:

mean: This is the predicted value for each observation based on the model. It’s essentially the point prediction (ŷ) for the corresponding input.
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:

Code

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='red',# 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.706

We see the same results on the test set that we did with the sklearn regression implementation.

Linear Regression with statsmodels

Data Loading

Data Splitting

X/Y Split

Adding Constants

Train Test Split

Model Selection and Training

Model Predictions and Evaluation

Linear Regression with `statsmodels`