Data Analysis

Python for Statistics: From Clean Data to Defensible Inference

GOSPELTRADER Research Desk · 11 September 2026 · 12 min read

Quick answer

For statistical inference in Python use pandas to clean, scipy.stats for hypothesis tests, and statsmodels for regression — statsmodels reports coefficients, standard errors, p-values, confidence intervals and diagnostics, which scikit-learn deliberately does not. Use scikit-learn only when prediction, not inference, is the goal.

Pick the right library for the question

The most common mistake in Python analysis is fitting a scikit-learn LinearRegression and then trying to report p-values that the object never produced. Inference and prediction are different tasks with different tools.

Goal	Library	Why
Data cleaning and reshaping	pandas	Vectorised joins, groupby, missing handling
Hypothesis tests	scipy.stats	t-tests, ANOVA, chi-square, non-parametrics
Regression with inference	statsmodels	SE, p-values, CIs, diagnostics, formula API
Predictive modelling	scikit-learn	Cross-validation, pipelines, regularisation
Time series	statsmodels / pmdarima	ARIMA, SARIMAX, stationarity tests
Visualisation	matplotlib / seaborn	Diagnostic and publication plots

Clean the data with an auditable pipeline

Write cleaning as chained, named steps rather than scattered edits, so the pipeline can be reviewed and re-run.

import pandas as pd, numpy as np

df = pd.read_csv("survey.csv")
print(df.info()); print(df.isna().sum()); print(df.duplicated().sum())

df = (df
      .drop_duplicates()
      .assign(q3r=lambda d: 6 - d["q3"],                       # reverse score
              gender=lambda d: d["gender"].map({1: "Male", 2: "Female"}))
      .query("age >= 18"))

df["satisfaction"] = df[["q1", "q2", "q3r", "q4"]].mean(axis=1)

# Outlier screen (z > 3.29)
z = (df["satisfaction"] - df["satisfaction"].mean()) / df["satisfaction"].std()
df["outlier"] = z.abs() > 3.29

Hypothesis testing with scipy

Always test the assumption first, then choose the parametric or non-parametric route, and always compute an effect size.

from scipy import stats

a = df.loc[df.gender == "Male", "satisfaction"]
b = df.loc[df.gender == "Female", "satisfaction"]

stats.shapiro(a); stats.shapiro(b)          # normality
stats.levene(a, b)                          # equal variance

t, p = stats.ttest_ind(a, b, equal_var=False)   # Welch's t-test

# Cohen's d (pooled SD)
n1, n2 = len(a), len(b)
sp = np.sqrt(((n1-1)*a.var(ddof=1) + (n2-1)*b.var(ddof=1)) / (n1+n2-2))
d = (a.mean() - b.mean()) / sp
print(f"t = {t:.2f}, p = {p:.3f}, d = {d:.2f}")

Regression with statsmodels, including diagnostics

The formula API mirrors R's syntax and gives a full inferential summary. Categorical predictors are dummy-coded automatically with C().

import statsmodels.formula.api as smf
import statsmodels.api as sm
from statsmodels.stats.outliers_influence import variance_inflation_factor

model = smf.ols("productivity ~ training_hours + experience + C(department)",
                data=df).fit()
print(model.summary())          # R², F, coefficients, p, 95% CI

# Heteroskedasticity-robust standard errors
robust = model.get_robustcov_results(cov_type="HC3")

# Assumption checks
sm.stats.diagnostic.het_breuschpagan(model.resid, model.model.exog)  # constant variance
stats.shapiro(model.resid)                                            # residual normality
sm.stats.durbin_watson(model.resid)                                   # independence

X = sm.add_constant(df[["training_hours", "experience"]])
[variance_inflation_factor(X.values, i) for i in range(X.shape[1])]   # VIF < 5

Key formulas behind the output

Knowing what the summary table computes prevents misreporting.

OLS estimator      : β̂ = (XᵀX)⁻¹ Xᵀ y
Coefficient of det.: R² = 1 - SSE/SST,  Adj. R² = 1 - (1-R²)(n-1)/(n-k-1)
t-statistic        : t = β̂ / SE(β̂),   df = n - k - 1
F-statistic        : F = (SSR/k) / (SSE/(n-k-1))
Cohen's d          : d = (M₁ - M₂) / s_pooled
Variance inflation : VIF_j = 1 / (1 - R²_j)

When to switch to scikit-learn

If the deliverable is an accurate prediction rather than a tested hypothesis — churn scoring, demand forecasting, credit risk — use scikit-learn with a train/test split, cross-validation and a pipeline that fits preprocessing inside each fold. Report RMSE, MAE, ROC-AUC or F1 instead of p-values, and never evaluate on data used to fit.

Frequently asked questions

Should I use statsmodels or scikit-learn for regression?

Use statsmodels when you need to interpret coefficients and report significance, which is the case for almost all academic and policy research. Use scikit-learn when the objective is out-of-sample predictive accuracy and you need cross-validation, regularisation and pipelines.

Is Python acceptable for a thesis or dissertation?

Yes, provided you report the same elements a reviewer expects: sample, assumptions tested, test statistic, degrees of freedom, exact p-value and effect size. Include the notebook or script as an appendix so the analysis is reproducible.

How do I handle missing data in pandas?

Diagnose first: report how much is missing and whether it appears random. Listwise deletion is defensible below roughly 5% missing on a variable. Above that, use multiple imputation (statsmodels MICE or sklearn IterativeImputer) and state the method and number of imputations in the write-up. Mean imputation understates variance and should be avoided in inferential work.

Which Python version and packages should a research setup use?

Python 3.11 or newer with pandas, numpy, scipy, statsmodels, matplotlib, seaborn and scikit-learn, installed in a per-project virtual environment with a pinned requirements.txt so the analysis can be reproduced exactly.

Want this analysis done for you?

We deliver cleaning, assumption testing, modelling and a reporting-ready write-up.

Start a project