xbe.at

Miscellany / Statistics and analysis

Estimating Causal Impact of Late Deliveries on Customer Satisfaction

8 min read, 14 sections

Contents
  1. Understanding Propensity Score Matching Fundamentals
  2. Data Preparation for Late Delivery Analysis
  3. Propensity Score Estimation
  4. Nearest Neighbor Matching Implementation
  5. Implementing Matching and Balance Assessment
  6. Average Treatment Effect Calculation
  7. Implementation of Full PSM Analysis Pipeline
  8. Real-world Example: E-commerce Late Delivery Analysis
  9. Visualization of Matching Quality
  10. Sensitivity Analysis Implementation
  11. Bootstrapping PSM for Uncertainty Estimation
  12. Advanced Matching Strategies
  13. Results and Statistical Inference
  14. Additional Resources

Understanding Propensity Score Matching Fundamentals

Propensity Score Matching (PSM) estimates treatment effects by creating matched pairs of treated and control subjects with similar characteristics. The propensity score represents the probability of receiving treatment based on observed covariates, helping reduce selection bias in observational studies.

import pandas as pd
import numpy as np
from sklearn.linear_model import LogisticRegression
from sklearn.preprocessing import StandardScaler

def calculate_propensity_scores(X, treatment):
    # Fit logistic regression to estimate propensity scores
    model = LogisticRegression(random_state=42)
    model.fit(X, treatment)
    
    # Calculate propensity scores
    propensity_scores = model.predict_proba(X)[:, 1]
    return propensity_scores

Data Preparation for Late Delivery Analysis

In our e-commerce scenario, we'll analyze how late deliveries impact customer satisfaction. We'll prepare a dataset containing customer features, delivery times, and satisfaction scores to demonstrate the practical application of PSM.

# Generate synthetic e-commerce data
np.random.seed(42)
n_samples = 1000

data = {
    'customer_age': np.random.normal(35, 10, n_samples),
    'purchase_amount': np.random.normal(100, 30, n_samples),
    'distance_km': np.random.normal(50, 20, n_samples),
    'is_prime_member': np.random.binomial(1, 0.3, n_samples),
    'late_delivery': np.random.binomial(1, 0.2, n_samples)
}

df = pd.DataFrame(data)
# Simulate satisfaction scores (1-5) with lower scores for late deliveries
df['satisfaction'] = np.where(
    df['late_delivery'] == 1,
    np.random.normal(3, 0.5, n_samples),
    np.random.normal(4, 0.5, n_samples)
).clip(1, 5)

Propensity Score Estimation

The first step in PSM involves estimating propensity scores using logistic regression. We'll use customer characteristics to predict the probability of experiencing a late delivery, creating a balanced comparison between affected and unaffected customers.

# Prepare features for propensity score calculation
features = ['customer_age', 'purchase_amount', 'distance_km', 'is_prime_member']
X = df[features]
treatment = df['late_delivery']

# Standardize features
scaler = StandardScaler()
X_scaled = scaler.fit_transform(X)
X_scaled = pd.DataFrame(X_scaled, columns=features)

# Calculate propensity scores
ps_scores = calculate_propensity_scores(X_scaled, treatment)
df['propensity_score'] = ps_scores

print("Propensity Score Summary:")
print(pd.DataFrame({'Propensity Scores': ps_scores}).describe())

Nearest Neighbor Matching Implementation

This implementation uses the nearest neighbor algorithm to match treated and control units based on propensity scores. We'll create a function that finds the closest match for each treated unit within a specified caliper distance.

