How to prove equality of standard errors for two-sample t-test and linear regression?

Question

Using classical linear model assumptions, we know that $$\frac{\hat \beta_j - \beta_j}{se(\hat \beta_j)} \sim t_{n-k-1}$$ meaning that the ratio of regression coefficients to their standard error follows a t-distribution with $n - k - 1$ degrees of freedom, where $k$ is the number of slope parameters.

This means that the T-score from a 2-sample T-test for samples that have a different number of observations and different sample variance, $$T = \frac{\bar x_1 - \bar x_2}{\sqrt{\frac{s_1^2}{n_1} + \frac{s_2^2}{n_2}}},$$

will be the same as $$T = \frac{\hat \beta_1}{se(\hat \beta_1)}$$ for linear regression $Y = \beta_0 + \beta_1D$ where $D$ is a binary dummy variable that takes the values of 0 or 1 depending on if an observation is in group 2 or not.

Since $\hat \beta_1$ will represent the difference in means between the two groups, this indicates that the standard errors must be the same for the T-score to be equal.

So my question is: How can I prove the standard errors identical in this case? In other words, how do I mathematically prove: $${\sqrt{\frac{s_1^2}{n_1} + \frac{s_2^2}{n_2}}} = se(\hat \beta_1)$$

Answer 1

I do not think you can prove it because it is not generally true!

Essentially, your first formula allows variances to differ across groups ("heteroskdasticity") while the standard s.e. formula assumes homoskedasticity.

As this answer explains, the left hand side of your final formula is the heteroskedasticity robust standard error, while the default output for the standard error will be the square root of the [2,2] element of $$ s^2(X'X)^{-1}, $$ with $s^2$ the sum of squared residuals of the regression divided by $n-2$ (see e.g. here).

Using results from this answer, we have $$ (X'X)^{-1}=\frac{1}{nn_1-n_1^2}\begin{pmatrix}n_1&-n_1\\-n_1&n\end{pmatrix}, $$ where $n_1$ is the number of observations for which $D_i=1$. Note that $n=n_0+n_1$ and hence $nn_1-n_1^2=n_0n_1$. Thus, the [2,2] element is $$ s^2n/n_0n_1, $$ and the expression generally does not agree with the heteroskedastic variance estimator. [If we replaced $s^2$ with the ML estimator that divides the sum of squared residuals by $n$ and if $n_0=n_1$ we would have a case in which the two standard errors agree.]

Borrowing and extending code from the first linked answer for illustration:

library(sandwich)
library(lmtest)

n <- 10
y <- rnorm(n)                   # some dependent variable
x1 <- rbinom(n, size=1, p=.4)   # a dummy regressor

n1 <- length(y[x1==1])          # the no. of y's belonging to D_i=1
n0 <- length(y[x1==0])          # the no. of y's belonging to D_i=0

reg <- lm(y~x1)
prepackaged <- coeftest(reg, vcov = vcovHC(reg, "HC0"))[2,2]^2 # square to get variance instead of standard error

s1.squared <- sum((y[x1==1] - mean(y[x1==1]))^2) # sum of squared residuals belonging to first group
s0.squared <- sum((y[x1==0] - mean(y[x1==0]))^2) # sum of squared residuals belonging to second group
sum.of.variances <- s1.squared/n1^2 + s0.squared/n0^2

all.equal(prepackaged, sum.of.variances) # TRUE

# homoskedastic versions, not identical:
prepackaged.homoskedastic <- summary(reg)$coefficients[2,2]^2          
(s1.squared + s0.squared)*n/((n-2)*n0*n1)         # by hand

all.equal(prepackaged, prepackaged.homoskedastic) # FALSE