---
title: "Statistical Analysis for Imaging Data"
task: ""
lineage_type: import
upstream_source: https://github.com/mims-harvard/ToolUniverse/blob/e2520a96/skills/tooluniverse-image-analysis/references/statistical_analysis.md
upstream_sha: e2520a96
imported_at: 2026-06-26
prompt_class: prompt
upstream_changes: accepted
author: upstream
validated: false
---

# Statistical Analysis for Imaging Data

Complete reference for statistical tests, effect sizes, power analysis, and regression modeling on microscopy measurement data.

---

## Table of Contents

1. [Descriptive Statistics](#descriptive-statistics)
2. [Normality Testing](#normality-testing)
3. [Two-Group Comparisons](#two-group-comparisons)
4. [Multiple Comparisons](#multiple-comparisons)
5. [Two-Way ANOVA](#two-way-anova)
6. [Effect Sizes](#effect-sizes)
7. [Power Analysis](#power-analysis)
8. [Regression Modeling](#regression-modeling)
9. [Model Comparison](#model-comparison)

---

## Descriptive Statistics

### Grouped Summary Statistics

```python
import pandas as pd
import numpy as np

def grouped_summary(df, group_cols, measure_col):
    """Calculate summary statistics by group.

    Returns DataFrame with Mean, SD, SEM, Median, Min, Max, N per group.
    """
    if isinstance(group_cols, str):
        group_cols = [group_cols]

    summary = df.groupby(group_cols)[measure_col].agg(
        Mean='mean',
        SD='std',
        Median='median',
        Min='min',
        Max='max',
        N='count'
    ).reset_index()

    summary['SEM'] = summary['SD'] / np.sqrt(summary['N'])

    return summary
```

### Percent Reduction

```python
def percent_reduction(df, group_col, measure_col, reference_group, comparison_group):
    """Calculate percent reduction of comparison vs reference.

    Returns: (percent_reduction, ref_mean, comp_mean)
    """
    ref_mean = df[df[group_col] == reference_group][measure_col].mean()
    comp_mean = df[df[group_col] == comparison_group][measure_col].mean()
    pct_reduction = ((ref_mean - comp_mean) / ref_mean) * 100
    return pct_reduction, ref_mean, comp_mean
```

### Relative Proportion

```python
def relative_proportion(df, group_col, measure_col, numerator_group, denominator_group):
    """Calculate relative proportion (as percentage) of one group vs another.

    Returns: proportion as percentage
    """
    num_mean = df[df[group_col] == numerator_group][measure_col].mean()
    den_mean = df[df[group_col] == denominator_group][measure_col].mean()
    return (num_mean / den_mean) * 100
```

---

## Normality Testing

### Shapiro-Wilk Test

```python
from scipy import stats

def shapiro_wilk_test(data):
    """Perform Shapiro-Wilk test for normality.

    Returns: (W_statistic, p_value)

    Interpretation:
    - p < 0.05: Data is NOT normally distributed
    - p >= 0.05: Data is normally distributed
    """
    stat, pvalue = stats.shapiro(data)
    return stat, pvalue

# Example usage
data = df[df['Condition'] == 'Control']['Measurement']
w_stat, p_val = shapiro_wilk_test(data)
print(f"Shapiro-Wilk W={w_stat:.4f}, p={p_val:.4f}")
```

---

## Two-Group Comparisons

### Independent T-Test

```python
def independent_ttest(group1, group2, equal_var=True):
    """Perform independent two-sample t-test.

    Args:
        group1, group2: array-like data
        equal_var: If True, use standard t-test. If False, use Welch's t-test

    Returns: (t_statistic, p_value)
    """
    t_stat, p_val = stats.ttest_ind(group1, group2, equal_var=equal_var)
    return t_stat, p_val
```

### Mann-Whitney U Test (non-parametric)

```python
def mann_whitney_test(group1, group2):
    """Perform Mann-Whitney U test (non-parametric alternative to t-test).

    Use when data is NOT normally distributed.

    Returns: (U_statistic, p_value)
    """
    u_stat, p_val = stats.mannwhitneyu(group1, group2, alternative='two-sided')
    return u_stat, p_val
```

---

## Multiple Comparisons

### Dunnett's Test

**CRITICAL**: Dunnett's test compares each treatment group to a single control group. R uses `multcomp::glht()` with `mcp(Strain_Ratio = "Dunnett")`. Python equivalent uses scipy.

```python
from scipy import stats

def dunnetts_test_scipy(df, group_col, value_col, control_group, alpha=0.05):
    """Dunnett's test using scipy.stats.dunnett (scipy >= 1.10).

    This is the preferred method - uses exact Dunnett distribution.

    Args:
        df: DataFrame
        group_col: Column with group labels
        value_col: Column with measurements
        control_group: Label of the control group
        alpha: Significance level (default 0.05)

    Returns: DataFrame with group, p_value, statistic, significant
    """
    groups = sorted(df[group_col].unique())
    control_data = df[df[group_col] == control_group][value_col].values
    treatment_groups = [g for g in groups if g != control_group]

    treatment_data = [df[df[group_col] == g][value_col].values for g in treatment_groups]

    # scipy.stats.dunnett: compare multiple treatment groups against control
    result = stats.dunnett(*treatment_data, control=control_data, alternative='two-sided')

    results = []
    for i, tg in enumerate(treatment_groups):
        results.append({
            'group': tg,
            'control': control_group,
            'p_value': result.pvalue[i],
            'statistic': result.statistic[i],
            'significant': result.pvalue[i] < alpha
        })

    return pd.DataFrame(results)
```

### Combined Dunnett's Test (Two Measures)

```python
def dunnett_area_circularity(df, group_col, area_col, circ_col, control_group, alpha=0.05):
    """Run Dunnett's test on both area and circularity.

    Returns: dict with:
        - area_results: Dunnett results for area
        - circ_results: Dunnett results for circularity
        - merged: Combined results
        - equivalent_in_both: groups NOT significant in EITHER (equivalent to control)
        - different_in_both: groups significant in BOTH (different from control)
    """
    area_dunnett = dunnetts_test_scipy(df, group_col, area_col, control_group, alpha)
    circ_dunnett = dunnetts_test_scipy(df, group_col, circ_col, control_group, alpha)

    # Merge results
    merged = area_dunnett[['group', 'p_value', 'significant']].merge(
        circ_dunnett[['group', 'p_value', 'significant']],
        on='group', suffixes=('_area', '_circ')
    )

    # Groups NOT significant in EITHER (equivalent to control in both)
    both_equiv = merged[~merged['significant_area'] & ~merged['significant_circ']]['group'].tolist()

    # Groups significant in both (different from control in both)
    both_diff = merged[merged['significant_area'] & merged['significant_circ']]['group'].tolist()

    return {
        'area_results': area_dunnett,
        'circ_results': circ_dunnett,
        'merged': merged,
        'equivalent_in_both': both_equiv,
        'different_in_both': both_diff,
    }
```

### Tukey HSD (all pairwise comparisons)

```python
from statsmodels.stats.multicomp import pairwise_tukeyhsd

def tukey_hsd_test(df, value_col, group_col, alpha=0.05):
    """Perform Tukey HSD test for all pairwise comparisons.

    Use when you want to compare ALL groups to EACH OTHER (not just vs control).

    Returns: Tukey HSD result object with summary table
    """
    result = pairwise_tukeyhsd(
        endog=df[value_col],
        groups=df[group_col],
        alpha=alpha
    )
    return result
```

---

## Two-Way ANOVA

```python
import statsmodels.api as sm
from statsmodels.formula.api import ols

def two_way_anova(df, dependent_var, factor1, factor2, anova_type=2):
    """Perform two-way ANOVA with interaction term.

    Args:
        df: DataFrame
        dependent_var: Column name of dependent variable
        factor1: First factor column name
        factor2: Second factor column name
        anova_type: Type of sum of squares (1, 2, or 3)
            - Type 1: Sequential (order matters)
            - Type 2: Hierarchical (recommended for balanced designs)
            - Type 3: Marginal (recommended for unbalanced designs)

    Returns: ANOVA table as DataFrame with columns:
        sum_sq, df, F, PR(>F) for each term including interaction
    """
    formula = f'{dependent_var} ~ C({factor1}) * C({factor2})'
    model = ols(formula, data=df).fit()
    anova_table = sm.stats.anova_lm(model, typ=anova_type)
    return anova_table


def extract_anova_interaction(anova_table, factor1, factor2):
    """Extract interaction F-statistic and p-value from ANOVA table.

    Returns: (F_statistic, p_value)
    """
    interaction_key = f'C({factor1}):C({factor2})'
    f_stat = anova_table.loc[interaction_key, 'F']
    p_val = anova_table.loc[interaction_key, 'PR(>F)']
    return f_stat, p_val

# Example usage
anova_result = two_way_anova(df, 'NeuN_count', 'Condition', 'Sex')
print(anova_result)

# Extract interaction
f_stat, p_val = extract_anova_interaction(anova_result, 'Condition', 'Sex')
print(f"Interaction F={f_stat:.3f}, p={p_val:.4f}")
```

---

## Effect Sizes

### Cohen's d

```python
def cohens_d(group1, group2):
    """Calculate Cohen's d using pooled standard deviation.

    This matches the standard formula: d = (mean1 - mean2) / sd_pooled
    where sd_pooled = sqrt(((n1-1)*s1^2 + (n2-1)*s2^2) / (n1+n2-2))

    NOTE: Uses pandas .std() which defaults to ddof=1 (sample std).

    Interpretation:
    - |d| < 0.2: Small effect
    - |d| = 0.2-0.5: Small to medium effect
    - |d| = 0.5-0.8: Medium to large effect
    - |d| > 0.8: Large effect

    Returns: Cohen's d value
    """
    n1, n2 = len(group1), len(group2)
    s1, s2 = group1.std(), group2.std()  # ddof=1 by default in pandas

    sd_pooled = np.sqrt(((n1 - 1) * s1**2 + (n2 - 1) * s2**2) / (n1 + n2 - 2))
    d = (group1.mean() - group2.mean()) / sd_pooled

    return d

# Example usage
control = df[df['Condition'] == 'Control']['Measurement']
treatment = df[df['Condition'] == 'Treatment']['Measurement']
effect_size = cohens_d(control, treatment)
print(f"Cohen's d = {effect_size:.3f}")
```

---

## Power Analysis

### Sample Size Calculation

```python
from statsmodels.stats.power import TTestIndPower

def power_analysis_sample_size(effect_size, alpha=0.05, power=0.8, alternative='two-sided'):
    """Calculate required sample size per group for a two-sample t-test.

    Args:
        effect_size: Cohen's d (can be negative, absolute value used)
        alpha: Significance level (Type I error rate)
        power: Desired statistical power (1 - Type II error rate)
        alternative: 'two-sided', 'larger', or 'smaller'

    Returns: Required sample size per group (rounded up to integer)

    Example:
        # To detect effect size of 0.8 with 80% power
        n = power_analysis_sample_size(0.8, alpha=0.05, power=0.8)
        print(f"Need {n} samples per group")
    """
    analysis = TTestIndPower()
    # Use absolute value of effect size for sample size calculation
    n = analysis.solve_power(
        effect_size=abs(effect_size),
        alpha=alpha,
        power=power,
        alternative=alternative
    )
    return int(np.ceil(n))
```

### Post-hoc Power Calculation

```python
def calculate_achieved_power(n_per_group, effect_size, alpha=0.05, alternative='two-sided'):
    """Calculate achieved power given sample size and effect size.

    Args:
        n_per_group: Sample size per group
        effect_size: Cohen's d
        alpha: Significance level
        alternative: 'two-sided', 'larger', or 'smaller'

    Returns: Achieved power (0-1)
    """
    analysis = TTestIndPower()
    power = analysis.solve_power(
        effect_size=abs(effect_size),
        nobs1=n_per_group,
        alpha=alpha,
        alternative=alternative
    )
    return power
```

---

## Regression Modeling

### Data Preparation for Co-culture Ratios

```python
def prepare_ratio_data(df, strain_col='StrainNumber', ratio_col='Ratio',
                       area_col='Area', exclude_strains=None):
    """Prepare co-culture ratio data for regression analysis.

    Converts ratio strings (e.g., "3:1") to frequency fractions.
    Filters out pure strains if requested.

    Args:
        df: DataFrame with swarming data
        strain_col: Column with strain identifiers
        ratio_col: Column with ratio strings
        area_col: Column with area measurements
        exclude_strains: List of strain IDs to exclude (e.g., pure strains)

    Returns: DataFrame with Frequency_rhlI column added
    """
    result = df.copy()

    if exclude_strains:
        result = result[~result[strain_col].isin(exclude_strains)]

    # Parse ratio into frequency
    # Ratio format: "rhlI_D:lasI_D" (e.g., "3:1" means 3 parts rhlI, 1 part lasI)
    ratio_parts = result[ratio_col].str.split(':', expand=True).astype(int)
    result['rhlI_D'] = ratio_parts[0]
    result['lasI_D'] = ratio_parts[1]
    result['Frequency_rhlI'] = result['rhlI_D'] / (result['rhlI_D'] + result['lasI_D'])

    return result
```

### Polynomial Regression

```python
def fit_polynomial_model(df, x_col, y_col, degree=2):
    """Fit polynomial regression model.

    Equivalent to R: lm(y ~ poly(x, degree, raw=TRUE))

    Args:
        df: DataFrame
        x_col: Predictor column name
        y_col: Response column name
        degree: Polynomial degree (2=quadratic, 3=cubic)

    Returns: dict with model, coefficients, R-squared, F-statistic, p-value,
             peak_frequency, peak_value, peak_ci
    """
    x = df[x_col].values
    y = df[y_col].values

    # Build design matrix for polynomial
    X_poly = np.column_stack([x**i for i in range(1, degree+1)])
    X = sm.add_constant(X_poly)

    model = sm.OLS(y, X).fit()

    result = {
        'model': model,
        'coefficients': model.params,
        'r_squared': model.rsquared,
        'adj_r_squared': model.rsquared_adj,
        'f_statistic': model.fvalue,
        'f_pvalue': model.f_pvalue,
        'aic': model.aic,
        'bic': model.bic,
        'summary': model.summary(),
    }

    # Find peak (maximum) using calculus
    if degree == 2:
        # y = b0 + b1*x + b2*x^2
        # Peak at x = -b1 / (2*b2)
        b1, b2 = model.params[1], model.params[2]
        peak_x = -b1 / (2 * b2)
    elif degree == 3:
        # y = b0 + b1*x + b2*x^2 + b3*x^3
        # Derivative: b1 + 2*b2*x + 3*b3*x^2 = 0
        b1, b2, b3 = model.params[1], model.params[2], model.params[3]
        discriminant = (2*b2)**2 - 4*(3*b3)*b1
        if discriminant >= 0:
            x1 = (-2*b2 + np.sqrt(discriminant)) / (2*3*b3)
            x2 = (-2*b2 - np.sqrt(discriminant)) / (2*3*b3)
            # Choose the one that's a maximum (second derivative < 0)
            candidates = [x1, x2]
            peak_x = None
            for cx in candidates:
                second_deriv = 2*b2 + 6*b3*cx
                if second_deriv < 0:
                    peak_x = cx
                    break
            if peak_x is None:
                peak_x = candidates[0]  # fallback
        else:
            peak_x = x.mean()  # no real critical points
    else:
        # For higher degrees, use numerical optimization
        from scipy.optimize import minimize_scalar
        poly_func = lambda xx: -sum(model.params[i] * xx**i for i in range(degree+1))
        opt = minimize_scalar(poly_func, bounds=(x.min(), x.max()), method='bounded')
        peak_x = opt.x

    # Predict at peak and get confidence interval
    X_peak = np.array([[1] + [peak_x**i for i in range(1, degree+1)]])
    peak_pred = model.get_prediction(X_peak)
    peak_value = peak_pred.predicted_mean[0]
    peak_ci = peak_pred.conf_int(alpha=0.05)[0]

    result['peak_x'] = peak_x
    result['peak_value'] = peak_value
    result['peak_ci_lower'] = peak_ci[0]
    result['peak_ci_upper'] = peak_ci[1]

    return result
```

### Natural Spline Regression

**CRITICAL**: This must match R's `lm(Area ~ ns(Frequency_rhlI, df=4))`.

```python
from patsy import dmatrix

def fit_natural_spline_model(df, x_col, y_col, spline_df=4):
    """Fit natural spline regression model.

    Equivalent to R: lm(y ~ ns(x, df=spline_df))
    Uses patsy's cr() with explicit quantile knots to match R's ns().

    CRITICAL: R's ns(x, df=N) places N-1 internal knots at equally-spaced
    quantiles (25th, 50th, 75th for df=4). patsy's cr(df=N) does NOT place
    knots at the same locations by default. You MUST specify knots explicitly.

    Args:
        df: DataFrame
        x_col: Predictor column name
        y_col: Response column name
        spline_df: Degrees of freedom for natural spline basis

    Returns: dict with model, R-squared, F-statistic, p-value,
             peak_frequency, peak_value, peak_ci
    """
    x = df[x_col].values
    y = df[y_col].values

    # Match R's ns() knot placement: df-1 internal knots at equally-spaced quantiles
    n_internal_knots = spline_df - 1
    quantile_pcts = np.linspace(100.0 / (n_internal_knots + 1),
                                100.0 * n_internal_knots / (n_internal_knots + 1),
                                n_internal_knots)
    knots = np.percentile(x, quantile_pcts)
    knot_str = ", ".join([str(k) for k in knots])

    # Create natural spline basis using patsy's cr() with explicit knots
    formula_str = f"cr({x_col}, knots=[{knot_str}]) - 1"
    X_spline = np.array(dmatrix(formula_str, df))
    X = sm.add_constant(X_spline)

    model = sm.OLS(y, X).fit()

    result = {
        'model': model,
        'r_squared': model.rsquared,
        'adj_r_squared': model.rsquared_adj,
        'f_statistic': model.fvalue,
        'f_pvalue': model.f_pvalue,
        'aic': model.aic,
        'bic': model.bic,
    }

    # Find peak by evaluating over fine grid
    x_grid = np.linspace(x.min(), x.max(), 1000)
    grid_df = pd.DataFrame({x_col: x_grid})
    X_grid_spline = np.array(dmatrix(formula_str, grid_df))
    X_grid = sm.add_constant(X_grid_spline)

    predictions = model.get_prediction(X_grid)
    pred_mean = predictions.predicted_mean
    pred_ci = predictions.conf_int(alpha=0.05)

    # Find peak
    max_idx = np.argmax(pred_mean)

    result['peak_x'] = x_grid[max_idx]
    result['peak_value'] = pred_mean[max_idx]
    result['peak_ci_lower'] = pred_ci[max_idx, 0]
    result['peak_ci_upper'] = pred_ci[max_idx, 1]
    result['predictions'] = pred_mean
    result['x_grid'] = x_grid
    result['ci_lower'] = pred_ci[:, 0]
    result['ci_upper'] = pred_ci[:, 1]

    return result
```

---

## Model Comparison

```python
def compare_regression_models(models_dict):
    """Compare multiple regression models.

    Args:
        models_dict: Dict mapping model_name -> result dict from fit_* functions

    Returns: DataFrame with model comparison metrics

    Example:
        models = {
            'quadratic': fit_polynomial_model(df, 'x', 'y', degree=2),
            'cubic': fit_polynomial_model(df, 'x', 'y', degree=3),
            'spline': fit_natural_spline_model(df, 'x', 'y', spline_df=4)
        }
        comparison = compare_regression_models(models)
        print(comparison)
    """
    comparison = []
    for name, result in models_dict.items():
        comparison.append({
            'model': name,
            'r_squared': result['r_squared'],
            'adj_r_squared': result['adj_r_squared'],
            'f_statistic': result['f_statistic'],
            'f_pvalue': result['f_pvalue'],
            'aic': result.get('aic'),
            'bic': result.get('bic'),
            'peak_x': result.get('peak_x'),
            'peak_value': result.get('peak_value'),
            'peak_ci_lower': result.get('peak_ci_lower'),
            'peak_ci_upper': result.get('peak_ci_upper'),
        })

    comparison_df = pd.DataFrame(comparison)
    comparison_df['best_r2'] = comparison_df['r_squared'] == comparison_df['r_squared'].max()
    comparison_df['best_aic'] = comparison_df['aic'] == comparison_df['aic'].min()
    comparison_df['best_bic'] = comparison_df['bic'] == comparison_df['bic'].min()

    return comparison_df.sort_values('r_squared', ascending=False)
```

---

## Answer Extraction Patterns

### Answer formatting

```python
def format_answer(value, question_type):
    """Format answer for reporting.

    Args:
        value: Numeric value to format
        question_type: Type of answer expected

    Returns: Formatted value
    """
    if question_type == "nearest_thousand":
        return int(round(value, -3))
    elif question_type == "percentage_int":
        return int(round(value))
    elif question_type == "percentage_2dec":
        return round(value, 2)
    elif question_type == "cohen_d":
        return round(value, 3)
    elif question_type == "statistic_3dec":
        return round(value, 3)
    elif question_type == "sample_size":
        return int(np.ceil(value))
    elif question_type == "count":
        return int(value)
    elif question_type == "r_squared":
        return round(value, 2)
    elif question_type == "p_value":
        if value < 0.0001:
            return f"{value:.2e}"
        else:
            return round(value, 4)
    else:
        return value
```

---

## Complete Example: Spline Regression Workflow

```python
import pandas as pd
import numpy as np

# Load data
df = pd.read_csv("Swarm_2.csv")

# Prepare ratio data
df_coculture = prepare_ratio_data(
    df,
    strain_col='StrainNumber',
    ratio_col='Ratio',
    area_col='Area',
    exclude_strains=['1', '98']  # Exclude pure strains
)

# Fit models
quadratic = fit_polynomial_model(df_coculture, 'Frequency_rhlI', 'Area', degree=2)
cubic = fit_polynomial_model(df_coculture, 'Frequency_rhlI', 'Area', degree=3)
spline = fit_natural_spline_model(df_coculture, 'Frequency_rhlI', 'Area', spline_df=4)

# Compare models
models = {'quadratic': quadratic, 'cubic': cubic, 'spline': spline}
comparison = compare_regression_models(models)
print(comparison)

# Best model by R-squared
best_model_name = comparison.iloc[0]['model']
best_model = models[best_model_name]

# Extract answers
print(f"Best model: {best_model_name}")
print(f"R-squared: {best_model['r_squared']:.4f}")
print(f"Peak frequency: {best_model['peak_x']:.4f}")
print(f"Peak area: {best_model['peak_value']:.1f}")
print(f"95% CI: [{best_model['peak_ci_lower']:.1f}, {best_model['peak_ci_upper']:.1f}]")
```