Hypothesis Testing

import pandas as pd
import numpy as np
import seaborn as sns
import matplotlib.pyplot as plt
from scipy import stats
from scipy.stats import (
    ttest_ind, ttest_rel, ttest_1samp,
    mannwhitneyu, wilcoxon, kruskal,
    chi2_contingency, fisher_exact,
    f_oneway, shapiro, levene, bartlett,
    spearmanr, pearsonr, kendalltau,
    kstest, normaltest, anderson
)
import warnings
warnings.filterwarnings('ignore')
 
# 1. Titanic Dataset
titanic = sns.load_dataset('titanic')
print(f"Titanic: {titanic.shape}")
 
# 2. Iris Dataset
iris = sns.load_dataset('iris')
print(f"Iris: {iris.shape}")
 
# 3. Tips Dataset
tips = sns.load_dataset('tips')
print(f"Tips: {tips.shape}")
 
# 4. Diamonds Dataset (sampled)
diamonds = sns.load_dataset('diamonds').sample(n=1000, random_state=42)
print(f"Diamonds: {diamonds.shape}")

Data Type?
├── Continuous
│   ├── Normal distribution → Parametric tests (t-test, ANOVA)
│   └── Non-normal distribution → Non-parametric tests (Mann-Whitney, Kruskal-Wallis)
│
└── Categorical
    ├── 2×2 table (small sample) → Fisher's Exact Test
    └── Otherwise → Chi-Square Test

# Titanic: Normality test for age distribution
ages = titanic['age'].dropna()
 
stat, p_value = shapiro(ages)
 
print("=== Shapiro-Wilk Normality Test ===")
print(f"Data: Titanic passenger ages (n={len(ages)})")
print(f"Statistic: {stat:.4f}")
print(f"p-value: {p_value:.4f}")
print(f"Conclusion: {'Follows normal distribution' if p_value >= 0.05 else 'Does not follow normal distribution'}")
 
# Visualization
fig, axes = plt.subplots(1, 2, figsize=(10, 4))
 
axes[0].hist(ages, bins=30, edgecolor='black', alpha=0.7)
axes[0].set_title('Age Distribution')
axes[0].set_xlabel('Age')
 
stats.probplot(ages, dist="norm", plot=axes[1])
axes[1].set_title('Q-Q Plot')
 
plt.tight_layout()
plt.show()

# Iris: Normality test for petal length
petal_length = iris['petal_length']
 
stat, p_value = normaltest(petal_length)
 
print("=== D'Agostino-Pearson Normality Test ===")
print(f"Data: Iris petal length (n={len(petal_length)})")
print(f"Statistic: {stat:.4f}")
print(f"p-value: {p_value:.4f}")
print(f"Conclusion: {'Follows normal distribution' if p_value >= 0.05 else 'Does not follow normal distribution'}")

# Tips: Normality test for tip amounts
tip_values = tips['tip']
 
# Test after standardization
tip_standardized = (tip_values - tip_values.mean()) / tip_values.std()
stat, p_value = kstest(tip_standardized, 'norm')
 
print("=== Kolmogorov-Smirnov Normality Test ===")
print(f"Data: Tips tip amounts (n={len(tip_values)})")
print(f"Statistic: {stat:.4f}")
print(f"p-value: {p_value:.4f}")
print(f"Conclusion: {'Follows normal distribution' if p_value >= 0.05 else 'Does not follow normal distribution'}")

# Diamonds: Normality test for prices
prices = diamonds['price']
 
result = anderson(prices, dist='norm')
 
print("=== Anderson-Darling Normality Test ===")
print(f"Data: Diamonds prices (n={len(prices)})")
print(f"Statistic: {result.statistic:.4f}")
print("\nCritical values by significance level:")
for i, (cv, sl) in enumerate(zip(result.critical_values, result.significance_level)):
    result_str = "Reject" if result.statistic > cv else "Accept"
    print(f"  {sl}%: Critical value = {cv:.4f} → H₀ {result_str}")

# Titanic: Compare age variance by survival status
survived_ages = titanic[titanic['survived'] == 1]['age'].dropna()
died_ages = titanic[titanic['survived'] == 0]['age'].dropna()
 
