Introduction

In this example, we simulate a study designed to test the effect of two different food types (A vs. B) on body weight. Our sample is split into two groups: a control group (food B) and a treatment group (food A). Each group has the same sample size (1000 individuals).

We artificially set a “true” treatment effect of 5 kg, meaning that on average, individuals in the treatment group gain 5 kg more than those in the control group (all else being equal). Random noise is then added to reflect individual variation.

Key Point:

The design flaw is that we later subdivide participants by their final weight (the outcome) into “Low,” “Medium,” and “High” categories and then compare the mean weight in each category for treatment vs. control. This introduces serious bias because we are stratifying on an outcome variable that is itself affected by the treatment.

Let’s walk through each step of the simulation and see how post-treatment stratification can create misleading results.

Step 1: Parameter Setup

# Sample sizes
n_per_group <- 1000  # individuals in each group

# True effect (in kilograms)
true_treatment_effect <- 5

# Standard deviation of random noise
sd_noise <- 3

n_per_group: number of individuals in each group (treatment vs. control).

true_treatment_effect: the true average weight gain (in kg) caused by the treatment (food A).

sd_noise: random noise to reflect natural variation in final weight.

Step 2: Generate Baseline Weights

We assume each individual’s baseline weight is normally distributed with a mean of 70 kg and a standard deviation of 10 kg, reflecting typical adult body weights.

# Create baseline weights
baseline_weight <- rnorm(2 * n_per_group, mean = 70, sd = 10)

# Treatment assignment: 0 = control, 1 = treatment
treatment <- rep(c(0, 1), each = n_per_group)

# Combine into initial data frame
study_data <- data.frame(
  ID = 1:(2 * n_per_group),
  Treatment = factor(treatment, labels = c("Control", "Treatment")),
  Baseline_Weight = baseline_weight
)

print(head(study_data))
##   ID Treatment Baseline_Weight
## 1  1   Control        63.77649
## 2  2   Control        72.67767
## 3  3   Control        67.40899
## 4  4   Control        94.17652
## 5  5   Control        79.26178
## 6  6   Control        75.72210
print(tail(study_data))
##        ID Treatment Baseline_Weight
## 1995 1995 Treatment        85.66696
## 1996 1996 Treatment        71.28655
## 1997 1997 Treatment        76.18052
## 1998 1998 Treatment        63.18923
## 1999 1999 Treatment        76.77395
## 2000 2000 Treatment        65.01933

Step 3: Generate Final Weights

Final weight is determined by: \[ Final\ Weight = Baseline\ Weight + (Treatment\ Effect \times Treatment\ Assignment) + Noise \]

If an individual is in the treatment group, they receive an extra true_treatment_effect kg on average.

We add random noise (sd_noise = 3) to reflect individual variation (e.g., genetics, lifestyle).

final_weight <- study_data$Baseline_Weight +
  treatment * true_treatment_effect +
  rnorm(2 * n_per_group, mean = 0, sd = sd_noise)

study_data$Final_Weight <- final_weight

You should see the distribution shifted upward compared to baseline, particularly in the treatment group.

Step 4: Correct (Unconditional) Analysis

The “correct” way to measure the average treatment effect (assuming no confounders and purely random assignment) is to simply compare the mean final weight of the treatment group vs. the control group:

# Compare means by treatment
correct_analysis <- aggregate(Final_Weight ~ Treatment, data = study_data, FUN = mean)
names(correct_analysis)[2] <- "Mean_Final_Weight"
correct_analysis
##   Treatment Mean_Final_Weight
## 1   Control          70.00291
## 2 Treatment          74.76488

We expect the difference in final weight between the groups to be close to 5 kg (our true_treatment_effect), though random noise may cause slight deviations.

Visualizing the Correct Analysis

We can illustrate this difference with a simple boxplot:

## # A tibble: 2 × 2
##   Treatment Mean_Final_Weight
##   <fct>                 <dbl>
## 1 Control                70.0
## 2 Treatment              74.8

Here we see the shift in final weight for those who received the treatment compared to the control group.

Step 5: (Flawed) Post-Treatment Stratification

