The OLS (Ordinary Least Squares) Paradox

Question

This is a concept in regression I have never been able to understand.

Suppose we have this model with an interaction term:

$$y = b_0 + b_1 \cdot \text{{age}} + b_2 \cdot \text{{gender}} + b_3 \cdot (\text{{age}} \cdot \text{{gender}})$$

where gender is a binary variable (1 for male and 0 for female).

Case 1: When fitting this model via OLS on the entire dataset, this will actually produce 2 mini models:

• For males (gender = 1):

$$y = b_0 + b_1 \cdot \text{{age}} + b_2 + b_3 \cdot \text{{age}}$$

$$y = (b_0' + b_2') + (b_1' + b_3') \cdot \text{{age}}$$

• For females (gender = 0):

$$y = b_0'' + b_1'' \cdot \text{{age}}$$

Case 2: Now, suppose we remove all males from the data and fit the model. Then, we remove all females from the data and fit the model again. The two resulting models will be:

• For males:

$$y = b_{0m} + b_{1m} \cdot \text{{age}}$$

• For females:

$$y = b_{0f} + b_{1f} \cdot \text{{age}}$$

For me, the paradox has always been trying to understand that: the model for Males in Case 1 is the same as for Males in Case 2, and the model for Females in Case 1 is the same for Females in Case 2 (i.e. the parameter estimates are the same). That is:

$$(b_0' + b_2') = b_{0m}$$ $$(b_1' + b_3') = b_{1m}$$ $$b_0'' = b_{0f}$$ $$b_1'' = b_{1f}$$

$$Var(b_0' + b_2') = Var(b_{0m})$$ $$Var(b_1' + b_3') = Var(b_{1m})$$ $$Var(b_0'') = Var(b_{0f})$$ $$Var(b_1'') = Var(b_{1f})$$

I have heard many intuitive explanations that try to logic out why this is true, but I have never been satisfied with these explanations. Therefore, I tried to write the estimating equations (based on either OLS or MLE) to try and see that if the estimating equations are same for Case 1 and Case 2, then the parameter estimates should also be the same for Case 1 and Case 2. (paradox)

I know that in OLS, the following estimating equations hold for linear regression ($n$ is the number of data points, $p$ is the number of predictor variables):

