Causal inference method for analyzing randomized control trial with covariates / pre intervention observations

Question

I've got a seemingly easy situation, which turns out to be a little more complex than originally thought.
Here's the Setup:
We have a randomized controlled trial. Test and Control groups are the same size with large n. We have many measures about the individuals and know that they differ in those. However, we don't use any of those in treatment assignment (no stratification, etc.) Also, we have historic data for the target measure. During a period (here: post period), we apply an intervention to the test treatment and expect it to affect the measure. We know, that there is no interaction between individuals, so each observation is independent.
Now, what would be the "best" approach to conclude if the intervention was successful?

1. Of course, the first and simplest thing you might do, is apply a simple t-Test on the intervention period group means.
2. Next, you might wonder if there are better analyses that yield a higher power / precision. For example, we know the measure before the intervention. This sounds like it carries some information that we could use. So you might do a difference in difference approach. Here, you could take the post and pre period difference for the measure for both treatments and compare those means (again t-Test).
3. Another possibility would be to apply a regression analysis. Here, you could regress the treatment and the pre period values on the post values. Again, making use of the information in the pre period for a better inference.
4. Finally, you could also add an interaction term to 3. between treatment and period. This seems to be a somewhat standard approach in econometrics.

Here is my question as Python code with simulated data

import pandas as pd
import numpy as np
import statsmodels.formula.api as smf

# Create some data
n = 100
np.random.seed(10)
x_pre = pd.Series(np.random.normal(10, 2, n), name="pre")
y_pre = pd.Series(np.random.normal(10, 2, n), name="pre")
intervention = np.random.normal(5, 1, n)

# Post and pre are correlated
x_post = x_pre * np.random.normal(2, 1, n) 
y_post = y_pre * np.random.normal(2, 1, n)  + intervention

# data to analysis format
x = pd.concat([x_pre, x_post], axis=1, keys=["pre", "post"])
y = pd.concat([y_pre, y_post], axis=1, keys=["pre", "post"])
x["test"] = 0
y["test"] = 1
df = pd.concat([x, y]).reset_index(drop=True)
print(df.sample(4))

           pre       post  test
17   10.270274  18.831519     0
77   11.241201  11.581746     0
80   13.970169  19.358396     0
114   9.374342  18.699756     1

Let's visualize the data:

import seaborn as sns
sns.scatterplot(x="pre", y="post", hue="test", data=df)

Now, let's compare the different approaches:

# center pre data 
df["pre_centered"] = df["pre"] - df["pre"].mean()

FORMULAE = [
    "post ~ test",  # t-Test on Outcome Means
    "I(post - pre) ~ test",  # t-Test on Diff-in-Diff Outcome Means
    "post ~ pre + test",  # Add covariate for pre, account for pre differences
    "post ~ pre + test + pre * test",  # " + interaction
    "post ~ pre_centered + test + pre_centered * test",  # " + center pre
]
results = [smf.ols(formula=f, data=df).fit() for f in FORMULAE]
for r in results:
    print(r.summary())

/e: Added the interaction regression with centered pre as suggested by Noah in the comments.

