Statistical Analysis Skill

Apply rigorous statistical methods to data questions. Cover descriptive statistics, hypothesis testing, correlation, regression, outlier detection, and trend analysis with proper interpretation.

Workflow

Step 1: Identify the Statistical Question

Map the business question to a statistical framework:

Step 2: Check Assumptions

Before applying any method, verify its assumptions:

import scipy.stats as stats
import numpy as np

# Normality check
stat, p = stats.shapiro(data)
print(f"Shapiro-Wilk: W={stat:.4f}, p={p:.4f}")
print(f"Normal: {'Yes' if p > 0.05 else 'No -- use non-parametric methods'}")

# Variance equality (for t-tests, ANOVA)
stat, p = stats.levene(group_a, group_b)
print(f"Levene's test: F={stat:.4f}, p={p:.4f}")
print(f"Equal variance: {'Yes' if p > 0.05 else 'No -- use Welch t-test'}")

# Sample size check
print(f"Sample sizes: n1={len(group_a)}, n2={len(group_b)}")
print(f"Sufficient for parametric test: {'Yes' if min(len(group_a), len(group_b)) >= 30 else 'Consider non-parametric'}")

Step 3: Apply the Method

Descriptive Statistics

def describe_distribution(series, name="Variable"):
    """Comprehensive descriptive statistics."""
    print(f"\n=== {name} ===")
    print(f"N:        {len(series):,}")
    print(f"Mean:     {series.mean():.4f}")
    print(f"Median:   {series.median():.4f}")
    print(f"Std Dev:  {series.std():.4f}")
    print(f"Skewness: {series.skew():.4f}")
    print(f"Kurtosis: {series.kurtosis():.4f}")
    print(f"IQR:      {series.quantile(0.75) - series.quantile(0.25):.4f}")
    print(f"\nPercentiles:")
    for p in [1, 5, 10, 25, 50, 75, 90, 95, 99]:
        print(f"  P{p:2d}: {series.quantile(p/100):.4f}")

Hypothesis Testing: Two-Sample Comparison

def compare_groups(group_a, group_b, alpha=0.05):
    """Compare two groups with appropriate test selection."""

    # Check normality
    _, p_norm_a = stats.shapiro(group_a[:5000])  # Shapiro limited to 5000
    _, p_norm_b = stats.shapiro(group_b[:5000])
    normal = p_norm_a > 0.05 and p_norm_b > 0.05

    if normal and len(group_a) >= 30 and len(group_b) >= 30:
        # Parametric: Welch's t-test (does not assume equal variance)
        stat, p = stats.ttest_ind(group_a, group_b, equal_var=False)
        test_name = "Welch's t-test"
    else:
        # Non-parametric: Mann-Whitney U
        stat, p = stats.mannwhitneyu(group_a, group_b, alternative='two-sided')
        test_name = "Mann-Whitney U"

    # Effect size (Cohen's d)
    pooled_std = np.sqrt((group_a.std()**2 + group_b.std()**2) / 2)
    cohens_d = (group_a.mean() - group_b.mean()) / pooled_std

    print(f"Test: {test_name}")
    print(f"Statistic: {stat:.4f}")
    print(f"p-value: {p:.6f}")
    print(f"Significant at alpha={alpha}: {'Yes' if p < alpha else 'No'}")
    print(f"Cohen's d: {cohens_d:.4f} ({'small' if abs(cohens_d) < 0.5 else 'medium' if abs(cohens_d) < 0.8 else 'large'})")
    print(f"Mean A: {group_a.mean():.4f}, Mean B: {group_b.mean():.4f}")
    print(f"Difference: {group_a.mean() - group_b.mean():.4f} ({(group_a.mean() - group_b.mean()) / group_b.mean():.2%})")

A/B Test Analysis (Proportions)

def ab_test_proportions(conversions_a, total_a, conversions_b, total_b, alpha=0.05):
    """Analyze an A/B test with conversion rate as the metric."""
    rate_a = conversions_a / total_a
    rate_b = conversions_b / total_b

    # Z-test for proportions
    pooled_p = (conversions_a + conversions_b) / (total_a + total_b)
    se = np.sqrt(pooled_p * (1 - pooled_p) * (1/total_a + 1/total_b))
    z = (rate_b - rate_a) / se
    p_value = 2 * (1 - stats.norm.cdf(abs(z)))

    # Confidence interval for the difference
    se_diff = np.sqrt(rate_a*(1-rate_a)/total_a + rate_b*(1-rate_b)/total_b)
    ci_lower = (rate_b - rate_a) - 1.96 * se_diff
    ci_upper = (rate_b - rate_a) + 1.96 * se_diff

    print(f"Control: {rate_a:.4%} ({conversions_a}/{total_a})")
    print(f"Treatment: {rate_b:.4%} ({conversions_b}/{total_b})")
    print(f"Lift: {(rate_b - rate_a) / rate_a:.2%}")
    print(f"95% CI for difference: [{ci_lower:.4%}, {ci_upper:.4%}]")
    print(f"p-value: {p_value:.6f}")
    print(f"Significant: {'Yes' if p_value < alpha else 'No'}")

