Data Science Statistical Analysis

Statistical Analysis

Purpose

Perform rigorous statistical analysis: descriptive statistics, hypothesis testing, regression, time series, Bayesian inference, and reproducible reporting.

Agent Protocol

Trigger

Exact user phrases: "statistical analysis", "hypothesis test", "t-test", "ANOVA", "chi-square", "regression", "linear regression", "logistic regression", "time series", "ARIMA", "forecast", "Bayesian", "confidence interval", "p-value", "effect size", "power analysis", "data exploration", "descriptive statistics", "correlation", "outlier detection".

Input Context

• Analysis goal (exploration, inference, prediction, forecasting)
• Data structure (rows, columns, types, missingness)
• Question to answer: what decision depends on this analysis?
• Assumptions about data (independence, distribution, measurement scale)
• Audience: technical or business stakeholders
• Reproducibility requirements (notebook, script, report)

Output Artifact

Statistical analysis results with interpretation, assumptions checked, and business implications.

Response Format

## Analysis Summary
Question: {what was asked}
Data: {N rows, N features, time period}
Method: {test/model, software, parameters}
Results: {estimate, CI, p-value, effect size}
Assumptions Checked: {normality, independence, homoscedasticity}
Interpretation: {what the result means}
Business Implication: {what to do about it}

Completion Criteria

• Data explored and cleaned (missingness, outliers, distributions)
• Appropriate test/model selected based on data and question
• Assumptions checked and violations addressed
• Results reported with effect size and confidence interval
• Interpretation translated to business context
• Reproducible code/notebook provided

Workflow

Step 1: Data Exploration

Summarize data: shape, types, missing values, summary statistics. Visualize: univariate distributions (histograms, box plots), bivariate relationships (scatter plots, correlation heatmaps), multivariate structure (pair plots, PCA). Check data quality: range violations, duplicates, logical inconsistencies, timestamp ordering.

Step 2: Question Mapping

Question Types

Question Refinement

Map vague business question to statistical formulation. "Are we doing better?" → "Is the mean conversion rate in period B higher than period A?" or "Is the trend slope positive and significant?"

Step 3: Model Selection

Method Selection Decision Tree

Data type of outcome (Y)?
├── Continuous (revenue, time, score)
│   ├── Normally distributed → Linear regression, t-test, ANOVA
│   ├── Skewed-positive → Log-transform + linear, or Gamma GLM
│   └── Bounded [0,1] → Beta regression, logit transform
├── Binary (converted, churned, clicked)
│   ├── Logistic regression (logit link)
│   └── Probit regression (normal CDF link)
├── Count (sessions, purchases, errors)
│   ├── Equidispersed → Poisson regression
│   └── Overdispersed → Negative binomial, quasi-Poisson
├── Ordinal (rating 1-5, tier A-B-C)
│   ├── Ordered logit (proportional odds)
│   └── Ordered probit
├── Categorical (brand choice, region)
│   ├── Multinomial logistic regression
│   └── Discriminant analysis
└── Time-to-event (survival, churn time)
    ├── Cox proportional hazards
    └── Kaplan-Meier + log-rank test

Step 4: Assumption Checking

Key Assumptions by Method

Handling Violations

Non-normal data: transform (log, Box-Cox, Yeo-Johnson) or use non-parametric alternative. Heteroscedasticity: use robust standard errors (HC0-HC4, Huber-White). Non-independence: use mixed effects models, clustered standard errors, or GEE. Outliers: robust regression (Huber, quantile), or trimmed estimators.

Step 5: Analysis and Reporting

Complete Report Structure

1. Question: what was asked and why
2. Data description: source, sample size, time period, preprocessing
3. Exploratory analysis: relevant visualizations and summary statistics
4. Model specification: chosen method, formula, assumptions
5. Results: estimate, SE, CI, test statistic, p-value, effect size
6. Diagnostics: assumption checks, sensitivity analysis
7. Interpretation: what the result means in business terms
8. Limitations: what could affect validity of conclusions

Reporting Best Practices