Here's the output (slightly shortened for brevity):

                               OLS Regression Results                            
    ==============================================================================
    Dep. Variable:                   post   R-squared:                       0.024
    Model:                            OLS   Adj. R-squared:                  0.019
    No. Observations:                 200   AIC:                             1524.
    Df Residuals:                     198   BIC:                             1531.
    Df Model:                           1                                         
    Covariance Type:            nonrobust                                         
    ==============================================================================
                     coef    std err          t      P>|t|      [0.025      0.975]
    ------------------------------------------------------------------------------
    Intercept     21.2926      1.088     19.572      0.000      19.147      23.438
    test           3.4092      1.539      2.216      0.028       0.375       6.443
    ==============================================================================
    Omnibus:                        2.489   Durbin-Watson:                   2.227
    Prob(Omnibus):                  0.288   Jarque-Bera (JB):                2.095
    Skew:                           0.223   Prob(JB):                        0.351
    Kurtosis:                       3.229   Cond. No.                         2.62
    ==============================================================================
    

                                OLS Regression Results                            
    ==============================================================================
    Dep. Variable:          I(post - pre)   R-squared:                       0.027
    Model:                            OLS   Adj. R-squared:                  0.022
    No. Observations:                 200   AIC:                             1502.
    Df Residuals:                     198   BIC:                             1509.
    Df Model:                           1                                         
    Covariance Type:            nonrobust                                         
    ==============================================================================
                     coef    std err          t      P>|t|      [0.025      0.975]
    ------------------------------------------------------------------------------
    Intercept     11.1337      1.029     10.822      0.000       9.105      13.163
    test           3.4296      1.455      2.357      0.019       0.560       6.299
    ==============================================================================
    Omnibus:                        4.666   Durbin-Watson:                   2.266
    Prob(Omnibus):                  0.097   Jarque-Bera (JB):                6.319
    Skew:                          -0.028   Prob(JB):                       0.0424
    Kurtosis:                       3.869   Cond. No.                         2.62
    ==============================================================================
    

                                OLS Regression Results                            
    ==============================================================================
    Dep. Variable:                   post   R-squared:                       0.167
    Model:                            OLS   Adj. R-squared:                  0.159
    No. Observations:                 200   AIC:                             1495.
    Df Residuals:                     197   BIC:                             1504.
    Df Model:                           2                                         
    Covariance Type:            nonrobust                                         
    ==============================================================================
                     coef    std err          t      P>|t|      [0.025      0.975]
    ------------------------------------------------------------------------------
    Intercept     -0.2797      3.841     -0.073      0.942      -7.855       7.295
    pre            2.1235      0.365      5.820      0.000       1.404       2.843
    test           3.4526      1.425      2.423      0.016       0.643       6.262
    ==============================================================================
    Omnibus:                       17.035   Durbin-Watson:                   2.287
    Prob(Omnibus):                  0.000   Jarque-Bera (JB):               34.674
    Skew:                          -0.391   Prob(JB):                     2.96e-08
    Kurtosis:                       4.884   Cond. No.                         56.4
    ==============================================================================
    

                                OLS Regression Results                            
    ==============================================================================
    Dep. Variable:                   post   R-squared:                       0.175
    Model:                            OLS   Adj. R-squared:                  0.163
    No. Observations:                 200   AIC:                             1495.
    Df Residuals:                     196   BIC:                             1508.
    Df Model:                           3                                         
    Covariance Type:            nonrobust                                         
    ==============================================================================
                     coef    std err          t      P>|t|      [0.025      0.975]
    ------------------------------------------------------------------------------
    Intercept     -5.4464      5.375     -1.013      0.312     -16.046       5.154
    pre            2.6321      0.520      5.064      0.000       1.607       3.657
    test          13.5859      7.526      1.805      0.073      -1.257      28.429
    pre:test      -0.9985      0.728     -1.371      0.172      -2.435       0.438
    ==============================================================================
    Omnibus:                       14.283   Durbin-Watson:                   2.289
    Prob(Omnibus):                  0.001   Jarque-Bera (JB):               24.704
    Skew:                          -0.375   Prob(JB):                     4.32e-06
    Kurtosis:                       4.549   Cond. No.                         145.
    ==============================================================================
    
                            OLS Regression Results                            
==============================================================================
Dep. Variable:                   post   R-squared:                       0.175
Model:                            OLS   Adj. R-squared:                  0.163
No. Observations:                 200   AIC:                             1495.
Df Residuals:                     196   BIC:                             1508.
Df Model:                           3                                         
Covariance Type:            nonrobust                                         
=====================================================================================
                        coef    std err          t      P>|t|      [0.025      0.975]
-------------------------------------------------------------------------------------
Intercept            21.2657      1.005     21.155      0.000      19.283      23.248
pre_centered          2.6321      0.520      5.064      0.000       1.607       3.657
test                  3.4528      1.422      2.429      0.016       0.649       6.256
pre_centered:test    -0.9985      0.728     -1.371      0.172      -2.435       0.438
==============================================================================
Omnibus:                       14.283   Durbin-Watson:                   2.289
Prob(Omnibus):                  0.001   Jarque-Bera (JB):               24.704
Skew:                          -0.375   Prob(JB):                     4.32e-06
Kurtosis:                       4.549   Cond. No.                         5.13
==============================================================================

Here's some more specific questions:

1. What are the differences between these methods?
2. Which one is the most appropriate for this case?
3. Why do the P values vary so widely?
4. How would you interpret 4. (the interaction term regression)?

P.S:: I've already read a lot of blogs and papers regarding the general topic before posting this. However, there seem to be diverging opinions. (e.g. some people write that you shouldn't apply regression methods for RCTs because the assumptions are not satisfied, some people think that in most cases this is fine). So basically, this has confused be more than it has helped. Moreover, I have found many theoretical and general examples but only very few applied ones and none was exactly my case.

/e: This paper does a pretty similar comparison of methods. Unfortunately, their data is different as they have two follow up measurements.

Answer 1

Turns out, the paper (Twisk, J., Bosman, L., Hoekstra, T., Rijnhart, J., Welten, M., & Heymans, M. (2018)) I mentioned before has a lot of the answers I was looking for. Also, the paper (McKenzie, D. (2012)) mentioned by @Dimitry has been helpful. I'll share some of my insights from studying them more thoroughly:

The kind of randomized control trial or experiment I'm referring to can often be found in a medical context. That's why there are a lot of papers in medical journals dealing with similar cases. It is often called a pre/post study or a repeated measurement study. Gliner, J. A., Morgan, G. A., & Harmon, R. J. (2003) is a good start for a concise overview.

So, how should you analyze the result of such an experiment? It would be totally fine to just take the group means for your post measurement and compare those with a simple t-Test. But is this always the best strategy?
The answer seems to be: No!

Why is that?
Well, even when you randomize your groups there will be baseline differences between them. Because in expectation, the difference in outcomes will only depend on your intervention in the test group, this seems not to be a big issue (especially when your n is high). But it is a problem for your Power! If there are stark differences between characteristics of your individuals which are correlated with your outcome, you will have a harder time finding the effect of the intervention. Just by chance there will be cases where your randomization produces very unequal groups. Imagine having 20 persons (10 male / 10 female) to randomize into two groups. If you end up with a test group of 10 f and control of 10 m and sex is related to your outcome, you will have a bad time looking at your results. Another aspect to consider is "regression to the mean": groups with a low (high) measure at baseline are more likely to increase (decrease) their measure in the post period. This might happen in the absence of any intervention effect!
Moreover, baseline differences don't even have to be significant in order to be problematic. Twisk et al. argue that this is a huge misunderstanding and you should always account for them.

One solution can be stratification. By stratifying, you make sure that your groups end up equal. You reduce uninformative grouping outcomes and thereby variance. This increases Power.
Another solution is to account for baseline differences when your pre period measure is related to the post measure. You can do so by using appropriate inference methods. While there has been some debate on whether this should be done, this is mostly settled (Twisk et al.). However, many people are unsure which method is appropriate to deal with baseline differences (I was one of those).

So, which method is best for taking baseline differences into account and increasing the Power of your experiment?

I've turned my code from above into a simulation script. This has helped me to make sense of the more theoretical concepts outlined by Twisk et al. and especially by McKenzie.

One of my mistakes in the original post, was not taking into account the correct structure of the data. Let's correct this. Here is how the data looks:

|  id |      pre |     post |   test |
|----:|---------:|---------:|-------:|
|  1  |  8.31908 |  1.06574 |      0 |
|  2  |  9.06334 | -9.53055 |      0 |
| 100 | 10.4542  | 47.5967  |      1 |
| 101 | 12.6218  |  3.11528 |      1 |

This is the "wide" data format and represents cross-sectional data (even though we have an underlying time component). We apply the following inference methods to it:

FORMULAE = [
    "post ~ test",                # 0a t-Test on Outcome Means
    "post ~ test + pre",          # 1a cross-sectional reg. control for baseline
    "I(post - pre) ~ test",       # 3a t-Test on outcome change
    "I(post - pre) ~ test + pre", # 3b cross-sec. reg. with control for baseline
]

I've named the formulas according to the Twisk et al. paper for direct comparison. However, they didn't include the simple t-Test (0a). It will be interesting to see how this most naive approach compares to the other though. While you might think that 3a controls for baseline effects, it does not! You still need to add the baseline as a covariate, thus we add 3b. Actually, 3b is analogous to 1a. (see Twisk et al. for the derivation) The coefficient for test will be the Average Treatment Effect (ATE) in all cases.

