Regression coefficients do not match conditional means In a nutshell, I want the regression coefficients of a model to match several differences in conditional means.
You can download the data from this repo.

I have a data set that has a dependent variable (Y) and three binary columns (T, X1 and X2).





Y
CONST
T
X1
X1T
X2
X2T




0
2.31252
1
1
0
0
1
1


1
-0.836074
1
1
1
1
1
1


2
-0.797183
1
0
0
0
1
0




I want to calculate the difference in the mean of Y for observations with T == 1 and those with T == 0 for each of the four possible combinations of X1 and X2:

*

*Mean difference given X1 == 0 and X2 == 0

*Mean difference given X1 == 0 and X2 == 1

*Mean difference given X1 == 1 and X2 == 0

*Mean difference given X1 == 1 and X2 == 1
I did this manually, but I cannot get the following model to match my results:
$$Y = \beta_0 + \beta_1 T + \beta_2 X_1 + \beta_3 X_1 T + \beta_4 X_2 + \beta_5 X_2 T + U$$
As per this post:

*

*$\hat{\beta_1}$ should match case 1

*$\hat{\beta_1} + \hat{\beta_5}$ should match case 2

*$\hat{\beta_1} + \hat{\beta_3}$ should match case 3

*$\hat{\beta_1} + \hat{\beta_3} + \hat{\beta_5}$ should match case 4

As can be seen in this jupyter notebook, I cannot get these two methods to match.
How come the linear regression results do not match the differences in conditional means?

It's not straightforward to estimate effects by hand. The math works out nicely only in "simple" cases.
One simple case is a saturated model. A saturated model includes all main effects and possible interactions, so it has a parameter for each unique combination of the predictors.
Your model is not saturated; it's missing a three-way interaction T * X1 * X2 as well as a two-way interaction between the covariates X1 * X2. By omitting interactions you impose (implicit) constraints on the remaining parameters, so the estimate of $\beta_1$ is a function of not only the sample means for groups {T=1, X1=0, X2=0} and {T=0, X1=0, X2=0} but all group means.
This is easier to understand in a model with one covariate in addition to the treatment:
$$
Y = \beta_0 + \beta_1 T + \beta_2 X
$$
In this model $\beta_1$ is the treatment effect for both group X = 0 and group X = 1 as omitting the T * X interaction imposes the constraint that the two groups have the same treatment effect. To estimate $\beta_1$ the regression takes into account all observations, not only those with X = 0. Your model is more complex, so the implicit constraints are harder to conceptualize. (For example, combine cases 1&4 and cases 2&3 to show that the average treatment effect is the same when X1 = X2 and when X1 ≠ X2.)

My original answer considered a different special case altogether. I include it for completeness.
The sample means for the treatment and control groups (t = 1 and t = 0, respectively) are the same as the marginal means (the population means) from a linear regression adjusted for covariates only if the data is balanced. A balanced dataset has the same number of observations for each unique combination of the predictors.
set.seed(1234)

n <- 50
y <- rnorm(2 * n)

# balanced data
t <- rep(c(0, 1), each = n)
x <- rep(c(0, 1), times = n)
table(t, x)
#>    x
#> t    0  1
#>   0 25 25
#>   1 25 25

coef(lm(y ~ t + x))["t"]
#>         t 
#> 0.5925825
mean(y[t == 1]) - mean(y[t == 0])
#> [1] 0.5925825

# imbalanced data
t <- sample(c(0, 1), 2 * n, replace = TRUE)
x <- sample(c(0, 1), 2 * n, replace = TRUE)
table(t, x)
#>    x
#> t    0  1
#>   0 22 33
#>   1 26 19

coef(lm(y ~ t + x))["t"]
#>         t 
#> 0.1194514
mean(y[t == 1]) - mean(y[t == 0])
#> [1] 0.1402377