Now we demonstrate the flaw: subdividing individuals based on their final weight (the outcome). Let’s define cutpoints for “Low,” “Medium,” and “High” final weights.

study_data$Category <- cut(
  study_data$Final_Weight,
  breaks = c(-Inf, 70, 80, Inf),
  labels = c("Low", "Medium", "High")
)

Notice that if the treatment does increase final weight, a larger fraction of treatment-group individuals may shift into “Medium” or “High,” while more control-group individuals remain “Low” or “Medium.” Once we condition on these categories, the randomization is broken because the categories are dependent on treatment status.

Summaries by Post-Treatment Category

# Compare mean final weight by Category and Treatment
flawed_analysis <- aggregate(Final_Weight ~ Treatment + Category,
                             data = study_data,
                             FUN = mean)

names(flawed_analysis)[3] <- "Mean_Final_Weight_by_Category"
flawed_analysis
##   Treatment Category Mean_Final_Weight_by_Category
## 1   Control      Low                      61.56383
## 2 Treatment      Low                      63.29160
## 3   Control   Medium                      74.68398
## 4 Treatment   Medium                      74.83197
## 5   Control     High                      85.60876
## 6 Treatment     High                      86.55938

In some categories, you might even see the control group has a higher mean final weight than the treatment group. This seems paradoxical given that we know the true effect is +5 kg. This paradox arises because once you slice by final weight categories—determined by the outcome itself—you lose the protection of randomization.

Step 6: Visualizing the Flawed Comparisons

Here’s a bar plot to illustrate the final mean weight within each post-treatment category, split by treatment:

Interpretation:

We artificially segmented the data based on a variable impacted by the treatment. This can easily mask or even reverse the true difference in some categories because we are no longer comparing like for like. The very act of slicing by final weight means we’re mixing up the effect of the treatment with each individual’s natural weight variance.

Conclusion

Correct Approach: Compare all individuals in the treatment group vs. all in the control group (i.e., unconditional means). This will yield a difference close to the true effect of 5 kg.

Flawed Approach: Stratify by final weight (the outcome) and compare treatment vs. control within each stratum. Doing this post-outcome conditioning breaks randomization and introduces significant bias, leading to results that can be highly misleading.

If you run this simulation multiple times (by changing set.seed() or removing it altogether), you will see that, although the “correct” analysis consistently shows around a 5 kg difference, the “flawed” analysis can vary wildly and often underestimates (or even reverses) the true effect in certain categories.

sessionInfo()
## R version 4.4.2 (2024-10-31 ucrt)
## Platform: x86_64-w64-mingw32/x64
## Running under: Windows 11 x64 (build 26100)
## 
## Matrix products: default
## 
## 
## locale:
## [1] LC_COLLATE=English_Canada.utf8  LC_CTYPE=English_Canada.utf8   
## [3] LC_MONETARY=English_Canada.utf8 LC_NUMERIC=C                   
## [5] LC_TIME=English_Canada.utf8    
## 
## time zone: America/Edmonton
## tzcode source: internal
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] ggplot2_3.5.1 dplyr_1.1.4  
## 
## loaded via a namespace (and not attached):
##  [1] vctrs_0.6.5       cli_3.6.3         knitr_1.49        rlang_1.1.4      
##  [5] xfun_0.49         generics_0.1.3    jsonlite_1.8.9    labeling_0.4.3   
##  [9] glue_1.8.0        colorspace_2.1-1  htmltools_0.5.8.1 sass_0.4.9       
## [13] scales_1.3.0      rmarkdown_2.29    grid_4.4.2        munsell_0.5.1    
## [17] evaluate_1.0.1    jquerylib_0.1.4   tibble_3.2.1      fastmap_1.2.0    
## [21] yaml_2.3.10       lifecycle_1.0.4   compiler_4.4.2    pkgconfig_2.0.3  
## [25] rstudioapi_0.17.1 farver_2.1.2      digest_0.6.37     R6_2.5.1         
## [29] utf8_1.2.4        tidyselect_1.2.1  pillar_1.10.0     magrittr_2.0.3   
## [33] bslib_0.8.0       withr_3.0.2       gtable_0.3.6      tools_4.4.2      
## [37] cachem_1.1.0