Outlier Detection

def detect_outliers(series, method='iqr'):
    """Detect outliers using IQR or Z-score method."""
    if method == 'iqr':
        q1 = series.quantile(0.25)
        q3 = series.quantile(0.75)
        iqr = q3 - q1
        lower = q1 - 1.5 * iqr
        upper = q3 + 1.5 * iqr
        outliers = series[(series < lower) | (series > upper)]
    elif method == 'zscore':
        z_scores = np.abs(stats.zscore(series.dropna()))
        outliers = series[z_scores > 3]

    print(f"Method: {method.upper()}")
    print(f"Outliers found: {len(outliers)} ({len(outliers)/len(series):.2%})")
    if len(outliers) > 0:
        print(f"Range: {outliers.min():.4f} to {outliers.max():.4f}")
    return outliers

Correlation Analysis

def correlation_analysis(df, method='pearson'):
    """Compute and interpret correlations."""
    corr_matrix = df.corr(method=method)

    # Find strong correlations
    pairs = []
    for i in range(len(corr_matrix.columns)):
        for j in range(i+1, len(corr_matrix.columns)):
            r = corr_matrix.iloc[i, j]
            col_a, col_b = corr_matrix.columns[i], corr_matrix.columns[j]
            # Statistical significance
            n = df[[col_a, col_b]].dropna().shape[0]
            t_stat = r * np.sqrt((n - 2) / (1 - r**2))
            p_val = 2 * (1 - stats.t.cdf(abs(t_stat), n - 2))
            pairs.append((col_a, col_b, r, p_val))

    pairs.sort(key=lambda x: abs(x[2]), reverse=True)
    print("Correlations (sorted by strength):")
    for col_a, col_b, r, p in pairs[:15]:
        sig = "*" if p < 0.05 else ""
        print(f"  {col_a} <-> {col_b}: r={r:.3f}{sig}  (p={p:.4f})")

Trend Analysis

from sklearn.linear_model import LinearRegression

def trend_analysis(dates, values):
    """Fit a linear trend and test significance."""
    x = np.arange(len(dates)).reshape(-1, 1)
    y = np.array(values)

    model = LinearRegression().fit(x, y)
    slope = model.coef_[0]
    r_squared = model.score(x, y)

    # Significance of the slope
    n = len(values)
    se = np.sqrt(np.sum((y - model.predict(x))**2) / (n - 2) / np.sum((x - x.mean())**2))
    t_stat = slope / se
    p_value = 2 * (1 - stats.t.cdf(abs(t_stat), n - 2))

    print(f"Trend slope: {slope:.4f} per period")
    print(f"R-squared: {r_squared:.4f}")
    print(f"p-value: {p_value:.6f}")
    print(f"Significant trend: {'Yes' if p_value < 0.05 else 'No'}")

Step 4: Interpret Results

Always provide interpretation in business terms:

## Statistical Summary

**Test**: [Name of test applied]
**Result**: [Statistically significant / not significant] at alpha = 0.05

**What this means**: [Plain language interpretation]
- The observed difference of [X%] is [likely real / could be due to chance].
- We are [95%] confident the true difference is between [CI lower] and [CI upper].
- The effect size is [small/medium/large], meaning [practical significance].

**Caveats**:
- [Assumption violations, if any]
- [Sample size limitations]
- [Multiple comparison concerns]

Output Format

Every statistical analysis should include:

1. The question being answered
2. The method chosen and why
3. Assumption checks and results
4. The test results with exact numbers (test statistic, p-value, CI, effect size)
5. Plain-language interpretation of what the numbers mean
6. Caveats and limitations

Best Practices

1. Report effect sizes, not just p-values. A statistically significant result may not be practically meaningful.
2. Use confidence intervals. They convey more information than binary significant/not-significant.
3. Check assumptions first. Using a parametric test on non-normal data can give misleading results.
4. Correct for multiple comparisons. If running many tests, apply Bonferroni or FDR correction.
5. State the alpha level upfront. Decide on 0.05, 0.01, or 0.10 before seeing the data.
6. Beware of large-sample significance. With millions of rows, even trivial differences become statistically significant. Focus on practical significance.
7. Never confuse correlation with causation. State explicitly that observational data shows association, not causation.
8. Report sample sizes. Always include N for each group. Small samples require different methods.
9. Use non-parametric methods when in doubt. They make fewer assumptions and are more robust.
10. Visualize alongside statistics. A histogram or box plot often reveals more than a p-value.