skills/statistical-analysis/references/effect_sizes_and_power.md
This document provides guidance on calculating, interpreting, and reporting effect sizes, as well as conducting power analyses for study planning.
Golden rule: ALWAYS report effect sizes alongside p-values.
Formula:
Interpretation (Cohen, 1988):
Context-dependent interpretation:
Python calculation:
import pingouin as pg
import numpy as np
# Independent t-test with effect size
result = pg.ttest(group1, group2, correction=False)
cohens_d = result['cohen_d'].values[0]
# (pingouin 0.6.0 renamed columns; on 0.5.x use 'p-val', 'cohen-d', 'CI95%', 'p-unc')
# Manual calculation
mean_diff = np.mean(group1) - np.mean(group2)
pooled_std = np.sqrt((np.var(group1, ddof=1) + np.var(group2, ddof=1)) / 2)
cohens_d = mean_diff / pooled_std
# Paired t-test
result = pg.ttest(pre, post, paired=True)
cohens_d = result['cohen_d'].values[0]
Confidence intervals for d:
import pingouin as pg
# compute_effsize_from_t returns only the point estimate;
# get the CI separately with compute_esci
d = pg.compute_effsize(group1, group2, eftype='cohen')
ci = pg.compute_esci(stat=d, nx=len(group1), ny=len(group2),
eftype='cohen', confidence=0.95)
Why use it: Cohen's d has slight upward bias with small samples (n < 20)
Formula: g = d × correction_factor, where correction_factor = 1 - 3/(4df - 1)
Python calculation:
# pg.ttest output has no Hedges' g column; compute it directly
hedges_g = pg.compute_effsize(group1, group2, eftype='hedges')
Use Hedges' g when:
When to use: When one group is a control with known variability
Formula: Δ = (M₁ - M₂) / SD_control
Use cases:
What it measures: Proportion of total variance explained by factor
Formula: η² = SS_effect / SS_total
Interpretation:
Limitation: In multi-factor designs each effect's η² shrinks as other factors are added (classical η² values sum to ≤ 1.0 by construction); it is partial η² that can sum to > 1.0 across factors
Python calculation:
import pingouin as pg
# One-way ANOVA (detailed=True is required for the SS column)
aov = pg.anova(dv='value', between='group', data=df, detailed=True)
eta_squared = aov['SS'][0] / aov['SS'].sum()
# Or read pingouin's np2 column, which is PARTIAL eta-squared:
partial_eta_sq = aov['np2'][0]
# np2 coincides with classical eta-squared only for one-way (single-factor) designs
What it measures: Proportion of variance explained by factor, excluding other factors
Formula: η²_p = SS_effect / (SS_effect + SS_error)
Interpretation: Same benchmarks as η²
When to use: Multi-factor ANOVA (standard in factorial designs)
Limitation: Across factors, partial η² values can sum to > 1.0 — they are not additive shares of total variance
Python calculation:
aov = pg.anova(dv='value', between=['factor1', 'factor2'], data=df)
# pingouin reports partial eta-squared by default
partial_eta_sq = aov['np2']
What it measures: Less biased estimate of population variance explained
Why use it: η² overestimates effect size; ω² provides better population estimate
Formula: ω² = (SS_effect - df_effect × MS_error) / (SS_total + MS_error)
Interpretation: Same benchmarks as η², but typically smaller values
Python calculation:
def omega_squared(aov_table):
ss_effect = aov_table.loc[0, 'SS']
ss_total = aov_table['SS'].sum()
ms_error = aov_table.loc[aov_table.index[-1], 'MS'] # Residual MS
df_effect = aov_table.loc[0, 'DF']
omega_sq = (ss_effect - df_effect * ms_error) / (ss_total + ms_error)
return omega_sq
What it measures: Effect size for ANOVA (analogous to Cohen's d)
Formula: f = √(η² / (1 - η²))
Interpretation:
Python calculation:
eta_squared = 0.06 # From ANOVA
cohens_f = np.sqrt(eta_squared / (1 - eta_squared))
Use in power analysis: Required for ANOVA power calculations
Interpretation:
Important notes:
Python calculation:
import pingouin as pg
# Pearson correlation with CI
result = pg.corr(x, y, method='pearson')
r = result['r'].values[0]
ci = result['CI95'].values[0] # pingouin 0.6.0 renamed CI95% to CI95
# Spearman correlation
result = pg.corr(x, y, method='spearman')
rho = result['r'].values[0]
What it measures: Proportion of variance in Y explained by model
Interpretation:
Context-dependent:
Python calculation:
from sklearn.metrics import r2_score
import statsmodels.api as sm
# Using statsmodels (add_constant adds the intercept column)
model = sm.OLS(y, sm.add_constant(X)).fit()
r_squared = model.rsquared
adjusted_r_squared = model.rsquared_adj
# Manual
r_squared = 1 - (SS_residual / SS_total)
Why use it: R² artificially increases when adding predictors; adjusted R² penalizes model complexity
Formula: R²_adj = 1 - (1 - R²) × (n - 1) / (n - k - 1)
When to use: Always report alongside R² for multiple regression
What it measures: Effect of one-SD change in predictor on outcome (in SD units)
Interpretation: Similar to Cohen's d
Python calculation:
from scipy import stats
# Standardize variables first
X_std = (X - X.mean()) / X.std()
y_std = (y - y.mean()) / y.std()
model = OLS(y_std, X_std).fit()
beta = model.params
What it measures: Effect size for individual predictors or model comparison
Formula: f² = (R²_AB - R²_A) / (1 - R²_AB)
Where:
Interpretation:
Python calculation:
# Compare two nested models
model_full = OLS(y, X_full).fit()
model_reduced = OLS(y, X_reduced).fit()
r2_full = model_full.rsquared
r2_reduced = model_reduced.rsquared
f_squared = (r2_full - r2_reduced) / (1 - r2_full)
What it measures: Association strength for χ² test (works for any table size)
Formula: V = √(χ² / (n × (k - 1)))
Where k = min(rows, columns)
Interpretation (benchmarks depend on df* = min(rows, columns) − 1):
| df* | Small | Medium | Large |
|---|---|---|---|
| 1 (2×2) | 0.10 | 0.30 | 0.50 |
| 2 | 0.07 | 0.21 | 0.35 |
| 3 | 0.06 | 0.17 | 0.29 |
For 2×2 tables: Use phi coefficient (φ)
Python calculation:
import numpy as np
from scipy.stats.contingency import association
# Cramér's V
cramers_v = association(contingency_table, method='cramer')
# Phi coefficient (2x2): |phi| equals Cramér's V for a 2x2 table.
# Caution: method='pearson' is Pearson's contingency coefficient, NOT phi.
a, b, c, d = np.asarray(contingency_table).ravel()
phi = (a * d - b * c) / np.sqrt((a + b) * (c + d) * (a + c) * (b + d)) # signed phi
For 2×2 contingency tables:
| Outcome + | Outcome - | |
|---|---|---|
| Exposed | a | b |
| Unexposed | c | d |
Odds Ratio: OR = (a/b) / (c/d) = ad / bc
Interpretation:
Risk Ratio: RR = (a/(a+b)) / (c/(c+d))
When to use:
Python calculation:
import numpy as np
from scipy import stats
import statsmodels.api as sm
# From contingency table
odds_ratio = (a * d) / (b * c)
# Fisher's exact test returns only the sample OR and a p-value (no CI)
table = np.array([[a, b], [c, d]])
oddsratio, pvalue = stats.fisher_exact(table)
# Odds-ratio confidence interval (scipy >= 1.10; conditional MLE estimate)
or_result = stats.contingency.odds_ratio(table)
ci = or_result.confidence_interval(confidence_level=0.95)
# From logistic regression
model = sm.Logit(y, X).fit()
odds_ratios = np.exp(model.params) # Exponentiate coefficients
ci = np.exp(model.conf_int()) # Exponentiate CIs
Rank-biserial correlation (r_rb): Effect size for Mann-Whitney U and Wilcoxon signed-rank tests (range −1 to 1; interpret |r_rb| roughly like r). Returned by pg.mwu and pg.wilcoxon as the RBC column.
Common-language effect size (CLES): Probability that a randomly sampled value from one group exceeds a randomly sampled value from the other (0.5 = no effect). Returned by pg.mwu as CLES.
r = z / √N: Classic effect size when a z approximation is reported for Mann-Whitney/Wilcoxon (small 0.10, medium 0.30, large 0.50).
Epsilon-squared (ε²): Effect size for Kruskal-Wallis: ε² = H × (n + 1) / (n² − 1).
Python calculation:
import numpy as np
import pandas as pd
import pingouin as pg
res = pg.mwu(group1, group2)
print(res[['U_val', 'p_val', 'RBC', 'CLES']])
# Kruskal-Wallis with epsilon-squared
df = pd.DataFrame({'value': np.concatenate([group1, group2, group3]),
'group': np.repeat(['a', 'b', 'c'],
[len(group1), len(group2), len(group3)])})
kw = pg.kruskal(df, dv='value', between='group')
H, n = kw['H'].values[0], len(df)
epsilon_sq = H * (n + 1) / (n**2 - 1)
What it measures: Ratio of evidence for alternative vs. null hypothesis
Interpretation:
For the full Jeffreys interpretation table and BF reporting language, see bayesian_statistics.md.
Python calculation:
import pingouin as pg
# Pingouin 0.5+: two-sided BF10 on independent t-tests; use BayesFactor/JASP/PyMC for full inference
result = pg.ttest(group1, group2, correction=False)
bf10 = result['BF10'].values[0]
When no analytic CI exists for an effect size (or its assumptions are doubtful), bootstrap one:
import numpy as np
from scipy import stats
def cohen_d(x, y):
nx, ny = len(x), len(y)
sp = np.sqrt(((nx - 1) * np.var(x, ddof=1) + (ny - 1) * np.var(y, ddof=1))
/ (nx + ny - 2))
return (np.mean(x) - np.mean(y)) / sp
boot = stats.bootstrap((group1, group2), cohen_d, n_resamples=9999,
method='BCa', rng=np.random.default_rng(42))
print(boot.confidence_interval) # 95% BCa CI by default
The same pattern works for any statistic (medians, correlations, rank-biserial, ...). Prefer method='BCa' and use at least 5000-10000 resamples.
Statistical power: Probability of detecting an effect if it exists (1 - β)
Conventional standards:
Four interconnected parameters (given 3, can solve for 4th):
Purpose: Determine required sample size before study
Steps:
Python implementation:
from statsmodels.stats.power import (
tt_ind_solve_power,
zt_ind_solve_power,
FTestAnovaPower,
NormalIndPower
)
# T-test power analysis
n_required = tt_ind_solve_power(
effect_size=0.5, # Cohen's d
alpha=0.05,
power=0.80,
ratio=1.0, # Equal group sizes
alternative='two-sided'
)
# ANOVA power analysis (kwarg is k_groups, not ngroups)
anova_power = FTestAnovaPower()
n_total = anova_power.solve_power(
effect_size=0.25, # Cohen's f
k_groups=3,
alpha=0.05,
power=0.80
)
# Returns the TOTAL sample size across all groups:
# f = 0.25, k = 3 -> ~158 total, i.e. ~53 per group
# Correlation power analysis
from pingouin import power_corr
n_required = power_corr(r=0.30, power=0.80, alpha=0.05)
⚠️ CAUTION: Post hoc power is controversial and often not recommended
Why it's problematic:
When it might be acceptable:
Better alternatives:
Purpose: Determine minimum detectable effect size given study parameters
When to use: After study is complete, to understand study's capability
Python implementation:
# What effect size could we detect with n=50 per group?
detectable_effect = tt_ind_solve_power(
effect_size=None, # Solve for this
nobs1=50,
alpha=0.05,
power=0.80,
ratio=1.0,
alternative='two-sided'
)
print(f"With n=50 per group, we could detect d ≥ {detectable_effect:.2f}")
T-test example:
"Group A (M = 75.2, SD = 8.5) scored significantly higher than Group B (M = 68.3, SD = 9.2), t(98) = 3.82, p < .001, d = 0.77, 95% CI [0.36, 1.18]."
ANOVA example:
"There was a significant main effect of treatment condition on test scores, F(2, 87) = 8.45, p < .001, η²p = .16. Post hoc comparisons using Tukey's HSD revealed..."
Correlation example:
"There was a moderate positive correlation between study time and exam scores, r(148) = .42, p < .001, 95% CI [.27, .55]."
Regression example:
"The regression model significantly predicted exam scores, F(3, 146) = 45.2, p < .001, R² = .48. Study hours (β = .52, p < .001) and prior GPA (β = .31, p < .001) were significant predictors."
Bayesian example: See bayesian_statistics.md (Reporting Bayesian Results) for Bayes Factor and posterior reporting templates.
| Analysis | Effect Size | Small | Medium | Large |
|---|---|---|---|---|
| T-test | Cohen's d | 0.20 | 0.50 | 0.80 |
| ANOVA | η², ω² | 0.01 | 0.06 | 0.14 |
| ANOVA | Cohen's f | 0.10 | 0.25 | 0.40 |
| Correlation | r, ρ | 0.10 | 0.30 | 0.50 |
| Regression | R² | 0.02 | 0.13 | 0.26 |
| Regression | f² | 0.02 | 0.15 | 0.35 |
| Chi-square (df* = 1, 2×2) | Cramér's V, φ | 0.10 | 0.30 | 0.50 |
| Chi-square (df* = 2) | Cramér's V | 0.07 | 0.21 | 0.35 |
| Chi-square (df* = 3) | Cramér's V | 0.06 | 0.17 | 0.29 |
Note: df* = min(rows, columns) − 1.