stat, p_value = levene(survived_ages, died_ages)
 
print("=== Levene's Test for Equal Variances ===")
print(f"Survived age variance: {survived_ages.var():.2f}")
print(f"Died age variance: {died_ages.var():.2f}")
print(f"Statistic: {stat:.4f}")
print(f"p-value: {p_value:.4f}")
print(f"Conclusion: {'Equal variance assumption satisfied' if p_value >= 0.05 else 'Equal variance assumption not satisfied'}")

# Iris: Compare sepal length variance by species
setosa = iris[iris['species'] == 'setosa']['sepal_length']
versicolor = iris[iris['species'] == 'versicolor']['sepal_length']
virginica = iris[iris['species'] == 'virginica']['sepal_length']
 
stat, p_value = bartlett(setosa, versicolor, virginica)
 
print("=== Bartlett's Test for Equal Variances ===")
print(f"Setosa variance: {setosa.var():.4f}")
print(f"Versicolor variance: {versicolor.var():.4f}")
print(f"Virginica variance: {virginica.var():.4f}")
print(f"Statistic: {stat:.4f}")
print(f"p-value: {p_value:.4f}")
print(f"Conclusion: {'Equal variance assumption satisfied' if p_value >= 0.05 else 'Equal variance assumption not satisfied'}")

# Tips: Test if mean tip is $3
# Situation: Restaurant manager claims "Our average tip is $3". Verify if this is true.
tip_values = tips['tip']
hypothesized_mean = 3.0
 
stat, p_value = ttest_1samp(tip_values, hypothesized_mean)
 
print("=== One-Sample t-test ===")
print(f"H₀: Mean tip = ${hypothesized_mean:.2f}")
print(f"H₁: Mean tip ≠ ${hypothesized_mean:.2f}")
print(f"\nSample mean: ${tip_values.mean():.2f}")
print(f"Sample standard deviation: ${tip_values.std():.2f}")
print(f"t-statistic: {stat:.4f}")
print(f"p-value: {p_value:.4f}")
print(f"\nConclusion: {'Mean tip differs from $3' if p_value < 0.05 else 'Mean tip can be considered equal to $3'}")

# Titanic: Compare fare by gender
# Situation: Analyze if there were differences in fares paid by males and females on the Titanic
male_fare = titanic[titanic['sex'] == 'male']['fare'].dropna()
female_fare = titanic[titanic['sex'] == 'female']['fare'].dropna()
 
# Test equal variance assumption
_, levene_p = levene(male_fare, female_fare)
equal_var = levene_p >= 0.05
 
# t-test (Welch's t-test if unequal variance)
stat, p_value = ttest_ind(male_fare, female_fare, equal_var=equal_var)
 
print("=== Independent Samples t-test ===")
print(f"H₀: Male mean fare = Female mean fare")
print(f"H₁: Male mean fare ≠ Female mean fare")
print(f"\nMale mean fare: ${male_fare.mean():.2f} (n={len(male_fare)})")
print(f"Female mean fare: ${female_fare.mean():.2f} (n={len(female_fare)})")
print(f"Difference: ${female_fare.mean() - male_fare.mean():.2f}")
print(f"\nEqual variance assumption: {'Satisfied (Student t-test)' if equal_var else 'Not satisfied (Welch t-test used)'}")
print(f"t-statistic: {stat:.4f}")
print(f"p-value: {p_value:.4f}")
print(f"\nConclusion: {'Significant difference in fare by gender' if p_value < 0.05 else 'No difference in fare by gender'}")

# Simulation: A/B test conversion rate (same users experience both UIs)
# Situation: Show 100 users old UI and new UI sequentially and measure click rates
np.random.seed(42)
n_users = 100
 
# Before: Old UI click rate
before = np.random.beta(2, 8, n_users)  # Mean about 20%
# After: New UI click rate (slight improvement)
after = before + np.random.normal(0.05, 0.03, n_users)
after = np.clip(after, 0, 1)
 
stat, p_value = ttest_rel(before, after)
 
