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