Equivalence of intercept beta test to one-sample t-test in linear regression model with categorical variable

I am testing my understanding of the equivalence between basic linear regression with categorical variables, and one-sample / independent samples t-tests. I don't think this corresponds to an existing question, though see here. Here is an example (comma separated) dataset I will use to demonstrate my problem:

response,cat_var
8.3,low
3.9,high
9.8,high
19.7,low
6.2,high
14.5,low
9.7,low
12.5,high
3.9,high
8.4,low
8.6,low
10.7,low
10.7,low
3.8,high
14.4,low
6.2,low
13.8,low
8.2,high
2.9,high
9.1,high
19.9,low
8.0,low
13.2,low
12.8,low
19.4,low
17.4,low
17.0,low
10.9,high
13.7,low
8.4,high
8.4,high
13.2,low
9.8,low
9.6,low
6.6,low
5.4,low
3.8,high
13.3,low
5.4,low
15.8,low
10.3,high
19.4,low


It consists of a continuous response variable and a categorical variable with two levels, 'high' and 'low'. Fitting this as a linear model using standard dummy coding results in a design matrix with a column of 1s (the intercept) and an indicator column, 0 for "high" and 1 for "low".

Code examples are in Python using Statsmodels, but that shouldn't be important here.

Fit the regression:

import numpy as np
import pandas as pd
import seaborn as sns
import statsmodels.formula.api as smf
import statsmodels.stats as sm_stats

fit = smf.ols(formula='response ~ cat_var', data=df).fit()
print(fit.summary())


OLS Regression Results

Dep. Variable:               response   R-squared:                       0.259
Method:                 Least Squares   F-statistic:                     13.99
Date:                Fri, 09 Jan 2015   Prob (F-statistic):           0.000575
Time:                        17:05:42   Log-Likelihood:                -117.87
No. Observations:                  42   AIC:                             239.7
Df Residuals:                      40   BIC:                             243.2
Df Model:                           1
Covariance Type:            nonrobust
==================================================================================
coef    std err          t      P>|t|      [95.0% Conf. Int.]
----------------------------------------------------------------------------------
Intercept          7.2929      1.097      6.649      0.000         5.076     9.510
cat_var[T.low]     5.0250      1.343      3.741      0.001         2.310     7.740
==============================================================================
Omnibus:                        2.740   Durbin-Watson:                   1.526
Prob(Omnibus):                  0.254   Jarque-Bera (JB):                1.557
Skew:                           0.165   Prob(JB):                        0.459
Kurtosis:                       2.117   Cond. No.                         3.23
==============================================================================


Now since the "high" group is coded with zeros in the design matrix, the Intercept parameter represents the mean of that group:

high_array = df['response'][df['cat_var']=='high']
low_array = df['response'][df['cat_var']=='low']

high_array.mean()


7.2928571428571436

The cat_var[T.low] coefficient represents the average difference between the groups. That is, it's the difference being tested in an independent groups t-test. Sure enough, if we run an independent groups t-test we get the same answer as the t-test of the regression coefficient:

t, p, df = sm_stats.weightstats.ttest_ind(low_array, high_array)
'Independent samples t-test: t({df}) = {t}, p = {p}'.format(
t=np.round(t, decimals=6), df=int(df), p=p)


'Independent samples t-test: t(40) = 3.740823, p = 0.0005752964234035196'

'Model values: t = {t}, p = {p}'.format(
t=np.round(fit.tvalues[1], decimals=6),
p=fit.pvalues[1])


'Model values: t = 3.740823, p = 0.0005752964234035218'

Great. That's the same (up to high precision), as expected.

The actual problem

Now I expect that the test of the intercept coefficient in the regression is equivalent to a one-sample t-test: that the mean of the "high" group is different from zero.

# using statsmodels, must first compute "descrstats":
d1 = sm_stats.weightstats.DescrStatsW(high_array)
t, p, df = d1.ttest_mean(0)
'One-sample t-test: t({df}) = {t}, p = {p}'.format(
t=np.round(t, decimals=6), df=int(df), p=p)


'One-sample t-test: t(13) = 8.632037, p = 9.6339663195598e-07'

'Model values: t = {t}, p = {p}'.format(
t=np.round(fit.tvalues[0], decimals=6),
p=fit.pvalues[0])


'Model values: t = 6.649277, p = 5.791887110561801e-08'

Hrm. Not the same. The test statistic in the regression model (6.64) is lower than the test statistic in the one-sample t-test (8.63). The value of the coefficient corresponds to the mean difference between the groups:

low_array.mean() - high_array.mean()


5.0249999999999977

So it's not that the coefficient is wrong – it's that the standard error of the coefficient (mean) is higher in the regression than in the one sample t-test. This results in a smaller test statistic.

What's going on here? I suspect that in the regression, the standard error of both coefficients depend on the total variance in the dataset. Whereas, in the one sample t-test, only the variance in the "high" category is taken into account (because this test doesn't know anything about the "low" category). For the independent samples t-test both groups enter into the comparison so the result is exactly equivalent to the regression.