print("=== Paired Samples t-test ===")
print(f"H₀: No change in conversion rate (new UI has no effect)")
print(f"H₁: Change in conversion rate (new UI has effect)")
print(f"\nBefore (old UI) mean: {before.mean():.4f} ({before.mean()*100:.1f}%)")
print(f"After (new UI) mean: {after.mean():.4f} ({after.mean()*100:.1f}%)")
print(f"Mean difference: {(after - before).mean():.4f} (+{(after - before).mean()*100:.1f}%p)")
print(f"\nt-statistic: {stat:.4f}")
print(f"p-value: {p_value:.4f}")
print(f"\nConclusion: {'New UI significantly improved' if p_value < 0.05 else 'No significant change'}")

# Diamonds: Compare prices by cut quality (Ideal vs Good)
# Situation: Compare price distributions of diamonds with Ideal vs Good cut quality
# Price data generally doesn't follow normal distribution, so use non-parametric test
ideal_price = diamonds[diamonds['cut'] == 'Ideal']['price']
good_price = diamonds[diamonds['cut'] == 'Good']['price']
 
stat, p_value = mannwhitneyu(ideal_price, good_price, alternative='two-sided')
 
print("=== Mann-Whitney U Test ===")
print(f"H₀: Ideal and Good cut have the same price distribution")
print(f"H₁: Price distributions are different")
print(f"\nIdeal median: ${ideal_price.median():,.2f} (n={len(ideal_price)})")
print(f"Good median: ${good_price.median():,.2f} (n={len(good_price)})")
print(f"\nU-statistic: {stat:,.2f}")
print(f"p-value: {p_value:.4f}")
print(f"\nConclusion: {'Price distributions significantly different' if p_value < 0.05 else 'No difference in price distributions'}")

# Tips: Compare lunch vs dinner tip rates (same waiters)
# Situation: Test if 50 waiters receive different tip rates at lunch vs dinner
np.random.seed(42)
n_waiters = 50
 
lunch_tip_rate = np.random.uniform(0.12, 0.22, n_waiters)
dinner_tip_rate = lunch_tip_rate + np.random.normal(0.02, 0.03, n_waiters)
 
stat, p_value = wilcoxon(lunch_tip_rate, dinner_tip_rate)
 
print("=== Wilcoxon Signed-Rank Test ===")
print(f"H₀: No difference between lunch and dinner tip rates")
print(f"H₁: Difference between lunch and dinner tip rates")
print(f"\nLunch tip rate median: {np.median(lunch_tip_rate)*100:.1f}%")
print(f"Dinner tip rate median: {np.median(dinner_tip_rate)*100:.1f}%")
print(f"\nW-statistic: {stat:.2f}")
print(f"p-value: {p_value:.4f}")
print(f"\nConclusion: {'Dinner tip rate significantly higher' if p_value < 0.05 else 'No difference'}")

# Iris: Compare petal width by species
# Situation: Test if petal width distributions differ among three iris species (setosa, versicolor, virginica)
setosa_pw = iris[iris['species'] == 'setosa']['petal_width']
versicolor_pw = iris[iris['species'] == 'versicolor']['petal_width']
virginica_pw = iris[iris['species'] == 'virginica']['petal_width']
 
stat, p_value = kruskal(setosa_pw, versicolor_pw, virginica_pw)
 
print("=== Kruskal-Wallis H Test ===")
print(f"H₀: All species have the same petal width distribution")
print(f"H₁: At least one species is different")
print(f"\nSetosa median: {setosa_pw.median():.2f}cm")
print(f"Versicolor median: {versicolor_pw.median():.2f}cm")
print(f"Virginica median: {virginica_pw.median():.2f}cm")
print(f"\nH-statistic: {stat:.4f}")
print(f"p-value: {p_value:.6f}")
print(f"\nConclusion: {'Significant difference among species' if p_value < 0.05 else 'No difference'}")

# Tips: Compare total bill amounts by day of week
# Situation: Analyze if customers' average payment amounts differ by day of the week
thur = tips[tips['day'] == 'Thur']['total_bill']
fri = tips[tips['day'] == 'Fri']['total_bill']
sat = tips[tips['day'] == 'Sat']['total_bill']
sun = tips[tips['day'] == 'Sun']['total_bill']
 
