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

Question

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
Model:                            OLS   Adj. R-squared:                  0.677
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?

Answer 1

You are fitting an additive model, so the fitted values will not match the sample means exactly. It takes a 4 degree of freedom model to match the cell means. An additive model only has 3 degrees of freedom. The additive model imposes the constraint that E[y|x1=1, x2=z] - E[y|x1=0, x2=z] does not depend on z. This constraint will generally not hold exactly in the sample means.