Simulate p-values with varying sample sizes or effect sizes

p_val_sim.py

#!/usr/bin/env python3

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


def simulate_p_values(effect_size, sample_sizes, iterations=200):
    """Simulate p-values with varying sample sizes"""
    mean_p_values = []

    # Create range of sample sizes to simulate through
    for n in sample_sizes:
        p_vals = []
        for _ in range(iterations):
            # Generate two samples with a known difference (effect_size)
            group_a = np.random.normal(0, 1, n)
            group_b = np.random.normal(effect_size, 1, n)

            # perform independent t-test
            _, p_val = stats.ttest_ind(group_a, group_b)
            p_vals.append(p_val)

        # Average the p-values across iterations for a smoother trend
        mean_p_values.append(np.mean(p_vals))

    return mean_p_values


def simulate_p_values_e(effect_size, sample_size, iterations=200):
    """Simulate p-values with varying effect sizes"""
    mean_p_values = []

    # Create range of sample sizes to simulate through
    for n in effect_sizes:
        p_vals = []
        for _ in range(iterations):
            # Generate two samples with a known difference (effect_size)
            group_a = np.random.normal(0, 1, sample_size)
            group_b = np.random.normal(n, 1, sample_size)

            # perform independent t-test
            _, p_val = stats.ttest_ind(group_a, group_b)
            p_vals.append(p_val)

        # Average the p-values across iterations for a smoother trend
        mean_p_values.append(np.mean(p_vals))

    return mean_p_values


# Parameters
# effect_size = 0.2  # A "small" signal
effect_size = 0.01
# sample_sizes = np.arange(10, 2001, 50)  # Sample sizes 10 to 2000, inc by 50
sample_sizes = np.arange(10, 1000001, 100000)
sample_size = 10000
effect_sizes = np.arange(10, 10000, 100) / 100000

# Run simulation
avg_p_values = simulate_p_values(effect_size, sample_sizes)
# avg_p_values = simulate_p_values_e(effect_sizes, sample_size)

# Plotting results
plt.figure(figsize=(10, 6))
plt.plot(
    sample_sizes,
    avg_p_values,
    marker='o',
    color='teal',
    label='Average p-value'
)
plt.axhline(
    y=0.05,
    color='red',
    linestyle='--',
    label='Significance threshold (alpha=0.05)'
)
plt.title(f'Relationship between sample size and p-value (effect size = {effect_size})')
plt.xlabel('Sample size (n per group)')
plt.ylabel('Mean p-value (log scale)')
plt.yscale('log')  # Help visualize over magnitudes
plt.grid(True, which='both', ls='-', alpha=0.5)
plt.legend()
plt.show()