stat, p_value = f_oneway(thur, fri, sat, sun)
 
print("=== One-Way ANOVA ===")
print(f"H₀: Average payment amounts are the same for all days")
print(f"H₁: At least one day is different")
print(f"\nAverage payment by day:")
print(f"  Thursday: ${thur.mean():.2f} (n={len(thur)})")
print(f"  Friday: ${fri.mean():.2f} (n={len(fri)})")
print(f"  Saturday: ${sat.mean():.2f} (n={len(sat)})")
print(f"  Sunday: ${sun.mean():.2f} (n={len(sun)})")
print(f"\nF-statistic: {stat:.4f}")
print(f"p-value: {p_value:.4f}")
print(f"\nConclusion: {'Significant difference by day → Post-hoc test needed' if p_value < 0.05 else 'No difference by day'}")

import statsmodels.api as sm
from statsmodels.formula.api import ols
 
# Tips: Effect of gender and smoking status on tips
# Situation: Analyze if tip amounts differ by gender and smoking status, and if there's a combination effect
model = ols('tip ~ C(sex) + C(smoker) + C(sex):C(smoker)', data=tips).fit()
anova_table = sm.stats.anova_lm(model, typ=2)
 
print("=== Two-Way ANOVA ===")
print(f"Dependent Variable: Tip amount")
print(f"Factor 1: Gender (sex)")
print(f"Factor 2: Smoking status (smoker)")
print(f"\n{anova_table.round(4)}")
print("\nInterpretation:")
for idx, row in anova_table.iterrows():
    if idx != 'Residual':
        sig = "Significant ***" if row['PR(>F)'] < 0.001 else \
              "Significant **" if row['PR(>F)'] < 0.01 else \
              "Significant *" if row['PR(>F)'] < 0.05 else "Not significant"
        print(f"  {idx}: p = {row['PR(>F)']:.4f} → {sig}")

# Titanic: Relationship between gender and survival
# Situation: Analyze if survival rates differed by gender during the Titanic sinking ("Women and children first" rule)
contingency = pd.crosstab(titanic['sex'], titanic['survived'])
print("Contingency Table:")
print(contingency)
print()
 
chi2, p_value, dof, expected = chi2_contingency(contingency)
 
print("=== Chi-Square Test of Independence ===")
print(f"H₀: Gender and survival are independent (no relationship)")
print(f"H₁: Gender and survival are associated")
print(f"\nChi-square statistic: {chi2:.4f}")
print(f"Degrees of freedom: {dof}")
print(f"p-value: {p_value:.6f}")
print(f"\nExpected frequencies (if independent, these would be expected):")
print(pd.DataFrame(expected,
                   index=contingency.index,
                   columns=contingency.columns).round(1))
print(f"\nConclusion: {'Gender and survival are strongly associated (higher female survival rate)' if p_value < 0.05 else 'Independent'}")

# Titanic: Test if cabin class distribution is uniform
# Situation: Check if Titanic passengers were evenly distributed across 1st, 2nd, 3rd class
observed = titanic['pclass'].value_counts().sort_index()
n = len(titanic)
expected = np.array([n/3, n/3, n/3])  # Expected uniform distribution
 
chi2, p_value = stats.chisquare(observed, expected)
 
print("=== Chi-Square Goodness of Fit Test ===")
print(f"H₀: Cabin classes are evenly distributed (each 33.3%)")
print(f"H₁: Not evenly distributed")
print(f"\nObserved frequencies:")
for cls, count in observed.items():
    print(f"  Class {cls}: {count} passengers ({count/n*100:.1f}%)")
print(f"\nExpected frequency (uniform distribution): {n/3:.0f} each (33.3%)")
print(f"\nChi-square statistic: {chi2:.4f}")
print(f"p-value: {p_value:.6f}")
print(f"\nConclusion: {'Not uniform - 3rd class is the majority' if p_value < 0.05 else 'Uniform distribution'}")