Round appropriately (2-3 significant digits for estimates). Use consistent terminology. State both statistical and practical significance. Include visualizations for key findings. Report uncertainty (CIs, prediction intervals) not just point estimates.

Step 6: Effect Sizes

Step 7: Reproducibility

Workflow Requirements

Script everything: analysis is code, not manual operations. Version control: analysis code + processed data. Dependency management: pinned package versions. Seed random generators for replicable sampling. Output pinning: save processed data and results.

Documentation

Document data sources and extraction queries. Document inclusion/exclusion criteria and filtering steps. Describe handling of missing data, outliers, and transformations. Justify modeling decisions.

Decision Trees

Hypothesis Test Selection

How many groups?
├── One group
│   ├── Compare to known value → One-sample t-test
│   └── Compare to distribution → Kolmogorov-Smirnov
├── Two groups
│   ├── Independent → Two-sample t-test or Mann-Whitney
│   ├── Paired/matched → Paired t-test or Wilcoxon signed-rank
│   └── Time-to-event → Log-rank test
└── Three+ groups
    ├── Independent → ANOVA or Kruskal-Wallis
    ├── Repeated measures → Repeated measures ANOVA or Friedman
    └── Time-to-event → Cox regression

Regression Type Selection

Outcome type?
├── Continuous
│   ├── Simple linear → Linear regression (OLS)
│   ├── With mixed effects → Mixed linear model (LMM)
│   └── Quantile of interest → Quantile regression
├── Binary → Logistic regression
├── Count → Poisson or negative binomial
├── Time-to-event → Cox proportional hazards
└── Multivariate → MANOVA or multivariate regression

Python Implementation Examples

Hypothesis Tests

import numpy as np
import pandas as pd
from scipy import stats

# Two-sample t-test (continuous, normal)
def two_sample_test(treatment, control, equal_var=True):
    t_stat, p_value = stats.ttest_ind(treatment, control, equal_var=equal_var)
    d = (np.mean(treatment) - np.mean(control)) / np.sqrt(
        (np.var(treatment, ddof=1) + np.var(control, ddof=1)) / 2
    )
    return {"statistic": t_stat, "p_value": p_value, "effect_size": d}

# Chi-square test (categorical)
def chi_square_test(observed, expected=None):
    if expected is None:
        chi2, p_value, dof, expected = stats.chi2_contingency(observed)
    else:
        chi2, p_value = stats.chisquare(observed, expected)
    # Cramér's V
    n = np.sum(observed)
    k = min(observed.shape)
    v = np.sqrt(chi2 / (n * (k - 1))) if n * (k - 1) > 0 else 0
    return {"statistic": chi2, "p_value": p_value, "effect_size": v}

# Mann-Whitney U (non-parametric two-group)
def mann_whitney_test(treatment, control):
    u_stat, p_value = stats.mannwhitneyu(treatment, control, alternative='two-sided')
    # Rank-biserial correlation as effect size
    n1, n2 = len(treatment), len(control)
    r = 1 - (2 * u_stat) / (n1 * n2)
    return {"statistic": u_stat, "p_value": p_value, "effect_size": r}

# ANOVA
def anova_test(*groups):
    f_stat, p_value = stats.f_oneway(*groups)
    # Eta-squared
    all_data = np.concatenate(groups)
    grand_mean = np.mean(all_data)
    ss_between = sum(len(g) * (np.mean(g) - grand_mean)**2 for g in groups)
    ss_total = sum((x - grand_mean)**2 for x in all_data)
    eta_sq = ss_between / ss_total
    return {"statistic": f_stat, "p_value": p_value, "effect_size": eta_sq}

Linear Regression with Diagnostics

import statsmodels.api as sm

