# Difference between regression p-values and t-test on residuals?

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

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()
fitted_model.summary()


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