# Small-scale clinical trial simulation
# Situation: Pilot study with 20 subjects. Compare cure rates between drug treatment group (10) and placebo group (10)
contingency = np.array([[8, 2],   # Treatment group: 8 cured, 2 not cured
                        [3, 7]])  # Control group: 3 cured, 7 not cured
 
odds_ratio, p_value = fisher_exact(contingency)
 
print("=== Fisher's Exact Test ===")
print("Situation: Small-scale clinical trial (n=20)")
print("\nContingency Table:")
print("         Cured  Not Cured")
print(f"Treatment     {contingency[0,0]}      {contingency[0,1]}")
print(f"Control       {contingency[1,0]}      {contingency[1,1]}")
print(f"\nTreatment cure rate: {contingency[0,0]/contingency[0].sum()*100:.0f}%")
print(f"Control cure rate: {contingency[1,0]/contingency[1].sum()*100:.0f}%")
print(f"\nOdds Ratio: {odds_ratio:.4f}")
print(f"p-value: {p_value:.4f}")
print(f"\nConclusion: {'Treatment effect exists' if p_value < 0.05 else 'No treatment effect'}")
print(f"\nInterpretation: Treatment group has {odds_ratio:.1f}x higher cure odds than control group")

from statsmodels.stats.contingency_tables import mcnemar
 
# Marketing campaign before/after purchase behavior change
# Situation: Track purchase status of 100 customers before and after campaign
# [Bought before/Bought after, Bought before/Didn't buy after]
# [Didn't buy before/Bought after, Didn't buy before/Didn't buy after]
table = np.array([[45, 15],   # Bought before & after / Bought before, didn't buy after
                  [35, 5]])   # Didn't buy before, bought after / Didn't buy before & after
 
result = mcnemar(table, exact=True)
 
print("=== McNemar's Test ===")
print("Situation: Purchase behavior change of 100 customers before/after campaign")
print("\nPaired Table:")
print("              Bought After  Didn't Buy After")
print(f"Bought Before      {table[0,0]}            {table[0,1]}")
print(f"Didn't Buy Before      {table[1,0]}             {table[1,1]}")
print(f"\nChange Analysis:")
print(f"  New purchase (Didn't buy → Bought): {table[1,0]} people")
print(f"  Churned (Bought → Didn't buy): {table[0,1]} people")
print(f"  No change: {table[0,0] + table[1,1]} people")
print(f"\np-value: {result.pvalue:.4f}")
print(f"\nConclusion: {'Campaign significantly changed purchase behavior (New purchases > Churn)' if result.pvalue < 0.05 else 'No significant change'}")

# Diamonds: Correlation between carat and price
# Situation: Analyze if there is a linear relationship between diamond carat (weight) and price
carat = diamonds['carat']
price = diamonds['price']
 
r, p_value = pearsonr(carat, price)
 
print("=== Pearson Correlation Analysis ===")
print(f"H₀: There is no correlation between carat and price (r = 0)")
print(f"H₁: There is correlation between carat and price (r ≠ 0)")
print(f"\nPearson r: {r:.4f}")
print(f"p-value: {p_value:.6f}")
print(f"Coefficient of determination (R²): {r**2:.4f} → {r**2*100:.1f}% of price variation explained by carat")
print(f"\nCorrelation strength interpretation:")
print(f"  |r| < 0.3: Weak correlation")
print(f"  0.3 <= |r| < 0.7: Moderate correlation")
print(f"  |r| >= 0.7: Strong correlation")
print(f"\nCurrent |r| = {abs(r):.4f} → Strong positive correlation (carat↑ → price↑)")

# Tips: Rank correlation between total bill and tip
# Situation: Check if higher bills tend to have higher tips (even if not exactly proportional)
total_bill = tips['total_bill']
tip = tips['tip']
 
rho, p_value = spearmanr(total_bill, tip)
r_pearson, _ = pearsonr(total_bill, tip)
 
print("=== Spearman Rank Correlation Analysis ===")
print(f"H₀: No monotonic relationship between bill and tip")
print(f"H₁: Monotonic relationship exists between bill and tip")
print(f"\nSpearman rho: {rho:.4f}")
print(f"(Comparison) Pearson r: {r_pearson:.4f}")
print(f"p-value: {p_value:.6f}")
print(f"\nConclusion: {'Significant monotonic relationship - higher bills tend to have higher tips' if p_value < 0.05 else 'No relationship'}")