For the upcoming methods, we need to adapt the data structure. This is what i didn't account for in my original post:

|      id |   test |   period |   value |
|--------:|-------:|---------:|--------:|
|       1 |      0 |        0 | 14.107  |
|       1 |      0 |        1 | -9.5547 |
|     100 |      1 |        0 |  8.9816 |
|     100 |      1 |        1 | 22.5591 |

Here, we really use the longitudinal / panel structure of the data. This is needed for the following methods:

FORMULAE = [
    "value ~ test + period + test * period",  # 2a panel regression with interaction
    "value ~ period + I(test * period)"       # 2c " without treatment covariate
]

These approaches can be helpful, when you have missing data. Subjects that have at least a baseline observation still contribute to the model. This is not the case with the previous approaches. Notice that 2a does not take baseline differences into account. Hence, 2c is introduced. (refer to Twisk et al. for more details) For 2a you need to calculate test + interaction coefficient for the ATE. For 2c the ATE is simply the interaction coefficient.

Here's the results. Cross-Sectional format data:

| formula                    |   auto_corr |      r_sq |   nobs |   df_resid |   df_model |   c_intercept |   p_intercept |   c_test |   p_test |     c_pre |         p_pre |
|:---------------------------|------------:|----------:|-------:|-----------:|-----------:|--------------:|--------------:|---------:|---------:|----------:|--------------:|
| post ~ test                |    0.505331 | 0.0163235 |    200 |        198 |          1 |    59.9287    |   6.83357e-56 |  5.15359 | 0.239359 | nan       | nan           |
| post ~ test + pre          |    0.505331 | 0.270734  |    200 |        197 |          2 |     0.0369226 |   0.519833    |  5.10506 | 0.195384 |   5.99582 |   1.25446e-07 |
| I(post - pre) ~ test       |    0.505331 | 0.0172487 |    200 |        198 |          1 |    49.94      |   8.34025e-47 |  5.14368 | 0.225567 | nan       | nan           |
| I(post - pre) ~ test + pre |    0.505331 | 0.209847  |    200 |        197 |          2 |     0.0369226 |   0.519833    |  5.10506 | 0.195384 |   4.99582 |   9.28722e-06 |

Panel format data:

| formula                               |   auto_corr |     r_sq |   nobs |   df_resid |   df_model |   c_intercept |   p_intercept |      c_test |     p_test |   c_period |    p_period |   c_test:period |   p_test:period |   c_i(test * period) |   p_i(test * period) |
|:--------------------------------------|------------:|---------:|-------:|-----------:|-----------:|--------------:|--------------:|------------:|-----------:|-----------:|------------:|----------------:|----------------:|---------------------:|---------------------:|
| value ~ test + period + test * period |    0.505331 | 0.713389 |    400 |        396 |          3 |       9.9887  |   2.01945e-08 |   0.0099174 |   0.923874 |    49.94   | 8.7505e-54  |         5.14368 |        0.237087 |            nan       |           nan        |
| value ~ period + I(test * period)     |    0.505331 | 0.713379 |    400 |        397 |          2 |       9.99366 |   2.26815e-14 | nan         | nan        |    49.935  | 1.78043e-65 |       nan       |      nan        |              5.15359 |             0.159908 |

What are the main insights?

1. When you have a pre/post experiment and a baseline for your measure, account for it!
2. How well the methods perform strongly depends on the (auto)correlation of the data. Especially the p-value varies greatly, while coefficients are somewhat more stable. With low correlation between pre and post (<0.1) there is almost no difference. With high correlation (>0.5) the methods differ strongly. (fits the main findings of McKenzie)
3. There are large power gains to be had when accounting for the baseline. Especially when the measure has high correlation over the time dimension. (see Kahan, B. C., Jairath, V., Doré, C. J., & Morris, T. P. (2014)
4. Method 1a seems to be a good choice all over.
5. You can and should (in many cases) add additional covariates in a similar fashion. However, adding the baseline is the most important one. (see Kahan et al. 2014)
6. All this only hold when you have randomized groups. In observational studies you must not control for the baseline like this! (see Twisk et al.)