def nearest_neighbor_matching(prop_scores, treatment, caliper=0.2):
    treated_indices = np.where(treatment == 1)[0]
    control_indices = np.where(treatment == 0)[0]
    
    # Calculate standard deviation of propensity scores
    prop_scores_std = np.std(prop_scores)
    caliper_distance = caliper * prop_scores_std
    
    matches = []
    used_control = set()
    
    for t_idx in treated_indices:
        treated_score = prop_scores[t_idx]
        
        # Calculate distances to all control units
        distances = np.abs(prop_scores[control_indices] - treated_score)
        
        # Find best match within caliper
        valid_matches = control_indices[distances < caliper_distance]
        valid_distances = distances[distances < caliper_distance]
        
        if len(valid_matches) > 0:
            # Find closest unused control unit
            unused_mask = ~np.isin(valid_matches, list(used_control))
            if any(unused_mask):
                best_idx = valid_matches[unused_mask][np.argmin(valid_distances[unused_mask])]
                matches.append((t_idx, best_idx))
                used_control.add(best_idx)
    
    return matches

Implementing Matching and Balance Assessment

After matching, we need to assess the quality of our matches by comparing the distribution of covariates between treated and control groups before and after matching.

def assess_balance(df, features, matches, treatment_col):
    matched_treated = [m[0] for m in matches]
    matched_control = [m[1] for m in matches]
    
    # Calculate standardized differences before and after matching
    results = []
    
    for feature in features:
        # Before matching
        treated_before = df[df[treatment_col] == 1][feature]
        control_before = df[df[treatment_col] == 0][feature]
        
        # After matching
        treated_after = df.iloc[matched_treated][feature]
        control_after = df.iloc[matched_control][feature]
        
        # Calculate standardized differences
        std_diff_before = (treated_before.mean() - control_before.mean()) / \
                         np.sqrt((treated_before.var() + control_before.var()) / 2)
        
        std_diff_after = (treated_after.mean() - control_after.mean()) / \
                        np.sqrt((treated_after.var() + control_after.var()) / 2)
        
        results.append({
            'Feature': feature,
            'Std_Diff_Before': std_diff_before,
            'Std_Diff_After': std_diff_after
        })
    
    return pd.DataFrame(results)

Average Treatment Effect Calculation

The Average Treatment Effect (ATE) represents the mean difference in outcomes between treated and control groups after matching. We'll implement a function to calculate this effect and its statistical significance.

def calculate_ate(df, matches, outcome_col):
    treated_outcomes = df.iloc[[m[0] for m in matches]][outcome_col]
    control_outcomes = df.iloc[[m[1] for m in matches]][outcome_col]
    
    # Calculate average treatment effect
    ate = np.mean(treated_outcomes - control_outcomes)
    
    # Calculate standard error and confidence interval
    se = np.std(treated_outcomes - control_outcomes) / np.sqrt(len(matches))
    ci_lower = ate - 1.96 * se
    ci_upper = ate + 1.96 * se
    
    # Perform t-test
    from scipy import stats
    t_stat, p_value = stats.ttest_ind(treated_outcomes, control_outcomes)
    
    return {
        'ATE': ate,
        'SE': se,
        'CI_Lower': ci_lower,
        'CI_Upper': ci_upper,
        'P_Value': p_value
    }

Implementation of Full PSM Analysis Pipeline

The complete pipeline combines all previous components into a cohesive analysis workflow. This implementation demonstrates the end-to-end process of propensity score matching for analyzing late delivery impacts on customer satisfaction.

def run_psm_analysis(df, features, treatment_col, outcome_col, caliper=0.2):
    # 1. Prepare data
    X = df[features]
    treatment = df[treatment_col]
    
    # 2. Calculate propensity scores
    scaler = StandardScaler()
    X_scaled = scaler.fit_transform(X)
    ps_scores = calculate_propensity_scores(pd.DataFrame(X_scaled, columns=features), treatment)
    
    # 3. Perform matching
    matches = nearest_neighbor_matching(ps_scores, treatment, caliper)
    
    # 4. Assess balance
    balance_stats = assess_balance(df, features, matches, treatment_col)
    
    # 5. Calculate treatment effect
    treatment_effect = calculate_ate(df, matches, outcome_col)
    
    return {
        'matches': matches,
        'balance_stats': balance_stats,
        'treatment_effect': treatment_effect,
        'propensity_scores': ps_scores
    }

Real-world Example: E-commerce Late Delivery Analysis