# Titanic: Relationship between cabin class and age
# Situation: Check if 1st class passengers tend to be older (ordinal vs continuous)
pclass = titanic['pclass'].dropna()
age = titanic['age'].dropna()
 
# Align indices
common_idx = pclass.index.intersection(age.index)
pclass_aligned = pclass.loc[common_idx]
age_aligned = age.loc[common_idx]
 
tau, p_value = kendalltau(pclass_aligned, age_aligned)
 
print("=== Kendall's Tau Correlation Analysis ===")
print(f"H₀: No relationship between cabin class and age")
print(f"H₁: Relationship exists between cabin class and age")
print(f"\nKendall tau: {tau:.4f}")
print(f"p-value: {p_value:.4f}")
print(f"\nConclusion: {'Significant relationship' if p_value < 0.05 else 'No relationship'}")
if tau < 0:
    print(f"Interpretation: tau < 0 means cabin class↓(1st class) → age↑ (older passengers in premium cabins)")

def cohens_d(group1, group2):
    n1, n2 = len(group1), len(group2)
    var1, var2 = group1.var(), group2.var()
    pooled_std = np.sqrt(((n1-1)*var1 + (n2-1)*var2) / (n1+n2-2))
    return (group1.mean() - group2.mean()) / pooled_std
 
# Titanic: Effect size of age difference between survivors vs non-survivors
# Situation: Relationship between age and survival. Even if p-value is significant, is it a meaningful difference?
survived_ages = titanic[titanic['survived'] == 1]['age'].dropna()
died_ages = titanic[titanic['survived'] == 0]['age'].dropna()
 
d = cohens_d(survived_ages, died_ages)
t_stat, p_value = ttest_ind(survived_ages, died_ages)
 
print("=== Effect Size Analysis ===")
print(f"Survivor mean age: {survived_ages.mean():.2f} years")
print(f"Non-survivor mean age: {died_ages.mean():.2f} years")
print(f"Difference: {abs(survived_ages.mean() - died_ages.mean()):.2f} years")
print(f"\nt-statistic: {t_stat:.4f}")
print(f"p-value: {p_value:.4f} → {'Significant' if p_value < 0.05 else 'Not significant'}")
print(f"Cohen's d: {d:.4f}")
print(f"\nEffect size interpretation:")
print(f"  |d| < 0.2: Small effect")
print(f"  0.2 <= |d| < 0.5: Medium effect")
print(f"  0.5 <= |d| < 0.8: Large effect")
print(f"  |d| >= 0.8: Very large effect")
print(f"\nCurrent: |d| = {abs(d):.4f} → ", end="")
if abs(d) >= 0.8:
    print("Very large effect")
elif abs(d) >= 0.5:
    print("Large effect")
elif abs(d) >= 0.2:
    print("Medium effect")
else:
    print("Small effect → Statistically significant but practical meaning is limited!")

def cramers_v(contingency_table):
    chi2 = chi2_contingency(contingency_table)[0]
    n = contingency_table.sum().sum()
    min_dim = min(contingency_table.shape) - 1
    return np.sqrt(chi2 / (n * min_dim))
 
# Titanic: Strength of association between gender and survival
contingency = pd.crosstab(titanic['sex'], titanic['survived'])
v = cramers_v(contingency)
chi2, p_value, _, _ = chi2_contingency(contingency)
 
print("=== Cramer's V (Association Strength) ===")
print(f"Chi-square = {chi2:.4f}, p-value = {p_value:.6f}")
print(f"Cramer's V = {v:.4f}")
print(f"\nInterpretation guidelines:")
print(f"  V < 0.1: Negligible")
print(f"  0.1 <= V < 0.3: Weak association")
print(f"  0.3 <= V < 0.5: Moderate association")
print(f"  V >= 0.5: Strong association")
print(f"\nCurrent: V = {v:.4f} → ", end="")
if v >= 0.5:
    print("Strong association → Gender strongly affected survival!")
