# The OLS (Ordinary Least Squares) Paradox

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}}$$

• These models might produce the same $\hat\beta$ (your first four equalities) but $\mathrm{Var}(\hat\beta)$ is most decidedly not the same - in particular because the $y_i$s that are actually in the model, $n$, $p$, and therefore $\hat\sigma$ differ. Commented Jan 13 at 12:30
• is it possible to see why this might be from a mathematical perspective? for example - for models with the interaction term, in case 1 vs case 2: by writing the estimating equations in matrix form, can we see that the parameter estimates will be the same but the variance estimates will be different? Commented Jan 13 at 16:00
• From a mathematical perspective, you see it directly from the dimension of X. In case 1, X is of dimension T1x4 whereas in case 2 X is of dimension T2/3x2. First, the different number of observations and regressors will lead to different estimates for the variance sig2. Second, you neglect the covariances between the intercept, age, gender, and the interaction term. Commented Jan 13 at 16:28
• thanks for this tip. can you please post an answer that shows some steps if you have time? Commented Jan 13 at 19:40
• oh wait, I think I understand. the values of $n$ and $p$ are enough to change the variance estimates in case1 vs case2? Commented Jan 13 at 19:46

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)

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.