I'm interested in testing if a particular (binary) feature is significant in explaining a target variable, after other factors have been considered.

Naively, I can do this by including the test feature in a regression model (alongside the other 'explaining' features), and then look at the coefficient and p-value associated with the test feature in the trained model.

Alternatively, I can build the regression model without including the test feature. Then I can then use that fitted model to make predictions for the target, then calculate residuals between the prediction and the true target values. Finally, I can t-test those residuals to see if there is a significant difference between the mean for group 0 and the mean for group 1.

The advantage of the second method is that if I'm happy with my fitted regression model, I can use it to test multiple different datasets, without having to refit a model each time.

My question is, are both the methods valid, and how can I expect the results (in particular the p-value) to differ in the two different scenarios?

I've included some code below to demonstrate the two methods:

from sklearn import datasets
import statsmodels.api as sm
import pandas as pd
from scipy.stats import ttest_ind
import numpy as np

iris = datasets.load_boston()
X = iris.data
y = iris.target

In the Boston dataset y is house price and X is a range of features to explain house prices. The 4th feature in X (X[:,3]) is in fact a binary variable - so we will use that as our test feature.

model = sm.OLS(y, X)
fitted_model = model.fit()

...which yeilds...

    Coef.   Std.Err.    t   P>|t|   [0.025  0.975]
x1  -0.0916 0.0343  -2.6747 0.0077  -0.1589 -0.0243
x2  0.0487  0.0144  3.3791  0.0008  0.0204  0.0770
x3  -0.0038 0.0644  -0.0586 0.9533  -0.1304 0.1228
x4  2.8564  0.9040  3.1597  0.0017  1.0802  4.6325
x5  -2.8808 3.3593  -0.8575 0.3916  -9.4812 3.7196
x6  5.9252  0.3091  19.1679 0.0000  5.3179  6.5326
x7  -0.0072 0.0138  -0.5229 0.6013  -0.0344 0.0199
x8  -0.9680 0.1957  -4.9475 0.0000  -1.3524 -0.5836
x9  0.1704  0.0667  2.5541  0.0109  0.0393  0.3016
x10 -0.0094 0.0039  -2.3930 0.0171  -0.0171 -0.0017
x11 -0.3924 0.1099  -3.5713 0.0004  -0.6083 -0.1765
x12 0.0150  0.0027  5.5611  0.0000  0.0097  0.0203
x13 -0.4170 0.0508  -8.2142 0.0000  -0.5167 -0.3172

So we can see x4 (our test feature) has a coefficient of 2.8564 with p-val of 0.0017.

Alternatively, we can use the t-test on residuals method.

# pull out test feature and cast to boolean
test_feature = X[:,3].astype(bool)

# make a design matrix without X[:,3]
X_without_test_feature = np.hstack([X[:,0:3], X[:,4:]])

model2 = sm.OLS(y, X_without_test_feature)
fitted_model2 = model2.fit()

y_pred = fitted_model2.predict(X_without_test_feature)
residuals = y-y_pred
group1 = residuals[test_feature]
group0 = residuals[~test_feature]

ttest_ind(group1, group0)

...which yields...

Ttest_indResult(statistic=3.00211771557846, pvalue=0.002814120078709997)

...indicating a difference between the average residuals of 3.002, with a p-value of 0.0028.

To repeat the question - are both the methods valid, and how can I expect the results (in particular the p-value) to differ in the two different scenarios? In this case we get very similar results - is this always likely to be the case?

Thanks in advance.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.