The Replication Crisis

"There are lies, damned lies, and statistics." — Often attributed to Benjamin Disraeli, but more accurately a commentary on the misuse of statistics

The Evidence

Reproducibility Project: Psychology

Other Replication Efforts

Causes of Non-Replication

P-Hacking

Common forms of p-hacking:

1. Testing multiple outcomes but reporting only the significant ones
2. Adding or removing covariates until $p < 0.05$
3. Running multiple statistical tests (t-test, ANOVA, regression) and choosing the one that works
4. Collecting more data until significance is reached (optional stopping)
5. Excluding outliers based on different criteria until significance
6. Analyzing subgroups until a significant effect appears

HARKing

HARKing is particularly insidious because:

• It is undetectable from the published paper alone
• It inflates the apparent confirmatory nature of research
• It makes the literature appear more theoretically driven than it actually is

The Garden of Forking Paths

Even a well-intentioned researcher can arrive at $p < 0.05$ through legitimate analytical choices, because the number of possible analyses is enormous and each choice has a small but real chance of producing a false positive.

Statistical Remedies

Pre-Registration

Pre-registration addresses:

• HARKing (hypotheses cannot be fabricated after seeing results)
• P-hacking (analysis plan is fixed before data collection)
• Publication bias (all registered studies are discoverable regardless of results)

Multi-Lab Replications

Advantages:

• Large aggregate sample sizes provide high statistical power
• Cross-cultural variation tests generalizability
• Independent teams eliminate single-lab bias
• Pre-registered protocols prevent analytical flexibility

Effect Sizes and Confidence Intervals

Bayesian Analysis

The Bayes factor compares the evidence for two hypotheses:

Open Science Reforms

Registered Reports

Open Data and Open Materials

Incentive Reform

Quantifying the Impact

Inflation of Effect Sizes

False Discovery Rate

Python Implementation

import numpy as np
from scipy import stats
import matplotlib.pyplot as plt

# --- Simulate p-hacking via optional stopping ---
def simulate_optional_stopping(true_effect=0, n_max=500, alpha=0.05, n_sims=10000):
    """
    Show how optional stopping inflates false positive rates.
    Test after each new observation; stop when p < alpha or n_max reached.
    """
    false_positives = 0
    sample_sizes = []
    
    for _ in range(n_sims):
        for n in range(10, n_max + 1, 5):
            data1 = np.random.normal(0, 1, n)
            data2 = np.random.normal(true_effect, 1, n)
            t_stat, p_val = stats.ttest_ind(data1, data2)
            
            if p_val < alpha:
                false_positives += 1
                sample_sizes.append(n)
                break
    
    fpr = false_positives / n_sims
    avg_n = np.mean(sample_sizes) if sample_sizes else n_max
    return fpr, avg_n

# Test at different true effects
print("=== Optional Stopping Simulation ===")
print(f"{'True Effect':<15} {'FPR (nominal α=0.05)':<25} {'Avg Sample Size':<15}")
print("-" * 55)

for true_eff in [0, 0.1, 0.2, 0.3]:
    fpr, avg_n = simulate_optional_stopping(true_effect=true_eff)
    print(f"{true_eff:<15.1f} {fpr:<25.3f} {avg_n:<15.1f}")

# --- Simulate p-hacking via multiple testing ---
def simulate_p_hacking(n_tests=20, n_sims=10000, alpha=0.05):
    """
    Show how testing multiple outcomes inflates false positive rate.
    """
    # No p-hacking: test only the primary outcome
    fp_no_hack = sum(
        1 for _ in range(n_sims) 
        if stats.ttest_ind(
            np.random.normal(0, 1, 100),
            np.random.normal(0, 1, 100)
        )[1] < alpha
    ) / n_sims
    
    # P-hacking: test n_tests outcomes, report smallest p-value
    fp_hack = 0
    for _ in range(n_sims):
        p_vals = [stats.ttest_ind(
            np.random.normal(0, 1, 100),
            np.random.normal(0, 1, 100)
        )[1] for _ in range(n_tests)]
        # Bonferroni-corrected threshold
        if min(p_vals) < alpha / n_tests:
            fp_hack += 1
    fp_hack_corrected = fp_hack / n_sims
    
    # P-hacking without correction
    fp_hack_uncorrected = sum(
        1 for _ in range(n_sims)
        if min([stats.ttest_ind(
            np.random.normal(0, 1, 100),
            np.random.normal(0, 1, 100)
        )[1] for _ in range(n_tests)]) < alpha
    ) / n_sims
    
    return fp_no_hack, fp_hack_uncorrected, fp_hack_corrected