We'll analyze the impact of late deliveries on customer satisfaction using our synthetic e-commerce dataset. This example demonstrates the practical application of PSM in a business context.

# Run complete PSM analysis
features = ['customer_age', 'purchase_amount', 'distance_km', 'is_prime_member']
results = run_psm_analysis(
    df=df,
    features=features,
    treatment_col='late_delivery',
    outcome_col='satisfaction',
    caliper=0.2
)

# Display results
print("Number of matched pairs:", len(results['matches']))
print("\nBalance Statistics:")
print(results['balance_stats'])
print("\nTreatment Effect Results:")
print(pd.DataFrame([results['treatment_effect']]))

Visualization of Matching Quality

Understanding the quality of matches is crucial for validating PSM results. We'll create visualizations to compare propensity score distributions and covariate balance before and after matching.

import matplotlib.pyplot as plt
import seaborn as sns

def plot_matching_diagnostics(df, ps_scores, matches, features):
    # Create figure with subplots
    fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(10, 12))
    
    # Plot propensity score distributions
    treated_ps = ps_scores[df['late_delivery'] == 1]
    control_ps = ps_scores[df['late_delivery'] == 0]
    
    sns.kdeplot(data=treated_ps, ax=ax1, label='Treated')
    sns.kdeplot(data=control_ps, ax=ax1, label='Control')
    ax1.set_title('Propensity Score Distributions')
    ax1.set_xlabel('Propensity Score')
    ax1.legend()
    
    # Plot standardized differences
    balance_stats = assess_balance(df, features, matches, 'late_delivery')
    
    balance_stats.plot(
        x='Feature',
        y=['Std_Diff_Before', 'Std_Diff_After'],
        kind='bar',
        ax=ax2
    )
    ax2.set_title('Standardized Differences Before and After Matching')
    ax2.axhline(y=0, color='black', linestyle='-', linewidth=0.5)
    ax2.set_xticklabels(balance_stats['Feature'], rotation=45)
    
    plt.tight_layout()
    plt.show()

Sensitivity Analysis Implementation

We'll implement Rosenbaum bounds sensitivity analysis to assess how robust our findings are to potential unmeasured confounding variables.

def sensitivity_analysis(df, matches, outcome_col, gamma_range=np.arange(1, 2.1, 0.1)):
    results = []
    
    for gamma in gamma_range:
        # Calculate Wilcoxon signed-rank test statistic
        treated_outcomes = df.iloc[[m[0] for m in matches]][outcome_col]
        control_outcomes = df.iloc[[m[1] for m in matches]][outcome_col]
        
        differences = treated_outcomes.values - control_outcomes.values
        abs_diff = np.abs(differences)
        ranks = stats.rankdata(abs_diff)
        
        # Calculate bounds
        p_upper = 1 / (1 + 1/gamma)
        p_lower = 1 / (1 + gamma)
        
        # Store results
        results.append({
            'Gamma': gamma,
            'P_Lower': p_lower,
            'P_Upper': p_upper
        })
    
    return pd.DataFrame(results)

# Example usage
sensitivity_results = sensitivity_analysis(df, results['matches'], 'satisfaction')
print("Sensitivity Analysis Results:")
print(sensitivity_results)

Bootstrapping PSM for Uncertainty Estimation

Bootstrapping allows us to estimate the uncertainty in our treatment effect estimates by resampling our matched pairs. This implementation provides confidence intervals and distribution of treatment effects.

def bootstrap_psm(df, features, treatment_col, outcome_col, n_iterations=1000):
    treatment_effects = []
    
    for _ in range(n_iterations):
        # Sample with replacement
        bootstrap_indices = np.random.choice(
            len(df), size=len(df), replace=True
        )
        bootstrap_df = df.iloc[bootstrap_indices].reset_index(drop=True)
        
        # Run PSM on bootstrap sample
        results = run_psm_analysis(
            bootstrap_df, 
            features, 
            treatment_col, 
            outcome_col
        )
        
        treatment_effects.append(results['treatment_effect']['ATE'])
    
    # Calculate confidence intervals
    ci_lower = np.percentile(treatment_effects, 2.5)
    ci_upper = np.percentile(treatment_effects, 97.5)
    
    return {
        'mean_effect': np.mean(treatment_effects),
        'ci_lower': ci_lower,
        'ci_upper': ci_upper,
        'effects_distribution': treatment_effects
    }

