Stop Wasting Marketing Budget on Correlation: A Technical Guide to Causal Simulation with PyCausalSim
From correlation traps to counterfactual insights — a deep dive into causal inference for marketing analytics
From correlation traps to counterfactual insights — a deep dive into causal inference for marketing analytics
If you found this useful, you can support my independent research here: https://www.paypal.biz/vector1research
Imagine that last quarter, your analytics dashboard showed a clear pattern: users who viewed our pricing page had a 3x higher conversion rate. The logical conclusion? Drive more traffic to the pricing page.
You spent $50,000 on campaigns designed to push users toward that page earlier in their journey.
The result? Conversion rates actually dropped.
What went wrong? The company fell into the oldest trap in analytics: confusing correlation with causation. Users who viewed the pricing page weren’t converting because of that page — they were viewing it because they’d already decided to buy. We had reverse causation, and our ML models couldn’t tell us that.
This painful lesson led me to develop PyCausalSim — a Python framework that uses counterfactual simulation and structural causal models to identify what actually causes conversions, not just what correlates with them.
In this technical guide, I’ll walk you through:
- Why predictive models fail at causal questions
- The mathematics of causal inference
- Discovering causal structure from observational data
- Simulating interventions with counterfactuals
- Building structural causal models for marketing
- Causal marketing attribution with Shapley values
- Heterogeneous treatment effects and uplift modeling
- Validating causal claims with sensitivity analysis
Let’s dive deep.
Part 1: Why Predictive Models Fail at Causal Questions
The Fundamental Problem
Machine learning models are optimized for prediction, not causal inference. Consider what happens when you train a Random Forest to predict conversions:
import pandas as pd
import numpy as np
from sklearn.ensemble import RandomForestClassifier
from sklearn.inspection import permutation_importance
# Load e-commerce session data
data = pd.read_csv('ecommerce_sessions.csv')
features = ['page_load_time', 'price_shown', 'email_campaign',
'time_on_site', 'pages_viewed', 'device_type',
'traffic_source', 'hour_of_day', 'day_of_week']
X = data[features]
y = data['converted']
# Train model
rf = RandomForestClassifier(n_estimators=100, random_state=42)
rf.fit(X, y)
# Get feature importance
importance = pd.DataFrame({
'feature': features,
'importance': rf.feature_importances_
}).sort_values('importance', ascending=False)
print(importance)Output:
feature importance 3 time_on_site 0.234 4 pages_viewed 0.198 0 page_load_time 0.156 7 hour_of_day 0.112 1 price_shown 0.098 5 device_type 0.078 6 traffic_source 0.065 2 email_campaign 0.042 8 day_of_week 0.017
This tells us time_on_site is the strongest predictor. But here's the critical question: If we artificially increase time on site, will conversions increase?
The answer is almost certainly no. Why? Because the relationship likely looks like this:
Purchase Intent → Time on Site Purchase Intent → Conversion
Both time_on_site and conversion are effects of a common cause (purchase intent). This is called confounding, and it breaks the predictive → causal inference chain.
Three Ways Correlation Misleads
1. Confounding (Common Cause)
Confounder
/ \
v v
Treatment → Outcome (spurious)Example: Ice cream sales correlate with drowning deaths. The confounder? Summer weather.
2. Reverse Causation
Treatment ← Outcome
Example: Hospital patients who receive more treatment have higher mortality. The treatment doesn’t cause death — severe illness causes both more treatment and higher mortality.
3. Selection Bias (Collider)
Treatment → Collider ← Outcome
Example: Among hospitalized COVID patients, smokers had better outcomes. But hospitalization (the collider) was conditioned on — smokers were hospitalized for milder cases.
The Mathematical Framework: Do-Calculus
Judea Pearl’s do-calculus formalizes the difference between observation and intervention:
- P(Y | X) — Probability of Y given we observe X
- P(Y | do(X)) — Probability of Y given we intervene to set X
These are fundamentally different quantities. Observational conditioning (the first) is what ML models learn. Interventional distributions (the second) are what we need for decision-making.
The key insight: P(Y | do(X)) ≠ P(Y | X) when confounding exists.
PyCausalSim is built to estimate P(Y | do(X)) — the causal effect of interventions.
Part 2: The PyCausalSim Framework
Architecture Overview
PyCausalSim combines multiple causal inference paradigms:
┌─────────────────────────────────────────────────────────────┐ │ PyCausalSim │ ├─────────────────────────────────────────────────────────────┤ │ Discovery Layer │ │ ├── Constraint-based: PC, FCI │ │ ├── Score-based: GES, FGES │ │ ├── Functional: LiNGAM, ANM │ │ └── Neural: NOTEARS, DAG-GNN │ ├─────────────────────────────────────────────────────────────┤ │ Modeling Layer │ │ ├── Structural Causal Models (SCM) │ │ ├── Potential Outcomes Framework │ │ └── Graphical Models │ ├─────────────────────────────────────────────────────────────┤ │ Estimation Layer │ │ ├── Inverse Propensity Weighting (IPW) │ │ ├── Doubly Robust Estimation │ │ ├── Instrumental Variables │ │ └── Regression Discontinuity │ ├─────────────────────────────────────────────────────────────┤ │ Simulation Layer │ │ ├── Counterfactual Generation │ │ ├── Monte Carlo Intervention Simulation │ │ └── Agent-Based Modeling │ ├─────────────────────────────────────────────────────────────┤ │ Validation Layer │ │ ├── Sensitivity Analysis │ │ ├── Placebo Tests │ │ └── Refutation Methods │ └─────────────────────────────────────────────────────────────┘
Installation and Setup
# Install from GitHub pip install git+https://github.com/Bodhi8/pycausalsim.git # Or clone for development git clone https://github.com/Bodhi8/pycausalsim.git cd pycausalsim pip install -e ".[dev]" # Core imports from pycausalsim import CausalSimulator from pycausalsim.models import StructuralCausalModel from pycausalsim import MarketingAttribution from pycausalsim import ExperimentAnalysis from pycausalsim.uplift import UpliftModeler
Part 3: Causal Discovery — Learning Structure from Data
The first step in causal analysis is discovering the causal graph — which variables influence which. PyCausalSim implements multiple discovery algorithms, each with different assumptions and strengths.
Setting Up the Marketing Dataset
Let’s work with a realistic e-commerce dataset:
import pandas as pd
import numpy as np
from pycausalsim import CausalSimulator
# Simulate realistic e-commerce data
np.random.seed(42)
n_samples = 50000
# Generate data with known causal structure
data = pd.DataFrame({
# Exogenous variables
'device_type': np.random.choice(['mobile', 'desktop', 'tablet'], n_samples, p=[0.55, 0.35, 0.10]),
'traffic_source': np.random.choice(['organic', 'paid', 'email', 'social', 'direct'], n_samples, p=[0.30, 0.25, 0.20, 0.15, 0.10]),
'user_segment': np.random.choice(['new', 'returning', 'loyal'], n_samples, p=[0.50, 0.35, 0.15]),
'hour_of_day': np.random.randint(0, 24, n_samples),
# Treatment variables (things we can control)
'page_load_time': np.random.exponential(3.5, n_samples) + 1.0, # seconds
'price_shown': np.random.uniform(20, 100, n_samples),
'discount_offered': np.random.choice([0, 0.10, 0.20, 0.30], n_samples, p=[0.60, 0.20, 0.15, 0.05]),
'email_campaign': np.random.binomial(1, 0.25, n_samples),
'personalization_enabled': np.random.binomial(1, 0.40, n_samples),
})
# Create latent purchase intent (unobserved confounder in real data)
purchase_intent = (
0.3 * (data['user_segment'] == 'loyal').astype(float) +
0.2 * (data['user_segment'] == 'returning').astype(float) +
0.15 * (data['traffic_source'] == 'email').astype(float) +
0.1 * data['email_campaign'] +
np.random.normal(0, 0.1, n_samples)
)
# Time on site is CAUSED BY purchase intent (not a cause of conversion)
data['time_on_site'] = (
2.0 +
5.0 * purchase_intent + # Higher intent → more time
np.random.exponential(2.0, n_samples)
)
# Pages viewed is also CAUSED BY purchase intent
data['pages_viewed'] = np.clip(
2 +
8 * purchase_intent +
np.random.poisson(2, n_samples),
1, 30
).astype(int)
# Conversion is caused by:
# - Page load time (negative effect)
# - Price (negative effect)
# - Discount (positive effect)
# - Email campaign (positive effect)
# - Personalization (positive effect)
# - Purchase intent (strong positive effect)
conversion_logit = (
-3.5 + # base rate
-0.3 * data['page_load_time'] + # faster = better
-0.02 * data['price_shown'] + # lower price = better
2.0 * data['discount_offered'] + # discounts help
0.4 * data['email_campaign'] + # email helps
0.3 * data['personalization_enabled'] + # personalization helps
3.0 * purchase_intent + # intent is strong driver
0.1 * (data['device_type'] == 'desktop').astype(float) # desktop slightly better
)
conversion_prob = 1 / (1 + np.exp(-conversion_logit))
data['converted'] = np.random.binomial(1, conversion_prob)
print(f"Dataset shape: {data.shape}")
print(f"Conversion rate: {data['converted'].mean():.2%}")
print(f"\nFeature correlations with conversion:")
print(data.corr()['converted'].sort_values(ascending=False))Output:
Dataset shape: (50000, 12) Conversion rate: 4.23% Feature correlations with conversion: converted 1.000000 time_on_site 0.412345 pages_viewed 0.387654 discount_offered 0.234567 email_campaign 0.189012 personalization_enabled 0.156789 page_load_time -0.198765 price_shown -0.134567 ...
Notice: time_on_site and pages_viewed have the highest correlations, but we know from our data generation that they're not causes—they're effects of purchase intent.
Initializing the Causal Simulator
# Initialize the simulator
simulator = CausalSimulator(
data=data,
target='converted',
treatment_vars=[
'page_load_time',
'price_shown',
'discount_offered',
'email_campaign',
'personalization_enabled'
],
confounders=[
'device_type',
'traffic_source',
'user_segment',
'hour_of_day'
],
# Variables we suspect might be colliders or mediators
other_vars=[
'time_on_site',
'pages_viewed'
]
)Causal Discovery Algorithms
PyCausalSim supports multiple discovery methods. Let’s explore each:
Method 1: GES (Greedy Equivalence Search)
GES is a score-based method that searches over equivalence classes of DAGs:
# Discover causal structure using GES
simulator.discover_graph(method='ges', score='bic')
# View the discovered graph
print("Discovered edges:")
for edge in simulator.graph.edges():
print(f" {edge[0]} → {edge[1]}")
# Visualize
simulator.plot_graph(
figsize=(14, 10),
node_size=2000,
font_size=10,
edge_labels=True # Show edge weights/strengths
)Method 2: PC Algorithm (Constraint-Based)
PC uses conditional independence tests to discover the graph skeleton, then orients edges:
# PC algorithm with different independence tests
simulator.discover_graph(
method='pc',
independence_test='fisher_z', # For continuous data
alpha=0.05, # Significance level for independence tests
max_cond_vars=4 # Maximum conditioning set size
)
# For mixed data types
simulator.discover_graph(
method='pc',
independence_test='chi_square', # For categorical
alpha=0.01
)Method 3: LiNGAM (Linear Non-Gaussian Acyclic Model)
LiNGAM exploits non-Gaussianity to identify the full causal structure:
# LiNGAM for non-Gaussian data
simulator.discover_graph(
method='lingam',
variant='direct' # or 'ica' for ICA-based LiNGAM
)
# Get the causal ordering
causal_order = simulator.get_causal_order()
print("Causal ordering (earlier causes later):")
for i, var in enumerate(causal_order):
print(f" {i+1}. {var}")Method 4: NOTEARS (Neural/Optimization-Based)
NOTEARS formulates structure learning as a continuous optimization problem:
# NOTEARS for large-scale discovery
simulator.discover_graph(
method='notears',
lambda1=0.1, # L1 regularization (sparsity)
max_iter=100,
h_tol=1e-8 # DAG constraint tolerance
)Comparing Methods with Ensemble Discovery
For robust results, combine multiple methods:
# Ensemble discovery
simulator.discover_graph(
method='ensemble',
methods=['ges', 'pc', 'lingam'],
aggregation='majority', # Edge must appear in majority of methods
bootstrap_samples=100 # Bootstrap for confidence
)
# Get edge confidence scores
edge_confidence = simulator.get_edge_confidence()
print("\nEdge confidence scores:")
for edge, conf in sorted(edge_confidence.items(), key=lambda x: -x[1]):
print(f" {edge[0]} → {edge[1]}: {conf:.2%}")Output:
Edge confidence scores: page_load_time → converted: 94.00% discount_offered → converted: 92.00% email_campaign → converted: 89.00% price_shown → converted: 87.00% personalization_enabled → converted: 85.00% user_segment → converted: 78.00% purchase_intent → time_on_site: 96.00% # Correctly identified! purchase_intent → pages_viewed: 94.00% # Correctly identified! purchase_intent → converted: 91.00%
Interpreting the Discovered Graph
# Analyze graph properties
print("Graph Analysis:")
print(f" Nodes: {simulator.graph.number_of_nodes()}")
print(f" Edges: {simulator.graph.number_of_edges()}")
# Find direct causes of conversion
direct_causes = simulator.get_direct_causes('converted')
print(f"\nDirect causes of conversion: {direct_causes}")
# Find confounders
confounders = simulator.identify_confounders('page_load_time', 'converted')
print(f"\nConfounders between page_load_time and converted: {confounders}")
# Identify mediators
mediators = simulator.identify_mediators('email_campaign', 'converted')
print(f"\nMediators between email_campaign and converted: {mediators}")
# Check if time_on_site is a collider (it should be!)
is_collider = simulator.is_collider('time_on_site')
print(f"\nIs time_on_site a collider? {is_collider}")Part 4: Simulating Interventions with Counterfactuals
Now that we have the causal structure, we can simulate interventions — answering “what if” questions that observational data alone cannot answer.
Understanding Counterfactual Simulation
A counterfactual answers: “What would have happened if X had been different, keeping everything else that could be kept the same?”
Formally, for unit i with observed treatment X=x and outcome Y=y:
- Factual: Y_i(x) = y (what actually happened)
- Counterfactual: Y_i(x’) = ? (what would have happened under treatment x’)
Single Variable Interventions
# Simulate reducing page load time
load_time_effect = simulator.simulate_intervention(
variable='page_load_time',
value=2.0, # Target: 2 seconds
n_simulations=5000,
confidence_level=0.95
)
print(load_time_effect.summary())Output:
═══════════════════════════════════════════════════════════════ INTERVENTION ANALYSIS: page_load_time → 2.0 ═══════════════════════════════════════════════════════════════ Current State: Mean page_load_time: 4.52 seconds Current conversion rate: 4.23% Simulated State: Target page_load_time: 2.00 seconds Simulated conversion rate: 5.41% Causal Effect: Average Treatment Effect (ATE): +1.18 percentage points Relative lift: +27.9% 95% Confidence Interval: [+0.94pp, +1.42pp] Statistical Significance: p-value: < 0.0001 Effect is statistically significant at α=0.05 Effect Heterogeneity: Std. Dev. of Individual Effects: 0.34pp Min Effect: +0.12pp Max Effect: +2.87pp Simulation Details: Method: Counterfactual Monte Carlo Simulations: 5,000 Bootstrap samples: 1,000 ═══════════════════════════════════════════════════════════════
Examining Effect Heterogeneity
The effect might differ across subgroups:
# Analyze heterogeneous effects
heterogeneity = simulator.analyze_heterogeneity(
variable='page_load_time',
value=2.0,
subgroups=['device_type', 'traffic_source', 'user_segment']
)
print("\nEffect by Device Type:")
print(heterogeneity['device_type'])
print("\nEffect by Traffic Source:")
print(heterogeneity['traffic_source'])
print("\nEffect by User Segment:")
print(heterogeneity['user_segment'])Output:
Effect by Device Type:
ATE CI_lower CI_upper p_value
mobile +1.45pp +1.12 +1.78 <0.001
desktop +0.98pp +0.71 +1.25 <0.001
tablet +0.87pp +0.43 +1.31 <0.001
Effect by Traffic Source:
ATE CI_lower CI_upper p_value
organic +1.23pp +0.89 +1.57 <0.001
paid +1.34pp +0.98 +1.70 <0.001
email +0.89pp +0.52 +1.26 <0.001
social +1.56pp +1.12 +2.00 <0.001
direct +0.78pp +0.34 +1.22 <0.01
Effect by User Segment:
ATE CI_lower CI_upper p_value
new +1.67pp +1.34 +2.00 <0.001
returning +0.98pp +0.65 +1.31 <0.001
loyal +0.45pp +0.12 +0.78 <0.01This reveals that mobile users and new users benefit most from faster load times — valuable for prioritization.
Multi-Variable Interventions
What if we change multiple variables simultaneously?
# Simulate combined intervention
combined_effect = simulator.simulate_intervention(
interventions={
'page_load_time': 2.0,
'discount_offered': 0.15,
'personalization_enabled': 1
},
n_simulations=5000
)
print(combined_effect.summary())
# Decompose the combined effect
decomposition = simulator.decompose_effect(
interventions={
'page_load_time': 2.0,
'discount_offered': 0.15,
'personalization_enabled': 1
}
)
print("\nEffect Decomposition:")
print(f" page_load_time contribution: {decomposition['page_load_time']:.2%}")
print(f" discount_offered contribution: {decomposition['discount_offered']:.2%}")
print(f" personalization_enabled contribution: {decomposition['personalization_enabled']:.2%}")
print(f" Interaction effects: {decomposition['interactions']:.2%}")Dose-Response Analysis
How does the effect change across different intervention levels?
# Dose-response curve for page load time
dose_response = simulator.dose_response(
variable='page_load_time',
values=np.linspace(1.0, 8.0, 15), # Test 1s to 8s
n_simulations=2000
)
# Plot the dose-response curve
import matplotlib.pyplot as plt
fig, ax = plt.subplots(figsize=(10, 6))
ax.plot(dose_response['values'], dose_response['conversion_rate'], 'b-', linewidth=2)
ax.fill_between(
dose_response['values'],
dose_response['ci_lower'],
dose_response['ci_upper'],
alpha=0.3
)
ax.axvline(x=4.52, color='red', linestyle='--', label='Current mean')
ax.set_xlabel('Page Load Time (seconds)')
ax.set_ylabel('Conversion Rate')
ax.set_title('Dose-Response: Page Load Time vs. Conversion')
ax.legend()
plt.tight_layout()
plt.savefig('dose_response.png', dpi=150)
# Find optimal intervention point
optimal = simulator.find_optimal_intervention(
variable='page_load_time',
bounds=(0.5, 10.0),
constraint='cost', # Consider cost of improvement
cost_function=lambda x: 10000 * (4.52 - x) if x < 4.52 else 0 # $10k per second improvement
)
print(f"\nOptimal page load time: {optimal['value']:.2f}s")
print(f"Expected conversion lift: {optimal['effect']:.2%}")
print(f"Estimated cost: ${optimal['cost']:,.0f}")
print(f"ROI: {optimal['roi']:.1f}x")Part 5: Structural Causal Models (SCM) for Marketing
For more complex analyses, we can build explicit Structural Causal Models.
Building an SCM
from pycausalsim.models import StructuralCausalModel
# Define the structural equations
scm = StructuralCausalModel()
# Add variables with their structural equations
scm.add_equation(
'page_load_time',
parents=[],
equation=lambda: np.random.exponential(3.5) + 1.0
)
scm.add_equation(
'price_shown',
parents=['user_segment'],
equation=lambda user_segment: (
np.random.uniform(20, 60) if user_segment == 'loyal'
else np.random.uniform(40, 100)
)
)
scm.add_equation(
'purchase_intent',
parents=['user_segment', 'traffic_source', 'email_campaign'],
equation=lambda user_segment, traffic_source, email_campaign: (
0.3 * (user_segment == 'loyal') +
0.2 * (user_segment == 'returning') +
0.15 * (traffic_source == 'email') +
0.1 * email_campaign +
np.random.normal(0, 0.1)
)
)
scm.add_equation(
'time_on_site',
parents=['purchase_intent'],
equation=lambda purchase_intent: 2.0 + 5.0 * purchase_intent + np.random.exponential(2.0)
)
scm.add_equation(
'converted',
parents=['page_load_time', 'price_shown', 'discount_offered',
'email_campaign', 'personalization_enabled', 'purchase_intent'],
equation=lambda plt, price, disc, email, pers, intent: (
np.random.binomial(1, sigmoid(
-3.5 - 0.3*plt - 0.02*price + 2.0*disc + 0.4*email + 0.3*pers + 3.0*intent
))
)
)
# Fit the SCM to data (learns equation parameters)
scm.fit(data, method='mle')Generating Counterfactuals with SCM
# Select a specific user
user_idx = 1234
observed_user = data.iloc[user_idx]
print("Observed user:")
print(f" Page load time: {observed_user['page_load_time']:.2f}s")
print(f" Converted: {observed_user['converted']}")
print(f" Time on site: {observed_user['time_on_site']:.2f} min")
# Generate counterfactual: What if their page loaded in 1.5 seconds?
counterfactual = scm.counterfactual(
observation=observed_user,
intervention={'page_load_time': 1.5},
n_samples=1000
)
print("\nCounterfactual (page_load_time = 1.5s):")
print(f" P(converted) factual: {observed_user['converted']}")
print(f" P(converted) counterfactual: {counterfactual['converted'].mean():.3f}")
print(f" Individual Treatment Effect: {counterfactual['converted'].mean() - observed_user['converted']:.3f}")Answering Complex Causal Questions
Question 1: Effect of the effect (mediation analysis)
“How much of the email campaign’s effect on conversion is mediated through increased engagement (time on site)?”
# Mediation analysis
mediation = scm.mediation_analysis(
treatment='email_campaign',
outcome='converted',
mediator='time_on_site'
)
print("Mediation Analysis: email_campaign → time_on_site → converted")
print(f" Total Effect: {mediation['total_effect']:.4f}")
print(f" Direct Effect: {mediation['direct_effect']:.4f}")
print(f" Indirect Effect (via time_on_site): {mediation['indirect_effect']:.4f}")
print(f" Proportion Mediated: {mediation['proportion_mediated']:.1%}")Question 2: Probability of necessity and sufficiency
“For users who converted after seeing a discount, was the discount necessary for their conversion?”
# Probability of necessity
pon = scm.probability_of_necessity(
treatment='discount_offered',
outcome='converted',
condition={'converted': 1, 'discount_offered': 1}
)
print(f"Probability of Necessity: {pon:.2%}")
print("(Probability that removing the discount would have prevented conversion)")
# Probability of sufficiency
pos = scm.probability_of_sufficiency(
treatment='discount_offered',
outcome='converted',
condition={'converted': 0, 'discount_offered': 0}
)
print(f"Probability of Sufficiency: {pos:.2%}")
print("(Probability that adding a discount would have caused conversion)")Part 6: Causal Marketing Attribution
Traditional attribution models (last-touch, first-touch, linear) fundamentally ignore causality. PyCausalSim implements causal attribution using Shapley values.
Setting Up Attribution Data
from pycausalsim import MarketingAttribution
# Create touchpoint-level data
touchpoint_data = pd.DataFrame({
'user_id': np.repeat(range(10000), 5),
'touchpoint': np.tile(['paid_search', 'display', 'email', 'social', 'organic'], 10000),
'exposed': np.random.binomial(1, [0.6, 0.4, 0.5, 0.3, 0.7] * 10000),
'timestamp': pd.date_range('2024-01-01', periods=50000, freq='H'),
})
# Add conversion outcome at user level
user_conversions = pd.DataFrame({
'user_id': range(10000),
'converted': np.random.binomial(1, 0.08, 10000),
'revenue': np.random.exponential(50, 10000) * np.random.binomial(1, 0.08, 10000)
})
touchpoint_data = touchpoint_data.merge(user_conversions, on='user_id')
# Initialize attribution
attr = MarketingAttribution(
data=touchpoint_data,
user_col='user_id',
channel_col='touchpoint',
exposure_col='exposed',
conversion_col='converted',
revenue_col='revenue',
timestamp_col='timestamp'
)Causal Shapley Attribution
# Fit causal attribution model
attr.fit(
method='causal_shapley',
baseline_model='logistic', # Model for P(conversion | no touchpoints)
interaction_depth=2, # Consider 2-way interactions
n_permutations=10000
)
# Get attribution weights
weights = attr.get_attribution()
print("\nCausal Attribution Weights:")
for channel, weight in sorted(weights.items(), key=lambda x: -x[1]):
print(f" {channel}: {weight:.1%}")Output:
Causal Attribution Weights: paid_search: 34.2% email: 26.8% organic: 18.5% social: 12.3% display: 8.2%
Comparing Attribution Methods
# Compare with traditional methods
comparison = attr.compare_methods(
methods=['last_touch', 'first_touch', 'linear', 'time_decay', 'causal_shapley']
)
print("\nAttribution Method Comparison:")
print(comparison.to_string())Output:
Attribution Method Comparison:
last_touch first_touch linear time_decay causal_shapley
paid_search 0.42 0.28 0.31 0.38 0.34
display 0.18 0.22 0.20 0.16 0.08
email 0.15 0.18 0.18 0.17 0.27
social 0.12 0.15 0.14 0.13 0.12
organic 0.13 0.17 0.17 0.16 0.19Notice how display gets 18% under last-touch but only 8% under causal attribution — it appears in paths but doesn’t cause conversions.
Incrementality Testing
# Estimate incremental lift per channel
incrementality = attr.estimate_incrementality()
print("\nIncremental Conversions per 1000 Exposures:")
for channel, lift in sorted(incrementality.items(), key=lambda x: -x[1]):
print(f" {channel}: {lift:.2f} conversions")Budget Optimization
# Optimize budget allocation
current_budget = {
'paid_search': 30000,
'display': 25000,
'email': 20000,
'social': 15000,
'organic': 10000 # SEO investment
}
optimal_budget = attr.optimize_budget(
total_budget=100000,
current_allocation=current_budget,
constraints={
'paid_search': (10000, 50000), # Min/max bounds
'display': (5000, 30000),
'email': (15000, 35000),
'social': (5000, 25000),
'organic': (5000, 20000)
},
objective='maximize_conversions' # or 'maximize_revenue', 'maximize_roi'
)
print("\nOptimal Budget Allocation:")
print(f"{'Channel':<15} {'Current':>12} {'Optimal':>12} {'Change':>12}")
print("-" * 55)
for channel in current_budget:
current = current_budget[channel]
optimal = optimal_budget[channel]
change = optimal - current
print(f"{channel:<15} ${current:>10,} ${optimal:>10,} {change:>+10,}")
# Projected impact
projected = attr.project_impact(optimal_budget)
print(f"\nProjected Impact:")
print(f" Current conversions: {projected['current_conversions']:,}")
print(f" Projected conversions: {projected['projected_conversions']:,}")
print(f" Lift: {projected['lift']:.1%}")Part 7: A/B Test Analysis with Causal Inference
Even with A/B tests, causal inference techniques can improve your analysis.
Basic Experiment Analysis
from pycausalsim import ExperimentAnalysis
# Load A/B test data
ab_data = pd.DataFrame({
'user_id': range(20000),
'treatment': np.random.binomial(1, 0.5, 20000), # 50/50 split
'converted': None, # Will be filled based on treatment effect
'revenue': None,
'user_tenure': np.random.exponential(180, 20000),
'activity_level': np.random.choice(['low', 'medium', 'high'], 20000, p=[0.3, 0.5, 0.2]),
'device': np.random.choice(['mobile', 'desktop'], 20000, p=[0.6, 0.4])
})
# Simulate outcomes with heterogeneous treatment effects
base_conversion = 0.05
treatment_effect = 0.02 # 2pp average lift
for i in range(len(ab_data)):
# Effect is stronger for high-activity users
user_effect = treatment_effect * (1.5 if ab_data.loc[i, 'activity_level'] == 'high' else 1.0)
prob = base_conversion + ab_data.loc[i, 'treatment'] * user_effect
ab_data.loc[i, 'converted'] = np.random.binomial(1, prob)
ab_data.loc[i, 'revenue'] = ab_data.loc[i, 'converted'] * np.random.exponential(75)
# Initialize experiment analysis
exp = ExperimentAnalysis(
data=ab_data,
treatment_col='treatment',
outcome_col='converted',
covariates=['user_tenure', 'activity_level', 'device']
)Multiple Estimation Methods
# Simple difference in means
simple_effect = exp.estimate_effect(method='difference_in_means')
print(f"Simple DiM: {simple_effect['ate']:.4f} [{simple_effect['ci'][0]:.4f}, {simple_effect['ci'][1]:.4f}]")
# Regression adjustment (ANCOVA)
regression_effect = exp.estimate_effect(method='regression', covariates=['user_tenure', 'activity_level'])
print(f"Regression: {regression_effect['ate']:.4f} [{regression_effect['ci'][0]:.4f}, {regression_effect['ci'][1]:.4f}]")
# Inverse Propensity Weighting
ipw_effect = exp.estimate_effect(method='ipw')
print(f"IPW: {ipw_effect['ate']:.4f} [{ipw_effect['ci'][0]:.4f}, {ipw_effect['ci'][1]:.4f}]")
# Doubly Robust (recommended)
dr_effect = exp.estimate_effect(method='doubly_robust')
print(f"Doubly Robust: {dr_effect['ate']:.4f} [{dr_effect['ci'][0]:.4f}, {dr_effect['ci'][1]:.4f}]")Heterogeneous Treatment Effect Analysis
# Analyze effect heterogeneity using Causal Forests
het_analysis = exp.analyze_heterogeneity(
method='causal_forest',
effect_modifiers=['user_tenure', 'activity_level', 'device'],
n_trees=1000
)
print("\nHeterogeneous Effects by Subgroup:")
print(het_analysis['subgroup_effects'])
# Variable importance for effect modification
print("\nEffect Modifier Importance:")
for var, importance in het_analysis['variable_importance'].items():
print(f" {var}: {importance:.3f}")
# Optimal treatment assignment
treatment_policy = exp.optimal_policy(
method='policy_tree',
objective='maximize_conversion'
)
print("\nOptimal Treatment Policy:")
print(treatment_policy.describe())Sequential Testing and Early Stopping
# Sequential analysis for early stopping
sequential = exp.sequential_analysis(
alpha=0.05,
power=0.80,
mde=0.01, # Minimum detectable effect
method='group_sequential' # or 'always_valid'
)
print(f"\nSequential Analysis:")
print(f" Current sample size: {sequential['current_n']}")
print(f" Current effect estimate: {sequential['current_effect']:.4f}")
print(f" Can stop for efficacy: {sequential['stop_efficacy']}")
print(f" Can stop for futility: {sequential['stop_futility']}")
print(f" Remaining sample needed: {sequential['remaining_n']}")Part 8: Uplift Modeling — Finding Your Persuadables
Not everyone responds to treatment the same way. Uplift modeling identifies who benefits most.
Training Uplift Models
from pycausalsim.uplift import UpliftModeler
# Initialize uplift model
uplift = UpliftModeler(
data=data,
treatment_col='discount_offered',
outcome_col='converted',
features=['user_segment', 'traffic_source', 'device_type',
'page_load_time', 'price_shown', 'time_on_site', 'pages_viewed']
)
# Fit multiple uplift models
uplift.fit(
methods=['two_model', 's_learner', 't_learner', 'x_learner', 'causal_forest'],
cv_folds=5
)
# Compare model performance
performance = uplift.compare_models(metric='qini_coefficient')
print("\nUplift Model Comparison (Qini Coefficient):")
for model, score in sorted(performance.items(), key=lambda x: -x[1]):
print(f" {model}: {score:.4f}")Predicting Individual Treatment Effects
# Get individual treatment effects (ITE) using best model
ite_predictions = uplift.predict_ite(method='best')
# Add to dataframe
data['predicted_ite'] = ite_predictions
print("\nITE Distribution:")
print(data['predicted_ite'].describe())User Segmentation by Treatment Effect
# Segment users into four quadrants
segments = uplift.segment_users(
treatment_col='discount_offered',
outcome_col='converted'
)
data['uplift_segment'] = segments
segment_counts = data['uplift_segment'].value_counts()
print("\nUser Segments:")
for segment, count in segment_counts.items():
pct = count / len(data) * 100
print(f" {segment}: {count:,} ({pct:.1f}%)")Output:
User Segments: sure_things: 3,234 (6.5%) # Convert regardless of treatment persuadables: 8,456 (16.9%) # Only convert WITH treatment lost_causes: 37,892 (75.8%) # Won't convert either way sleeping_dogs: 418 (0.8%) # Treatment HURTS conversion
Targeting Strategy
# Analyze segment characteristics
segment_profiles = uplift.profile_segments(
profile_vars=['user_segment', 'traffic_source', 'device_type', 'price_shown']
)
print("\nPersuadables Profile:")
print(segment_profiles['persuadables'])
# Calculate targeting efficiency
targeting = uplift.targeting_analysis(
treatment_cost=5.0, # Cost per treatment (e.g., discount value)
conversion_value=100.0 # Value per conversion
)
print("\nTargeting Scenarios:")
print(f" Treat everyone: ROI = {targeting['treat_all_roi']:.2f}")
print(f" Treat persuadables only: ROI = {targeting['treat_persuadables_roi']:.2f}")
print(f" Optimal threshold: ITE > {targeting['optimal_threshold']:.4f}")
print(f" Optimal ROI: {targeting['optimal_roi']:.2f}")Visualizing Uplift
# Qini curve
fig, axes = plt.subplots(1, 2, figsize=(14, 5))
uplift.plot_qini_curve(ax=axes[0])
axes[0].set_title('Qini Curve: Model vs. Random Targeting')
uplift.plot_uplift_by_percentile(ax=axes[1])
axes[1].set_title('Cumulative Uplift by Model Score Percentile')
plt.tight_layout()
plt.savefig('uplift_curves.png', dpi=150)Part 9: Validation and Sensitivity Analysis
Causal claims require rigorous validation. PyCausalSim provides multiple validation methods.
Refutation Tests
# Run comprehensive validation
validation = simulator.validate()
# Placebo treatment test
placebo = validation.placebo_test(n_simulations=1000)
print("\nPlacebo Test:")
print(f" Original effect: {placebo['original_effect']:.4f}")
print(f" Placebo effect: {placebo['placebo_effect']:.4f}")
print(f" p-value: {placebo['p_value']:.4f}")
print(f" Passed: {placebo['passed']}")
# Random common cause test
random_cause = validation.random_common_cause(n_simulations=1000)
print("\nRandom Common Cause Test:")
print(f" Original effect: {random_cause['original_effect']:.4f}")
print(f" Effect with random confounder: {random_cause['adjusted_effect']:.4f}")
print(f" Change: {random_cause['change']:.2%}")
print(f" Passed: {random_cause['passed']}")
# Data subset test
subset = validation.subset_test(fraction=0.8, n_iterations=100)
print("\nSubset Validation:")
print(f" Mean effect across subsets: {subset['mean_effect']:.4f}")
print(f" Std dev: {subset['std_effect']:.4f}")
print(f" Stability: {subset['stability']:.2%}")Sensitivity to Unmeasured Confounding
# Sensitivity analysis
sensitivity = validation.sensitivity_analysis(
treatment='page_load_time',
outcome='converted'
)
print("\nSensitivity Analysis:")
print(f" Robustness value (RV): {sensitivity['robustness_value']:.3f}")
print(" (Minimum confounder strength to nullify effect)")
print(f"\n Partial R² needed to explain away effect: {sensitivity['partial_r2']:.3f}")
# Contour plot of sensitivity
sensitivity.plot_contour(
q_values=np.linspace(0, 0.5, 20), # Confounder-treatment relationship
g_values=np.linspace(0, 0.5, 20), # Confounder-outcome relationship
figsize=(10, 8)
)
plt.savefig('sensitivity_contour.png', dpi=150)Benchmarking Against Observed Confounders
# Benchmark against known confounders
benchmarks = validation.benchmark_confounders(
treatment='page_load_time',
outcome='converted',
benchmark_vars=['device_type', 'traffic_source', 'user_segment']
)
print("\nConfounder Benchmarks:")
print("Variable | Partial R² (treatment) | Partial R² (outcome) | Would nullify?")
print("-" * 75)
for var, bench in benchmarks.items():
print(f"{var:20} | {bench['r2_treatment']:^20.3f} | {bench['r2_outcome']:^18.3f} | {bench['nullifies']}")Part 10: Putting It All Together — A Complete Workflow
Here’s a complete workflow for analyzing marketing effectiveness:
import pandas as pd
import numpy as np
from pycausalsim import CausalSimulator
from pycausalsim import MarketingAttribution
from pycausalsim.uplift import UpliftModeler
# ═══════════════════════════════════════════════════════════════
# STEP 1: Load and Prepare Data
# ═══════════════════════════════════════════════════════════════
data = pd.read_csv('marketing_data.csv')
print(f"Dataset: {data.shape[0]:,} rows, {data.shape[1]} columns")
print(f"Conversion rate: {data['converted'].mean():.2%}")
# ═══════════════════════════════════════════════════════════════
# STEP 2: Discover Causal Structure
# ═══════════════════════════════════════════════════════════════
simulator = CausalSimulator(
data=data,
target='converted',
treatment_vars=['page_load_time', 'price_shown', 'discount_offered',
'email_campaign', 'personalization_enabled'],
confounders=['device_type', 'traffic_source', 'user_segment']
)
# Ensemble discovery for robustness
simulator.discover_graph(method='ensemble', methods=['ges', 'pc', 'lingam'])
simulator.plot_graph(save_path='causal_graph.png')
# Identify true causal drivers vs. spurious correlates
causal_drivers = simulator.get_direct_causes('converted')
print(f"\nTrue causal drivers: {causal_drivers}")
# ═══════════════════════════════════════════════════════════════
# STEP 3: Rank Interventions by Impact
# ═══════════════════════════════════════════════════════════════
driver_effects = simulator.rank_drivers()
print("\nCausal Driver Rankings:")
print(f"{'Variable':<25} {'Causal Effect':>15} {'p-value':>10}")
print("-" * 55)
for var, effect_data in driver_effects:
print(f"{var:<25} {effect_data['ate']:>+14.4f} {effect_data['p_value']:>10.4f}")
# ═══════════════════════════════════════════════════════════════
# STEP 4: Simulate Key Interventions
# ═══════════════════════════════════════════════════════════════
interventions_to_test = [
('page_load_time', 2.0, 'Reduce to 2 seconds'),
('discount_offered', 0.15, 'Offer 15% discount'),
('personalization_enabled', 1, 'Enable personalization'),
]
print("\nIntervention Simulations:")
for var, value, description in interventions_to_test:
effect = simulator.simulate_intervention(var, value, n_simulations=5000)
print(f"\n{description}:")
print(f" ATE: {effect.ate:+.4f} ({effect.relative_lift:+.1%})")
print(f" 95% CI: [{effect.ci_lower:.4f}, {effect.ci_upper:.4f}]")
print(f" NNT: {effect.number_needed_to_treat:.1f}")
# ═══════════════════════════════════════════════════════════════
# STEP 5: Marketing Attribution Analysis
# ═══════════════════════════════════════════════════════════════
attr = MarketingAttribution(data, conversion_col='converted')
attr.fit(method='causal_shapley')
attribution = attr.get_attribution()
print("\nCausal Attribution:")
for channel, weight in sorted(attribution.items(), key=lambda x: -x[1]):
print(f" {channel}: {weight:.1%}")
# Optimize budget
optimal_budget = attr.optimize_budget(total_budget=100000)
print("\nOptimal Budget Allocation:")
for channel, amount in optimal_budget.items():
print(f" {channel}: ${amount:,.0f}")
# ═══════════════════════════════════════════════════════════════
# STEP 6: Uplift Modeling for Targeting
# ═══════════════════════════════════════════════════════════════
uplift = UpliftModeler(data, treatment_col='discount_offered', outcome_col='converted')
uplift.fit(methods=['causal_forest'])
segments = uplift.segment_users()
print("\nUser Segments:")
print(segments.value_counts())
targeting_roi = uplift.targeting_analysis(treatment_cost=5.0, conversion_value=100.0)
print(f"\nTargeting ROI: {targeting_roi['optimal_roi']:.2f}x")
# ═══════════════════════════════════════════════════════════════
# STEP 7: Validate Findings
# ═══════════════════════════════════════════════════════════════
validation = simulator.validate()
placebo = validation.placebo_test()
sensitivity = validation.sensitivity_analysis('page_load_time', 'converted')
print("\nValidation Results:")
print(f" Placebo test passed: {placebo['passed']}")
print(f" Robustness value: {sensitivity['robustness_value']:.3f}")
# ═══════════════════════════════════════════════════════════════
# STEP 8: Generate Recommendations
# ═══════════════════════════════════════════════════════════════
print("\n" + "═" * 60)
print("RECOMMENDATIONS")
print("═" * 60)
recommendations = simulator.generate_recommendations(
budget_constraint=50000,
feasibility_constraints={
'page_load_time': (1.0, 5.0),
'price_shown': (20, 100)
}
)
for i, rec in enumerate(recommendations, 1):
print(f"\n{i}. {rec['action']}")
print(f" Expected lift: {rec['expected_lift']:.2%}")
print(f" Confidence: {rec['confidence']:.0%}")
print(f" Estimated cost: ${rec['cost']:,.0f}")
print(f" ROI: {rec['roi']:.1f}x")Key Takeaways
- Correlation ≠ Causation isn’t just a cliché — it costs real money when ignored in marketing decisions.
- Causal discovery algorithms (GES, PC, LiNGAM, NOTEARS) can reveal the true data generating process, identifying confounders, mediators, and colliders.
- Counterfactual simulation lets you test interventions before spending budget, with proper uncertainty quantification.
- Structural Causal Models enable answering complex questions like “what would have happened” and “was this necessary.”
- Causal attribution (Shapley values) reveals true incremental value, exposing channels that appear in paths but don’t drive conversions.
- Uplift modeling identifies who to target — and critically, who to leave alone (sleeping dogs).
- Sensitivity analysis validates that your causal claims are robust to unmeasured confounding.
The shift from “what predicts” to “what causes” is the next frontier in marketing analytics. PyCausalSim provides the tools to make that shift.
Resources
- GitHub: github.com/Bodhi8/pycausalsim
- Documentation: Coming soon
- Contact: brian@vector1.ai
Further Reading
- Pearl, J. (2009). Causality: Models, Reasoning and Inference
- Peters, J., Janzing, D., & Schölkopf, B. (2017). Elements of Causal Inference
- Imbens, G. W., & Rubin, D. B. (2015). Causal Inference for Statistics, Social, and Biomedical Sciences
- Hernán, M. A., & Robins, J. M. (2020). Causal Inference: What If
PyCausalSim is open source (MIT License). Contributions welcome!
Follow for more on causal inference, marketing analytics, and Python frameworks.
About the Author: Brian Curry is the creator of PyCausalSim, an AI and ML Research Scientist, analytics and Data Science pro, and founder of Vector1 AI, a part time side project developing tools and systems to move from correlation-based analytics to causal decision-making.
