statistical-analysis

# Statistical Analysis ## Overview Conduct hypothesis tests (t-tests, ANOVA, chi-square), regression, correlation, and Bayesian analyses with systematic assumption checking, effect sizes, and APA-style reporting. Define the estimand and sampling units, diagnose model limitations, and report effect estimates with uncertainty. Numerical screens cannot certify assumptions or scientific validity. ## When to Use This Skill Use this skill when: - Conducting statistical hypothesis tests (t-tests, ANOVA, chi-square, non-parametric) - Performing regression or correlation analyses - Running Bayesian statistical analyses - Checking statistical assumptions and diagnostics - Calculating effect sizes and conducting power analyses - Reporting statistical results in APA format - Analyzing experimental or observational data for research --- ## Installation Targets Pingouin 0.7.0, SciPy 1.18.1, statsmodels 0.15.0, PyMC 6.3.2 and ArviZ 1.3.0 (reviewed 2026-10-01). Create a separate environment; do not add mutually incompatible scientific stacks to a shared project environment. ```bash uv venv --python 3.12 .venv-statistics uv pip install --python .venv-statistics/bin/python "pingouin==0.7.0" "scipy==1.18.1" "statsmodels==0.15.0" pandas matplotlib seaborn # Optional Bayesian examples: uv pip install --python .venv-statistics/bin/python "pymc==6.3.2" "arviz==1.3.0" ``` On Windows use `.venv-statistics/Scripts/python.exe`. Lock the resolved environment for reproduction. - Pingouin 0.7 fixes erroneous mixed/repeated-measures corrections, categorical ANOVA, and other results; do not reproduce an old analysis by merely changing its reported version. See the [official changelog](https://pingouin-stats.org/changelog.html). - Pingouin 0.6+ names are `p_val`, `cohen_d`, `CI95`, `p_unc`. `ttest` reports **absolute** Cohen's d; calculate a signed estimate explicitly. Its `power` column is observed power and should not be used to interpret a null result. - `correction='auto'` chooses Welch by unequal **sample sizes**, not by variance evidence. Set `correction=True` for prespecified Welch inference. - PyMC 6 returns an xarray `DataTree`. ArviZ 1 uses `ci_prob` in `summary`/`plot_dist` and `prob` in `hdi`; specify 95% intervals explicitly. See [Bayesian guidance](references/bayesian_statistics.md). - Pingouin omits `BF10` for one-sided t-tests; posterior directional probabilities are not Bayes factors. A two-sided `BF10` uses its own model/prior and is not the Bayesian counterpart of a Welch unequal-variance model. For model-specific APIs (OLS, GLM, ARIMA), see the **statsmodels** skill. For PyMC workflows, see the **pymc** skill. --- ## Analysis Workflow Use the following workflow, documenting prespecified decisions and exploratory departures. 1. **Frame the question before touching the data.** State the hypothesis, the outcome and predictor variables, and the design (independent vs. paired, number of groups). Specify the target effect, sampling unit, dependence, contrasts and multiplicity family before outcome-driven selection. Label adaptive exploration explicitly. 2. **Inspect the data.** Per group: n, mean, SD, median, missing values. Plot the raw data (histograms or box plots) before any test. Unequal group sizes, missingness, floor/ceiling effects, and outliers all change what test is appropriate — surface them to the user rather than silently working around them. 3. **Select the test** using the quick reference below, or `references/test_selection_guide.md` for designs beyond the basics (counts, time-to-event, reliability, factorial). 4. **Check assumptions** with `scripts/assumption_checks.py`. Interpret plots and screens alongside the design and estimand; do not switch tests automatically at a diagnostic p-value threshold. Report justified sensitivity analyses and any departures from the plan. 5. **Run the test** and always compute the effect size alongside it — a p-value assesses incompatibility with the null model; the effect estimate and its uncertainty support judgments about practical importance. 6. **Report** using the APA templates below, including descriptives, exact statistics, effect sizes with CIs, and the assumption checks performed. If the user only needs one step (e.g., "how many participants do I need?"), jump straight to that section — but still confirm the design assumptions the calculation rests on. --- ## Test Selection Guide ### Quick Reference: Choosing the Right Test Use `references/test_selection_guide.md` for comprehensive guidance (counts, survival, reliability, factorial designs). Quick reference: **Group comparisons:** Independent continuous means usually use prespecified Welch's t-test or Welch's ANOVA. Pooled Student/ANOVA inference needs justified equal variance. Paired means use a paired t-test on aligned within-person differences. Repeated or clustered observations require a model of that dependence. **Rank questions:** Mann-Whitney or Kruskal-Wallis can address distribution/rank contrasts; they are not automatic tests of medians. Wilcoxon signed-rank requires meaningful, symmetrically distributed paired differences. Friedman addresses blocked/repeated rank contrasts. See the reference for ties, small samples, and exact/permutation choices. **Relationships:** Pearson measures linear association; Spearman measures monotone rank association. Choose by the scientific question, not a normality pretest. Use a suitable regression likelihood for continuous, binary, count or time-to-event outcomes. **Bayesian alternatives:** Specify the likelihood, prior and comparison explicitly. Posterior estimation does not automatically yield a Bayes factor for a point null. --- ## Assumption Checking **Always check assumptions before interpreting test results**, and report the checks — reviewers look for them. Use the bundled `scripts/assumption_checks.py` module. Add the skill’s `scripts/` directory to `PYTHONPATH` (the skill root alone is insufficient): ```python from assumption_checks import comprehensive_assumption_check # Within-group outlier screens, normality and variance diagnostics; with plots results = comprehensive_assumption_check( data=df, value_col='score', group_col='group', # Independent groups; paired tests use difference scores alpha=0.05 ) ``` For targeted checks, import individual functions: ```python from assumption_checks import ( check_normality, # Shapiro-Wilk + Q-Q plot + histogram check_normality_per_group, check_homogeneity_of_variance, # Levene's test + box plots check_linearity, # scatter + residual plot for simple regression check_regression_diagnostics, # full OLS diagnostics (see Regression below) detect_outliers # IQR or z-score methods ) result = check_normality(data=df['score'], name='Test Score', alpha=0.05, plot=True) print(result['interpretation']) print(result['recommendation']) ``` ### What to Do When Assumptions Are Violated - A non-significant Shapiro/Levene result does not verify normality/equal variance. There is no universal n=30 robustness rule; tails, skewness, imbalance, dependence and the target estimator matter. - For paired t-tests diagnose differences; for regression inspect residuals. Rank tests change the question and retain assumptions. - Prespecify Welch for independent mean comparisons; use Games-Howell pairwise comparisons when appropriate to a heteroscedastic one-way design. Regression HC3 addresses heteroscedasticity, not clustering, confounding or a wrong mean model. - The helper rejects infinities, undersized or degenerate samples, reports omitted NaNs, retains original outlier positions, and rejects missing group labels. No binary Shapiro verdict is returned above 5000 observations, where its p-value may be inaccurate. - Grouped output uses `outliers_per_group`; legacy `is_normal`/`is_homogeneous` and the `Normal` table label mean only non-rejection. Independence is established through design review, never this helper. Pass `ordered=True` to regression diagnostics only for meaningful row order; Durbin-Watson is descriptive without a universal cutoff. See [assumptions and diagnostics](references/assumptions_and_diagnostics.md). --- ## Running Statistical Tests Primary libraries: - **pingouin**: user-friendly tests that return effect sizes by default — prefer it for standard tests - **scipy.stats**: core statistical tests - **statsmodels**: regression, diagnostics, power analysis - **pymc** + **arviz**: Bayesian modeling and diagnostics ### T-Test with Complete Reporting ```python import pingouin as pg # Prespecified independent mean comparison; arrays must contain finite observations result = pg.ttest(group_a, group_b, correction=True) # Pingouin >= 0.6 column names t_stat = result['T'].values[0] dof = result['dof'].values[0] p_value = result['p_val'].values[0] cohens_d = pg.compute_effsize(group_a, group_b, eftype='cohen') # signed A - B ci_lower, ci_upper = result['CI95'].values[0] # CI for the mean difference p_text = "p < .001" if p_value < .001 else f"p = {p_value:.3f}" print(f"Welch t({dof:.2f}) = {t_stat:.2f}, {p_text}, d = {cohens_d:.2f}") ``` ### ANOVA with Post-Hoc Tests ```python import pingouin as pg aov = pg.anova(dv='score', between='group', data=df, detailed=True) print(aov) # Effect size: partial eta-squared eta_p2 = aov['np2'].values[0] # A prespecified all-pairs family: Tukey HSD controls family-wise error posthoc = pg.pairwise_tukey(dv='score', between='group', data=df) print(posthoc) # Hedges' g, not Cohen's d # Alternative prespecified unequal-variance workflow: welch = pg.welch_anova(dv='score', between='group', data=df) games_howell = pg.pairwise_gameshowell(dv='score', between='group', data=df) ``` ### Linear Regression with Diagnostics ```python import statsmodels.api as sm from assumption_checks import check_regression_diagnostics X = sm.add_constant(X_predictors) # Add intercept model = sm.OLS(y, X).fit() print(model.summary()) # 4-panel residual plot + Shapiro-Wilk, Breusch-Pagan, Durbin-Watson, VIF diag = check_regression_diagnostics(model) print(diag['interpretation']) print(diag['vif']) # If HC3 was specified for independent errors, fit/report that inference explicitly robust = model.get_robustcov_results('HC3') ``` ### Bayesian T-Test ```python import pymc as pm import arviz as az import numpy as np with pm.Model() as model: # Priors mu1 = pm.Normal('mu_group1', mu=0, sigma=10) mu2 = pm.Normal('mu_group2', mu=0, sigma=10) sigma = pm.HalfNormal('sigma', sigma=10) # Likelihood y1 = pm.Normal('y1', mu=mu1, sigma=sigma, observed=group_a) y2 = pm.Normal('y2', mu=mu2, sigma=sigma, observed=group_b) # Derived quantity diff = pm.Deterministic('difference', mu1 - mu2) trace = pm.sample(2000, tune=1000, chains=4, cores=1, random_seed=42, nuts_sampler='pymc') # ArviZ 1.x defaults to 89% intervals; request 95% explicitly for reporting print(az.summary(trace, var_names=['difference'], ci_prob=0.95)) # Direct probability statement (this is what one-sided questions become) prob_greater = np.mean(trace.posterior['difference'].values > 0) print(f"P(mu1 > mu2 | data) = {prob_greater:.3f}") # ArviZ 1.x removed az.plot_posterior; use plot_dist (on 0.x, plot_posterior still works) az.plot_dist(trace, var_names=['difference'], ci_prob=0.95) ``` This is a shared-variance normal model, not Welch's model or a Bayes-factor calculation. Choose priors using meaningful units and external knowledge; inspect prior predictions and prior sensitivity. Short synthetic API smoke tests do not establish convergence; inspect finite R-hat/ESS/MCSE, divergences and posterior predictions before reporting. --- ## Effect Sizes **Effect sizes quantify magnitude; p-values measure incompatibility with a specified null model under its assumptions.** A p-value does not prove that an effect exists or measure its practical importance. Report an effect estimate and uncertainty for every test. See `references/effect_sizes_and_power.md` for the full guide. ### Quick Reference: Common Effect Sizes | Test | Effect Size | Small | Medium | Large | |------|-------------|-------|--------|-------| | T-test | Cohen's d | 0.20 | 0.50 | 0.80 | | ANOVA | η²_p | 0.01 | 0.06 | 0.14 | | Correlation | r | 0.10 | 0.30 | 0.50 | | Regression | R² | 0.02 | 0.13 | 0.26 | | Chi-square (2×2) | Cramér's V | 0.10 | 0.30 | 0.50 | Benchmarks are conventions, not laws — a "small" effect can matter enormously (drug side effects) and a "large" one can be trivial. Interpret in context. ### Calculating Effect Sizes Pingouin returns effect sizes with its tests (`cohen_d` magnitude from `pg.ttest`, `np2` from `pg.anova`, `hedges` from `pg.pairwise_tukey`; `r` from `pg.corr` is already an effect size). ### Confidence Intervals for Effect Sizes Report a CI for the effect size to show its precision; use analyzed sample sizes after the declared missing-data policy. The following is an approximate pooled-d interval, not a heteroscedastic standardized-effect guarantee. Use `pg.compute_esci` (note: `pg.compute_effsize_from_t` returns only the point estimate — it does **not** return a CI): ```python import pingouin as pg d = pg.compute_effsize(group_a, group_b, eftype='cohen') ci_lower, ci_upper = pg.compute_esci(stat=d, nx=len(group_a), ny=len(group_b), eftype='cohen', confidence=0.95) print(f"d = {d:.2f}, 95% CI [{ci_lower:.2f}, {ci_upper:.2f}]") ``` --- ## Power Analysis ### A Priori Power Analysis (Study Planning) Determine required sample size before data collection. These are independent, equal-variance planning models; unequal variances, clustering, attrition or a different estimand require appropriate design-specific calculations or simulation: ```python from statsmodels.stats.power import tt_ind_solve_power, FTestAnovaPower # T-test: What n per group is needed to detect d = 0.5? n_required = tt_ind_solve_power( effect_size=0.5, alpha=0.05, power=0.80, ratio=1.0, alternative='two-sided' ) import math print(f"Required n per group: {math.ceil(n_required)}") # One-way ANOVA: What n is needed to detect Cohen's f = 0.25? # Notes: the parameter is k_groups; effect_size is Cohen's f (f = sqrt(eta2/(1-eta2))); # and solve_power returns the TOTAL sample size, not n per group. import math anova_power = FTestAnovaPower() n_total = anova_power.solve_power( effect_size=0.25, k_groups=3, alpha=0.05, power=0.80 ) per_group = math.ceil(n_total / 3) print(f"Balanced design: {3 * per_group} total ({per_group} per group)") ``` ### Sensitivity Analysis (Post-Study) Determine what effect size the study could detect: ```python # With n=50 per group, what effect could we detect at 80% power? detectable_d = 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"Study could detect d >= {detectable_d:.2f}") ``` **Note**: Post-hoc "observed power" (computing power from the observed effect) is circular and misleading — it is a deterministic function of the p-value. If a study is done and someone asks about power, run a sensitivity analysis instead. See `references/effect_sizes_and_power.md` for detailed guidance. --- ## Reporting Results Follow `references/reporting_standards.md` for APA style. Every report needs: 1. **Descriptive statistics**: M, SD, n for all groups/variables