elif v >= 0.3:
    print("Moderate association")
elif v >= 0.1:
    print("Weak association")
else:
    print("Negligible")

from statsmodels.stats.multitest import multipletests
 
# Multiple A/B test results
# Situation: Testing 5 UI elements simultaneously. Which ones have real effects?
p_values = [0.03, 0.04, 0.01, 0.08, 0.002]
test_names = ['Button Color', 'Headline', 'CTA Position', 'Image', 'Price Display']
 
# Bonferroni correction
rejected, corrected_p, _, _ = multipletests(p_values, method='bonferroni')
 
print("=== Multiple Testing Correction (Bonferroni) ===")
print(f"Number of tests: {len(p_values)}")
print(f"Corrected significance level: 0.05 / {len(p_values)} = {0.05/len(p_values):.3f}")
print(f"\n{'Test':<15} {'Original p-value':<18} {'Corrected p-value':<18} {'Conclusion'}")
print("-" * 70)
for name, p, cp, rej in zip(test_names, p_values, corrected_p, rejected):
    result = "Significant" if rej else "Not significant"
    print(f"{name:<15} {p:<18.4f} {min(cp, 1.0):<18.4f} {result}")

# Benjamini-Hochberg correction
rejected_bh, corrected_p_bh, _, _ = multipletests(p_values, method='fdr_bh')
 
print("=== Multiple Testing Correction (Benjamini-Hochberg FDR) ===")
print(f"\n{'Test':<15} {'Original p-value':<18} {'Corrected p-value':<18} {'Conclusion'}")
print("-" * 70)
for name, p, cp, rej in zip(test_names, p_values, corrected_p_bh, rejected_bh):
    result = "Significant" if rej else "Not significant"
    print(f"{name:<15} {p:<18.4f} {cp:<18.4f} {result}")
 
print(f"\nComparison:")
print(f"  Significant with Bonferroni: {sum(rejected)}")
print(f"  Significant with FDR(BH): {sum(rejected_bh)}")
print(f"\n→ FDR is less conservative, allowing more discoveries")
print(f"   However, about 5% of these may be false positives")

from statsmodels.stats.power import TTestIndPower
 
power_analysis = TTestIndPower()
 
# Scenario: Effect size 0.3, Power 80%, Significance level 5%
# Situation: "Expecting the new UI to improve conversion rate moderately (d=0.3).
#           How many people are needed to detect this effect with 80% probability?"
effect_size = 0.3
alpha = 0.05
power = 0.8
 
n = power_analysis.solve_power(effect_size=effect_size,
                                alpha=alpha,
                                power=power,
                                ratio=1.0,
                                alternative='two-sided')
 
print("=== Sample Size Calculation ===")
print(f"Target effect size (Cohen's d): {effect_size} (medium effect)")
print(f"Significance level (α): {alpha}")
print(f"Target power (1-β): {power} (80% probability of detecting effect)")
print(f"\nRequired sample size: {n:.0f} per group")
print(f"Total required: {n*2:.0f} people")
 
# Required sample sizes for various effect sizes
print("\nRequired sample sizes by effect size (Power 80%):")
for es, desc in [(0.2, 'small effect'), (0.3, 'medium effect'), (0.5, 'large effect'), (0.8, 'very large effect')]:
    n = power_analysis.solve_power(effect_size=es, alpha=0.05, power=0.8, ratio=1.0)
    print(f"  d = {es} ({desc}): {n:.0f} per group (total {n*2:.0f})")

Q: Comparing means of two groups?
├── Same subjects before/after? → Paired t-test (normal) / Wilcoxon (non-normal)
└── Different subjects? → Independent t-test (normal) / Mann-Whitney (non-normal)

Q: Comparing 3+ groups?
├── Normal distribution? → One-way ANOVA
└── Non-normal distribution? → Kruskal-Wallis

Q: Relationship between two categorical variables?
├── Any expected frequency < 5? → Fisher's exact
└── All expected frequencies >= 5? → Chi-square

Q: Relationship between two continuous variables?
├── Linear relationship? → Pearson r
└── Monotonic relationship? → Spearman ρ