2
$\begingroup$

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


df = pd.read_csv('my_file.txt')

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

OLS Regression Results

Dep. Variable:               response   R-squared:                       0.259
Model:                            OLS   Adj. R-squared:                  0.241
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.

Is this about right?

$\endgroup$
  • 2
    $\begingroup$ Looking closely at the reported degrees of freedom in the tests will answer your question. $\endgroup$ – whuber Jan 9 '15 at 18:11
  • $\begingroup$ Thanks @whuber. The df for both t-tests in the regression would be 40. Correct? The df for the one-sample t-test is 13, meaning that the test of the intercept is using a critical value based on all the data. I also surmise from this answer that the standard errors are computed from the entire design matrix, hence different to the one-sample case. Still don't fully understand why it holds for the independent-samples case (2nd coef), though. $\endgroup$ – tsawallis Jan 12 '15 at 13:46
  • $\begingroup$ I guess because the independent samples t-test is using the full design matrix (or exactly equivalent), whereas the one-sample t-test does not. Hence why standard errors are the same for the t-test, different for the one-sample? $\endgroup$ – tsawallis Jan 12 '15 at 13:47

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.