$$\hat{y} = \mathbf{x}^T \hat{\beta}$$ $$\hat{B} = (X'X)^{-1}X'Y$$ $$Var(\hat{\beta}) = \hat{\sigma}^2 (X'X)^{-1}$$ $$\hat{\sigma} = \sqrt{\frac{1}{n - p}\sum_{i=1}^{n}(y_i - \hat{y}_i)^2}$$

And we can write the likelihood for this regression model as :

$$ L(\beta | y, X) = \prod_{i=1}^{n} \frac{1}{\sqrt{2\pi\sigma^2}} \exp\left(-\frac{(y_i - (b_0 + b_1 \cdot \text{{age}}_i + b_2 \cdot \text{{gender}}_i + b_3 \cdot (\text{{age}}_i \cdot \text{{gender}}_i)))^2}{2\sigma^2}\right) $$

But it is still not clear to me that the estimates from Case 1 will be the same as the estimates from Case 2. Is it possible to see that only from the estimating equations perspective (OLS or MLE), is it possible to see that the parameter estimates and variance of the parameter estimates will be identical whether estimated from Case 1 or Case 2?

• Note 1: (My attempt at) Variance Derivation for OLS

  $$ (AB)^T = B^T A^T $$

  $$ (A^{-1})^T = (A^T)^{-1} $$

  $$ (A^T)^T = A $$

$$\begin{align*} \text{Var}(\hat{B}) &= (X^tX)^{-1} X^t \text{Var}(y) [(X^tX)^{-1} X^t]^T \\ &= \hat{\sigma}^2 (X^tX)^{-1} X^t [(X^t)^t ((X^tX)^{-1})^t] \\ &= \hat{\sigma}^2 (X^tX)^{-1} X X ((X^tX)^t)^{-1} \\ &= \hat{\sigma}^2 (X^tX)^{-1} X X (X^tX)^{-1} \\ &= \hat{\sigma}^2 (X^tX)^{-1} (X^tX) (X^tX)^{-1} &=\hat{\sigma}^2 (X'X)^{-1} \end{align*}$$

• Note 2: In the absence of an interaction term, Case 1 (fitting the model on the entire dataset) and Case 2 (fitting the model on each gender separately) will produce different estimates for the slopes and intercepts (i.e. paradox does not happen):

$$y = b_0 + b_1 \cdot \text{{age}} + b_2 \cdot \text{{gender}} $$

Answer 1

You can show the equivalence between the $\hat\beta$ of these models even more directly if you parametrize your case 1 model like so (removing the overall intercepts):

$$ y=[\text{sex}=f]\cdot b_{0f}+[\text{sex}=m]\cdot b_{0m}+[\text{sex}=f]\cdot \text{age}\cdot b_{1f}+[\text{sex}=m]\cdot \text{age}\cdot b_{1m} $$

The estimated values of these $b_{x}$ match one-to-one with your second case and the model will still return identical predictions to your original parametrization. An example in R:

## Simulate some data
set.seed(1)
sex <- rep(0:1, each=20)
age <- sample(18:65, size=length(sex), replace=TRUE)

y <- 10 + 1.5*age + 3*sex - 0.5*age*sex + rnorm(length(sex), sd=5)

## Your first model
case1 <- lm(y ~ age*sex)

## Alternative parametrization
bf <- 0+(sex==0)
bm <- 0+(sex==1)

case1b <- lm(y ~ 0 + bf + bm + bf:age + bm:age)

## Just to confirm that these are indeed equal predictions, including variances
all.equal(
  predict(case1, newdata=list(sex = c(0, 1), age = c(15, 33)), se.fit=TRUE),
  predict(case1b, newdata=list(bf = c(1, 0), bm = c(0, 1), age = c(15, 33)), se.fit=TRUE)
)
> TRUE

If you then compare the coefficients to the separate case 2 models, you don't even need to recalculate anything.

## Case 2 - subset outside the model or `predict()` will break
fi <- which(sex==0)
yf <- y[fi]
agef <- age[fi]

mi <- which(sex==1)
ym <- y[mi]
agem <- age[mi]

case2f <- lm(yf ~ agef)
case2m <- lm(ym ~ agem)

cbind(coef(case1b)[c(1, 3, 2, 4)], c(coef(case2f), coef(case2m)))
>          [,1]   [,2]
> bf     10.369 10.369
> bf:age  1.481  1.481
> bm     20.048 20.048
> bm:age  0.880  0.880

However...

all.equal(
  predict(case1b, newdata=list(bf = c(0, 0), bm = c(1, 1), age = c(15, 33)), se.fit=TRUE),
  predict(case2m, newdata=list(agem = c(15, 33)), se.fit=TRUE)
)
> [1] "Component “se.fit”: Mean relative difference: 0.03439107" ...

The reason for this is that the $y_i$ vector that went into the first model is not the same length (here twice as long) as in the other models. For that reason, even if $\bar y$ is the same, $\sum(y_i-\bar y)^2$ is not. This will cause the residual error $\hat\sigma$ between the models to differ, and thereby $\mathrm{Var}(\hat\beta)$ as well. We can manually recalculate some of these quantities (but please note that the real computations are implemented much more efficiently & safer):

## Calculate design variance (X'X)^-1
xpxi <- function(x) solve(t(x) %*% x)
## Calculate residual error sigma
residual_error <- function(o) sqrt(sum(o$residuals**2) / o$df.residual)

## Combine the above into standard errors for beta
## Compare this also to the printout of `summary(model)`
model_se <- function(m) sqrt(residual_error(m)**2 * diag(xpxi(model.matrix(m))))

cbind(model_se(case1b)[c(2,4)], model_se(case2m))
>          [,1]   [,2]
> bm     4.1099 4.2513
> bm:age 0.0845 0.0874

It is the length of $y_i$ (i.e $n$ and the dimensions of $X$) that is driving most of these differences - another small discrepancy exists in the number of parameters between the two cases. As an even more extreme example you could take data from one sex and replicate it, the mean will stay the same but look what happens to the variance estimates:

yx <- rep(ym, each=2)
agex <- rep(agem, each=2)

case3 <- lm(yx ~ agex)
cbind(model_se(case2m), model_se(case3))
>       [,1]   [,2]
> Int 4.2513 2.9259
> age 0.0874 0.0601

The point in your original example is that the residual error is calculated from the entire model, it is the weighted sum of male and female residuals in case 1.

Something I didn't go into here is that the covariance of your parameters will change as well, not just the diagonal elements of $\mathrm{Var}(\hat\beta)$, so all of this will also apply to contrasts and hypothesis tests between model parameters.