# How do I interpret the coefficients of a statsmodels multiple regression model with 2 dichotomous variables?

I have a continuous dependent variable Y and 2 dichotomous, crossed grouping factors forming 4 groups: A1, A2, B1, and B2. I am looking for the main effects of either factor, so I fit a linear model without an interaction with statsmodels.formula.api.ols Here's a reproducible example:

np.random.seed(12312)

means = {
'A1': 5,
'A2': 6,
'B1': 3,
'B2': 4
}

N = 20
var = .85

y = []
x1 = []
x2 = []

for k, v in means.items():
y.append(np.random.normal(loc=v, scale=var, size=N))
x1.append([int(k[0]=='A') for i in range(N)])
x2.append([int(k[1]=='1') for i in range(N)])

y = np.concatenate(y)
x1 = np.concatenate(x1)
x2 = np.concatenate(x2)

data = np.stack([y,x1,x2], axis=1)

df = pd.DataFrame(data, columns=['y','x1','x2'])
df.loc[:, 'x1'] = df.x1.astype(int); df.loc[:, 'x2'] = df.x2.astype(int)

lm = ols('y ~ x1 + x2', data=df).fit()


And here is the results summary given by print(lm.summary()):

                            OLS Regression Results
==============================================================================
Dep. Variable:                      y   R-squared:                       0.686
Method:                 Least Squares   F-statistic:                     83.93
Date:                Tue, 03 Dec 2019   Prob (F-statistic):           4.53e-20
Time:                        11:51:53   Log-Likelihood:                -98.488
No. Observations:                  80   AIC:                             203.0
Df Residuals:                      77   BIC:                             210.1
Df Model:                           2
Covariance Type:            nonrobust
==============================================================================
coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
Intercept      4.0107      0.164     24.518      0.000       3.685       4.336
x1             2.1321      0.189     11.288      0.000       1.756       2.508
x2            -1.2013      0.189     -6.360      0.000      -1.577      -0.825
==============================================================================
Omnibus:                        0.639   Durbin-Watson:                   2.121
Prob(Omnibus):                  0.727   Jarque-Bera (JB):                0.761
Skew:                           0.185   Prob(JB):                        0.684
Kurtosis:                       2.699   Cond. No.                         3.19
==============================================================================



I can see that the main effects I generated are there, but I am not sure how exactly to interpret the coefficients. Intuitively, the intercept term should be precisely the mean of the reference category (x1=0; x2=0), but looking at the group means, it is not:

x1 x2
0  0   4.090842
1   2.729360
1  0   6.062789
1   5.021698


And the difference (between the intercept and the mean) is even more pronounced when I work with the real data.

Since I cannot interpret the intercept coefficient, I am not sure whether the other two coefficients represent group differences in relation to the reference category.

I noticed, that when an interaction is included (e.g. lm = ols('y ~ x1 * x2', data=df).fit(), the intercept coefficient becomes precisely the mean of the reference category, and all other coefficients correspond to group differences. So what are the coefficients when the interaction is not included?