def regression_with_diagnostics(X, y):
    X = sm.add_constant(X)
    model = sm.OLS(y, X).fit()
    
    # Diagnostics
    residuals = model.resid
    fitted = model.fittedvalues
    
    # Homoscedasticity: Breusch-Pagan test
    bp_test = sm.stats.diagnostic.het_breuschpagan(residuals, X)
    
    # Normality: Jarque-Bera test
    jb_test = sm.stats.jarque_bera(residuals)
    
    # Multicollinearity: VIF
    from statsmodels.stats.outliers_influence import variance_inflation_factor
    vif = [variance_inflation_factor(X.values, i) for i in range(X.shape[1])]
    
    # Durbin-Watson (autocorrelation)
    dw = sm.stats.durbin_watson(residuals)
    
    return {
        "model": model,
        "r_squared": model.rsquared,
        "adj_r_squared": model.rsquared_adj,
        "coeff": model.params,
        "p_values": model.pvalues,
        "dw_statistic": dw,
        "bp_pvalue": bp_test[1],
        "jb_pvalue": jb_test[1],
        "max_vif": max(vif),
        "significant": model.pvalues.iloc[1:].lt(0.05).any()
    }

Time Series Analysis

Components and Decomposition

from statsmodels.tsa.seasonal import seasonal_decompose
from statsmodels.tsa.stattools import adfuller, acf, pacf

# Decompose time series into trend + seasonal + residual
def decompose_series(series, period=7):
    result = seasonal_decompose(series, model='additive', period=period)
    return {
        "trend": result.trend,
        "seasonal": result.seasonal,
        "residual": result.resid,
        "strength_trend": 1 - np.var(result.resid.dropna()) / np.var(result.trend.dropna()),
        "strength_seasonality": 1 - np.var(result.resid.dropna()) / np.var(result.seasonal.dropna())
    }

# Stationarity test: ADF
def check_stationarity(series):
    result = adfuller(series.dropna(), autolag='AIC')
    return {
        "adf_statistic": result[0],
        "p_value": result[1],
        "critical_values": result[4],
        "is_stationary": result[1] < 0.05
    }

ARIMA Model Selection

arima_workflow:
  step_1_check_stationarity:
    - "Run ADF test on series"
    - "If non-stationary: difference once, re-test"
    - "If still non-stationary: log transform + difference"
    - "d = number of differences needed"
  
  step_2_identify_p_and_q:
    - "Plot ACF: significant spikes at lags → p component"
    - "Plot PACF: significant spikes at lags → q component"
    - "If ACF decays gradually, PACF cuts off after p → AR(p)"
    - "If PACF decays gradually, ACF cuts off after q → MA(q)"
    - "Seasonal: use seasonal differencing, inspect seasonal lags"
  
  step_3_fit_and_evaluate:
    - "Fit candidate models (p-1, p, p+1) × (q-1, q, q+1)"
    - "Compare AIC / BIC — lower is better"
    - "Check residuals: should be white noise (Ljung-Box test)"
    - "Residuals should have no significant ACF spikes"
  
  step_4_forecast:
    - "Generate point forecasts with prediction intervals"
    - "Evaluate out-of-sample: rolling window CV or holdout"
    - "Metrics: MAE, RMSE, MAPE, MASE"

from statsmodels.tsa.arima.model import ARIMA
from statsmodels.stats.diagnostic import acorr_ljungbox

def fit_and_diagnose_arima(series, order, seasonal_order=None):
    model = ARIMA(series, order=order, seasonal_order=seasonal_order)
    fitted = model.fit()
    
    # Residual diagnostics
    residuals = fitted.resid
    lb_test = acorr_ljungbox(residuals, lags=[10, 20, 30], return_df=True)
    
    return {
        "model": fitted,
        "aic": fitted.aic,
        "bic": fitted.bic,
        "residuals_white_noise": lb_test['lb_pvalue'].gt(0.05).all(),
        "forecast": fitted.get_forecast(steps=30)
    }

Bayesian Inference

Bayesian A/B Test with PyMC

import pymc as pm