print("\n=== Multiple Testing Simulation (20 outcomes) ===")
fp1, fp2, fp3 = simulate_p_hacking()
print(f"Single test (no hacking):           FPR = {fp1:.3f}")
print(f"20 tests, no correction (hacking):  FPR = {fp2:.3f}")
print(f"20 tests, Bonferroni corrected:      FPR = {fp3:.3f}")

# --- Funnel plot of published vs all studies ---
np.random.seed(42)
true_effect = 0.3
true_se = np.random.uniform(0.05, 0.3, 200)
true_effects = np.random.normal(true_effect, 0.1, 200)
observed_effects = np.random.normal(true_effects, true_se)
p_values = 2 * (1 - stats.norm.cdf(np.abs(observed_effects / true_se)))

# "Published" studies (significant only)
published = p_values < 0.05

fig, axes = plt.subplots(1, 2, figsize=(14, 6))

# All studies
axes[0].scatter(observed_effects, true_se, s=20, alpha=0.5, c='gray')
axes[0].axvline(x=true_effect, color='red', linestyle='--', label='True effect')
axes[0].set_xlabel('Observed Effect Size')
axes[0].set_ylabel('Standard Error')
axes[0].set_title(f'All Studies (n={len(observed_effects)})')
axes[0].invert_yaxis()
axes[0].legend()
axes[0].grid(True, alpha=0.3)

# Published only (publication bias)
axes[1].scatter(observed_effects[published], true_se[published], s=30, alpha=0.7, c='blue')
axes[1].scatter(observed_effects[~published], true_se[~published], s=15, alpha=0.3, c='gray', 
                label='Not published')
axes[1].axvline(x=true_effect, color='red', linestyle='--', label='True effect')
axes[1].set_xlabel('Observed Effect Size')
axes[1].set_ylabel('Standard Error')
axes[1].set_title(f'Published Only (n={published.sum()})')
axes[1].invert_yaxis()
axes[1].legend()
axes[1].grid(True, alpha=0.3)

plt.suptitle('Publication Bias: Funnel Plot Distortion', fontsize=14)
plt.tight_layout()
plt.savefig('publication_bias.png', dpi=150)
plt.show()

# Effect size inflation
print(f"\n=== Effect Size Inflation ===")
print(f"True effect:                          {true_effect:.3f}")
print(f"Mean of all observed effects:         {np.mean(observed_effects):.3f}")
print(f"Mean of published (significant) effects: {np.mean(observed_effects[published]):.3f}")
print(f"Inflation factor:                     {np.mean(observed_effects[published])/true_effect:.2f}x")

# --- Bayes Factor computation (simplified using JZS prior) ---
def bayes_factor_ttest(x, y, r=np.sqrt(2)/2):
    """
    Compute Bayes Factor for independent samples t-test using JZS prior.
    Rouder et al. (2009) method.
    """
    n1, n2 = len(x), len(y)
    t_stat, _ = stats.ttest_ind(x, y)
    df = n1 + n2 - 2
    
    # JZS Bayes Factor
    BF01 = (1 + t_stat**2 / df)**(-(n1 + n2) / 2) * \
           (1 + t_stat**2 / (df * (1 + r**2)))**(df / 2 + 0.5) * \
           np.sqrt(1 / (1 + r**2))
    
    return 1 / BF01  # BF10 (evidence for alternative)

print("\n=== Bayes Factor vs p-Value ===")
print(f"{'Scenario':<30} {'p-value':<10} {'BF10':<10} {'Interpretation':<20}")
print("-" * 70)

scenarios = [
    ("Large effect, n=20", np.random.normal(0.8, 1, 20), np.random.normal(0, 1, 20)),
    ("Small effect, n=20", np.random.normal(0.2, 1, 20), np.random.normal(0, 1, 20)),
    ("Large effect, n=200", np.random.normal(0.3, 1, 200), np.random.normal(0, 1, 200)),
    ("No effect, n=200", np.random.normal(0, 1, 200), np.random.normal(0, 1, 200)),
]

for name, x, y in scenarios:
    t_stat, p_val = stats.ttest_ind(x, y)
    bf = bayes_factor_ttest(x, y)
    
    if bf > 100:
        interp = "Extreme for H1"
    elif bf > 10:
        interp = "Strong for H1"
    elif bf > 3:
        interp = "Moderate for H1"
    elif bf > 1/3:
        interp = "Anecdotal"
    elif bf > 1/10:
        interp = "Moderate for H0"
    else:
        interp = "Strong for H0"
    
    print(f"{name:<30} {p_val:<10.4f} {bf:<10.2f} {interp:<20}")

Key Takeaways

Next Steps

• -> Systematic Review Methodology
• -> Meta-Analysis
• -> Bayesian Statistics