Advanced Matching Strategies

This implementation extends basic PSM with optimal matching and genetic matching algorithms to improve covariate balance and treatment effect estimation accuracy.

def optimal_matching(prop_scores, treatment, max_ratio=1):
    from scipy.optimize import linear_sum_assignment
    
    treated_indices = np.where(treatment == 1)[0]
    control_indices = np.where(treatment == 0)[0]
    
    # Create cost matrix
    cost_matrix = np.zeros((len(treated_indices), len(control_indices)))
    for i, t_idx in enumerate(treated_indices):
        for j, c_idx in enumerate(control_indices):
            cost_matrix[i, j] = abs(prop_scores[t_idx] - prop_scores[c_idx])
    
    # Solve assignment problem
    row_ind, col_ind = linear_sum_assignment(cost_matrix)
    
    # Create matches
    matches = list(zip(treated_indices[row_ind], control_indices[col_ind]))
    
    return matches

Results and Statistical Inference

Implementation of comprehensive statistical testing and inference procedures for PSM results, including multiple hypothesis testing corrections and effect size calculations.

def statistical_inference(df, matches, outcome_col, covariates):
    from scipy import stats
    import statsmodels.api as sm
    
    # Extract matched samples
    treated_indices = [m[0] for m in matches]
    control_indices = [m[1] for m in matches]
    
    # Prepare data for regression
    X = df.iloc[treated_indices + control_indices][covariates]
    y = df.iloc[treated_indices + control_indices][outcome_col]
    treatment = np.concatenate([
        np.ones(len(treated_indices)),
        np.zeros(len(control_indices))
    ])
    
    # Add treatment indicator to covariates
    X = sm.add_constant(X)
    X = pd.concat([X, pd.Series(treatment, name='treatment')], axis=1)
    
    # Fit regression model
    model = sm.OLS(y, X).fit()
    
    # Calculate effect size (Cohen's d)
    treated_outcomes = df.iloc[treated_indices][outcome_col]
    control_outcomes = df.iloc[control_indices][outcome_col]
    
    pooled_std = np.sqrt(
        (np.var(treated_outcomes) + np.var(control_outcomes)) / 2
    )
    cohens_d = (np.mean(treated_outcomes) - np.mean(control_outcomes)) / pooled_std
    
    return {
        'regression_results': model.summary(),
        'effect_size': cohens_d,
        'treatment_coef': model.params['treatment'],
        'p_value': model.pvalues['treatment']
    }

Additional Resources

More in Statistics and analysis

  1. Feature scaling is NOT always necessary9 min

    Feature scaling is a common preprocessing step in machine learning, but it's not always required. In this presentation, we'll explore scenarios where feature...

  2. Detecting Anomalies and Outliers in Data10 min

    Outliers are data points that significantly differ from other observations. Flagging outliers is crucial for data analysis and model performance.

  3. Feature Scaling Techniques for Outlier-Robust Data6 min

    Feature scaling is crucial when dealing with datasets containing outliers. This presentation explores various techniques to scale features effectively...

  4. Cosine Similarity in Vector Stores7 min

    Cosine similarity is a metric used to measure the similarity between two non-zero vectors in a multi-dimensional space. It calculates the cosine of the angle...

  5. Independent and Dependent Variables in Predictive Modeling8 min

    Independent variables are features used as input to predict an outcome. They are the factors that we manipulate or control in an experiment or analysis to...

  6. Correlation vs Causation Understanding the Difference14 min

    Correlation measures the statistical relationship between two variables, indicating how they move together. It quantifies the strength and direction of their...