def bayesian_ab_test(clicks_a, impressions_a, clicks_b, impressions_b):
    with pm.Model() as ab_model:
        # Priors
        p_a = pm.Beta("p_a", alpha=1, beta=1)
        p_b = pm.Beta("p_b", alpha=1, beta=1)
        
        # Likelihood
        obs_a = pm.Binomial("obs_a", n=impressions_a, p=p_a, observed=clicks_a)
        obs_b = pm.Binomial("obs_b", n=impressions_b, p=p_b, observed=clicks_b)
        
        # Derived quantities
        lift = pm.Deterministic("lift", (p_b - p_a) / p_a)
        effect_greater_0 = pm.Deterministic("prob_positive", pm.math.gt(lift, 0))
        
        # Sample
        trace = pm.sample(2000, tune=1000, chains=4, random_seed=42)
    
    return {
        "p_a_posterior": trace.posterior.p_a.values.flatten(),
        "p_b_posterior": trace.posterior.p_b.values.flatten(),
        "lift_posterior": trace.posterior.lift.values.flatten(),
        "prob_better_than_a": np.mean(trace.posterior.lift.values > 0),
        "hdi_lift": pm.hdi(trace.posterior.lift, hdi_prob=0.94)
    }

Bayesian Linear Regression

def bayesian_regression(X, y):
    with pm.Model() as reg_model:
        # Priors
        alpha = pm.Normal("alpha", mu=0, sigma=10)
        beta = pm.Normal("beta", mu=0, sigma=5, shape=X.shape[1])
        sigma = pm.HalfNormal("sigma", sigma=5)
        
        # Linear model
        mu = alpha + pm.math.dot(X, beta)
        
        # Likelihood
        likelihood = pm.Normal("y", mu=mu, sigma=sigma, observed=y)
        
        # Sample
        trace = pm.sample(2000, tune=1000, chains=4, random_seed=42)
    
    return {
        "alpha_samples": trace.posterior.alpha.values.flatten(),
        "beta_samples": trace.posterior.beta.values,
        "sigma_samples": trace.posterior.sigma.values.flatten(),
        "r2_samples": pm.r2_score(y, trace.posterior.alpha + 
                                   pm.math.dot(X, trace.posterior.beta.T)).values
    }

Dimensionality Reduction and Clustering

PCA Implementation

from sklearn.decomposition import PCA
from sklearn.preprocessing import StandardScaler

def pca_analysis(X, n_components=None, standardize=True):
    if standardize:
        X_scaled = StandardScaler().fit_transform(X)
    else:
        X_scaled = X
    
    pca = PCA(n_components=n_components or min(X.shape[0], X.shape[1], 50))
    X_pca = pca.fit_transform(X_scaled)
    
    # Scree plot data
    explained_var = pca.explained_variance_ratio_
    cumulative_var = np.cumsum(explained_var)
    n_for_90pct = np.argmax(cumulative_var >= 0.90) + 1
    n_for_95pct = np.argmax(cumulative_var >= 0.95) + 1
    
    # Feature loadings
    loadings = pd.DataFrame(
        pca.components_.T,
        index=X.columns,
        columns=[f"PC{i+1}" for i in range(pca.n_components_)]
    )
    
    return {
        "explained_variance_ratio": explained_var,
        "cumulative_variance": cumulative_var,
        "n_components_90pct": n_for_90pct,
        "n_components_95pct": n_for_95pct,
        "loadings": loadings,
        "transformed": X_pca
    }

K-Means with Optimal K Selection

from sklearn.cluster import KMeans
from sklearn.metrics import silhouette_score

def optimal_kmeans(X, k_range=(2, 15), random_state=42):
    inertias = []
    silhouettes = []
    
    for k in range(k_range[0], k_range[1] + 1):
        kmeans = KMeans(n_clusters=k, n_init=10, random_state=random_state)
        labels = kmeans.fit_predict(X)
        inertias.append(kmeans.inertia_)
        silhouettes.append(silhouette_score(X, labels))
    
    # Elbow point (rate of change)
    deltas = np.diff(inertias)
    delta_deltas = np.diff(deltas)
    elbow_k = np.argmax(delta_deltas) + k_range[0] + 1
    
    # Best silhouette
    best_k = np.argmax(silhouettes) + k_range[0]
    
    return {"inertias": inertias, "silhouettes": silhouettes, 
            "elbow_k": elbow_k, "best_silhouette_k": best_k}

