The Illusion of Data: Unpacking Simpson’s Paradox Through Python Simulations

October 2026
By Zhenguo Zhang

In the realm of data science and statistical analysis, numbers are often revered as objective truths. However, a notorious statistical phenomenon known as Simpson’s paradox proves that without proper context, aggregated data can tell a story that is diametrically opposed to reality. Simpson’s paradox occurs when a trend or association that appears in different groups of data reverses itself when these groups are combined into a single, pooled dataset.

To demystify this counterintuitive concept, data scientist Zhenguo Zhang recently published a comprehensive, reproducible walkthrough on his blog, leveraging Python within an R Markdown environment. Utilizing popular data science libraries—including NumPy, Pandas, Statsmodels, and Plotnine—Zhang simulates an educational dataset of 20,000 students to demonstrate how a confounding variable can completely flip statistical outcomes.

Main Facts: The Anatomy of a Statistical Reversal

At its core, Simpson’s paradox exposes the dangers of ignoring confounding variables. Zhang’s simulation centers on a hypothetical academic scenario: evaluating the success rates of two study programs, designated as Program A and Program B, while factoring in the difficulty of their respective exams.

The simulation establishes a stark disparity in how students from each program choose their test paths:

• Program A students predominantly opt for the rigorous route, with 90% taking the hard exam and only 10% taking the easy exam.
• Program B students lean heavily toward the path of least resistance, with 90% choosing the easy exam and only 10% attempting the hard exam.

When examining the pass rates within specific subgroups, Program A consistently outperforms Program B regardless of the test difficulty:

• On the Easy Exam, Program A achieves an 88.6% pass rate compared to Program B’s 79.7%.
• On the Hard Exam, Program A maintains a 39.3% pass rate against Program B’s 31.0%.

Yet, when all 20,000 student records are pooled together into an aggregate overview, the narrative completely inverts. Program B boasts an overall pass rate of 74.8%, while Program A languishes at just 44.2%.

Chronology of the Simulation Workflow

To replicate these findings and understand how the statistical inversion manifests, Zhang structured his analysis as a sequential, reproducible workflow.

1. Environment Setup and Tooling

The analysis was executed using an R Markdown document that bridged R and Python via the reticulate package. This setup allowed for seamless execution of Python scripts within a shared session. The foundational libraries utilized include:

• NumPy: For setting a pseudo-random seed (rng = np.random.default_rng(2026)) and independently simulating binomial pass/fail outcomes (rng.binomial(1, p)).
• Pandas: For structuring, cleaning, and aggregating the dataset of 20,000 student records.
• Statsmodels: For fitting logistic regression models to evaluate the odds ratios before and after introducing a confounder.
• Plotnine: A Python implementation of the Grammar of Graphics, used to build a faceted visualization illustrating the paradox.

2. Dataset Generation

Zhang established fixed group sizes and pre-determined pass probabilities to construct the raw data frame:

• Easy / Program A: 1,000 students, 90% pass probability
• Easy / Program B: 9,000 students, 80% pass probability
• Hard / Program A: 9,000 students, 40% pass probability
• Hard / Program B: 1,000 students, 30% pass probability

3. Regression Modeling and Visualization

First, a baseline logistic regression model tested passed ~ C(program) without controls. Next, a second model incorporated the confounding variable: passed ~ C(program) + C(difficulty). Finally, data aggregation scripts and plotnine graphics were deployed to map out the side-by-side performance metrics across Easy, Hard, and Combined categories.

Supporting Data and Logistic Regression Results

The empirical evidence of the paradox is laid bare when comparing the unadjusted and adjusted logistic regression outputs generated by Statsmodels.

Unadjusted Model (Program Only)

When regressing student success solely against their program membership, the model yields a positive coefficient for Program B:

===================================================================================
                      coef    std err          z      P>|z|      [0.025      0.975]
-----------------------------------------------------------------------------------
Intercept          -0.2322      0.020   -11.534      0.000      -0.272      -0.193
C(program)[T.B]     1.3197      0.031    43.140      0.000       1.260       1.380
===================================================================================

The coefficient of 1.3197 on the log-odds scale indicates a statistically significant advantage for students enrolled in Program B.

Adjusted Model (Adding Exam Difficulty)

However, once exam difficulty is introduced as a covariate, the regression coefficients undergo a dramatic shift:

=========================================================================================
                            coef    std err          z      P>|z|      [0.025      0.975]
-----------------------------------------------------------------------------------------
Intercept                 1.8555      0.059     31.210      0.000       1.739       1.972
C(program)[T.B]          -0.4760      0.058     -8.185      0.000      -0.590      -0.362
C(difficulty)[T.Hard]    -2.2809      0.059    -38.892      0.000      -2.396      -2.166
=========================================================================================

Holding exam difficulty constant, the coefficient for Program B drops to -0.4760. The statistical reality is revealed: when comparing students taking exams of identical difficulty, Program A students are actually more likely to pass.

Weighted Averages Breakdown

The mathematical engine driving this paradox is the weighted arithmetic mean. Because the overall pass rate is a composite of subgroup performances weighted by sample size proportions, the uneven distribution skews the aggregate:

• Program A Overall Rate: $(0.10 times 0.90) + (0.90 times 0.40) = 44.22%$
• Program B Overall Rate: $(0.90 times 0.80) + (0.10 times 0.30) = 74.79%$

Even though Program A outperforms Program B at every granular level, Program A’s heavy concentration in the difficult exam tier drags its aggregate score down.

Methodological Implications: The Danger of Confounders

Zhang’s simulation offers critical lessons for data analysts, researchers, and decision-makers across industries—from healthcare and education to economics and tech.

1. The Perils of Aggregation: Relying purely on top-level metrics can lead to deeply flawed conclusions. If a school administrator only looked at overall pass rates, they would mistakenly praise Program B while penalizing Program A.
2. Confounding Variables: Exam difficulty acts as a classic confounder. It correlates with program choice while directly impacting the success probability. Failing to control for it distorts the true relationship between variables.
3. Causation vs. Correlation: A statistical reversal highlights the fact that simply observing data is not enough to establish causality. Analysts must maintain rigorous assumptions about how variables interact.
4. The Risk of Over-Adjustment: While controlling for confounders is necessary, blindly conditioning on every available variable—such as accidentally controlling for a collider variable—can introduce new statistical biases. Subgroup slicing is not a universal panacea and must be guided by sound domain knowledge and causal frameworks.

Ultimately, Zhenguo Zhang’s Python-based exploration serves as a timely reminder that data literacy requires looking beneath the surface. In a world increasingly driven by automated analytics and pooled datasets, understanding the mechanics of Simpson’s paradox is essential to ensuring that data illuminates the truth rather than obscuring it.