Missing Data Handling

Missing Data Mechanisms

MICE (Multiple Imputation by Chained Equations)

from sklearn.experimental import enable_iterative_imputer
from sklearn.impute import IterativeImputer
from sklearn.ensemble import RandomForestRegressor

def mice_imputation(df, n_imputations=5, max_iter=10):
    imputer = IterativeImputer(
        estimator=RandomForestRegressor(n_estimators=100, n_jobs=-1),
        max_iter=max_iter,
        random_state=42,
        sample_posterior=True,  # Add noise for proper multiple imputation
        initial_strategy='median'
    )
    
    # Generate multiple imputed datasets
    imputed_datasets = []
    for i in range(n_imputations):
        imputed = imputer.fit_transform(df)
        imputed_datasets.append(pd.DataFrame(imputed, columns=df.columns))
    
    return imputed_datasets

# Pooling results (Rubin's rules)
def pool_estimates(estimates, variances):
    # Within-imputation variance
    m = len(estimates)
    w = np.mean(variances)
    # Between-imputation variance
    b = np.var(estimates, ddof=1)
    # Total variance
    total_var = w + (1 + 1/m) * b
    return np.mean(estimates), np.sqrt(total_var)

Bootstrapping

def bootstrap_ci(data, statistic=np.mean, n_bootstrap=10000, ci_level=0.95):
    n = len(data); stats = np.zeros(n_bootstrap)
    for i in range(n_bootstrap):
        stats[i] = statistic(np.random.choice(data, size=n, replace=True))
    alpha = 1 - ci_level
    return {"statistic": statistic(data),
            "ci": (np.percentile(stats, 100*alpha/2), np.percentile(stats, 100*(1-alpha/2))),
            "se": np.std(stats)}

Decision Trees (continued)

Missing Data Strategy

Missing data mechanism?
├── MCAR (< 5% missing)
│   └── Listwise deletion (simple, unbiased)
├── MCAR (5-15% missing)
│   ├── Numeric → Median imputation
│   └── Categorical → Mode imputation
├── MAR
│   ├── < 10% variables missing → MICE (5-10 imputations)
│   ├── Time series → Forward fill / interpolation
│   └── ML prediction → Regression imputation
└── MNAR
    ├── Sensitivity analysis (bounds for missing values)
    └── Pattern-mixture models (stratify by missing pattern)

Time Series Forecasting Method

Data characteristics?
├── Strong trend, no seasonality
│   ├── Linear trend → Holt's linear (double exponential smoothing)
│   └── Exponential trend → Holt-Winters (triple, no season)
├── Trend + seasonality
│   ├── Stable seasonality → Holt-Winters additive
│   └── Growing seasonality → Holt-Winters multiplicative
├── No clear pattern, many data points
│   ├── Univariate → ARIMA/SARIMA (auto-arima for selection)
│   └── Multivariate → VAR, Prophet, or LSTM
└── Irregular events, holidays
    └── Prophet (handles holidays, changepoints, outliers)

Rules

• Explore data before modeling — summary stats + visualizations first
• Check assumptions before choosing a test — wrong test gives wrong answer
• Report effect sizes and confidence intervals, not just p-values
• Address missing data explicitly — document method and rationale
• Use robust methods when assumptions are violated (robust SEs, non-parametric)
• Script everything for reproducibility — no manual steps
• Visualize results for stakeholders — tables are for reference, charts tell the story
• Document data inclusion/exclusion criteria
• Run sensitivity analyses to test robustness of conclusions
• Limit conclusions to what the data and method can support
• Validate stationarity before fitting time series models
• Use ACF/PACF for ARIMA order selection — avoid p-hacking orders
• Apply Rubin's rules for pooling multiple imputation results
• Prefer cross-validation over single train/test split for model evaluation
• Use bootstrap for non-standard metrics or small samples
• Diagnose residuals — a good model has white-noise residuals

References

• references/statistical-analysis-fundamentals.md — Statistical Analysis Fundamentals
• references/statistical-analysis-advanced.md — Statistical Analysis Advanced Topics