Below is a self-contained, reproducible Python script that uses matplotlib to visualize experimental and control observations. It includes data loading, preprocessing (NaN handling and robust winsorization), statistical estimation (bootstrap confidence intervals), and three plots: distribution, mean with error intervals, and scatter correlation with a bootstrap-based regression band.
Usage:
- If you have a CSV, provide either:
- Wide format: columns "control" and "experiment" (paired).
- Long format: columns "group" (values: "control" or "experiment") and "value". Optional "subject_id" for pairing.
- If no input is provided, the script generates synthetic reproducible data.
- Example:
- python visualize_groups.py --input your_data.csv --outdir figures
- python visualize_groups.py (uses synthetic data)
Required: Python 3.9+, numpy, pandas, matplotlib. Optional: scipy for KDE. If SciPy is not available, the script falls back to histograms without KDE.
Script (save as visualize_groups.py):
import argparse
import sys
from pathlib import Path
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
Optional SciPy
try:
from scipy.stats import gaussian_kde
SCIPY_AVAILABLE = True
except Exception:
SCIPY_AVAILABLE = False
def set_seed(seed: int = 42):
np.random.seed(seed)
def winsorize_by_mad(x: np.ndarray, zmax: float = 3.0) -> np.ndarray:
"""
Robust winsorization based on MAD. Caps values beyond zmax * MAD from the median.
"""
x = np.asarray(x).astype(float)
median = np.median(x)
mad = np.median(np.abs(x - median))
if mad == 0:
return x # No variability; nothing to winsorize
z = (x - median) / (1.4826 * mad) # 1.4826 scales MAD to be comparable to std under normality
x_w = x.copy()
x_w[z > zmax] = median + zmax * 1.4826 * mad
x_w[z < -zmax] = median - zmax * 1.4826 * mad
return x_w
def bootstrap_mean_ci(x: np.ndarray, n_boot: int = 5000, alpha: float = 0.05, seed: int = 123):
"""
Bootstrap mean and 95% CI (percentile) for a 1D array.
Returns: mean, ci_low, ci_high, se_boot
"""
x = np.asarray(x).astype(float)
x = x[~np.isnan(x)]
rng = np.random.default_rng(seed)
boots = rng.choice(x, (n_boot, x.size), replace=True).mean(axis=1)
mean = x.mean()
se_boot = boots.std(ddof=1)
ci_low = np.quantile(boots, alpha / 2)
ci_high = np.quantile(boots, 1 - alpha / 2)
return mean, ci_low, ci_high, se_boot
def fisher_r_ci(r: float, n: int, alpha: float = 0.05):
"""
Fisher z-transformation CI for Pearson correlation r with sample size n.
Returns (ci_low, ci_high). Assumes |r|<1 and n>3.
"""
r = np.clip(r, -0.999999, 0.999999)
z = np.arctanh(r)
se = 1 / np.sqrt(n - 3)
z_low = z - 1.96 * se
z_high = z + 1.96 * se
return np.tanh(z_low), np.tanh(z_high)
def load_data(csv_path: Path | None):
"""
Load data from CSV, or generate synthetic paired data if csv_path is None.
Returns: control, experiment, paired (bool)
"""
if csv_path is None:
# Synthetic, reproducible paired data
set_seed(42)
n = 200
control = np.random.normal(loc=50, scale=10, size=n)
effect = 5.0
experiment = control + effect + np.random.normal(loc=0, scale=8, size=n)
# Inject some missing values
miss_idx = np.random.choice(n, size=int(0.05 * n), replace=False)
control[miss_idx[: len(miss_idx)//2]] = np.nan
experiment[miss_idx[len(miss_idx)//2:]] = np.nan
return control, experiment, True
df = pd.read_csv(csv_path)
cols = {c.lower(): c for c in df.columns}
if "control" in cols and "experiment" in cols:
control = df[cols["control"]].to_numpy(dtype=float)
experiment = df[cols["experiment"]].to_numpy(dtype=float)
return control, experiment, True
if "group" in cols and "value" in cols:
# Long format; attempt pairing by subject_id if present
gcol = cols["group"]
vcol = cols["value"]
df = df[[gcol, vcol] + ([cols["subject_id"]] if "subject_id" in cols else [])].copy()
df[gcol] = df[gcol].str.lower().str.strip()
df = df[df[gcol].isin(["control", "experiment"])]
if "subject_id" in cols:
# Pivot to paired wide
wide = df.pivot_table(index=cols["subject_id"], columns=gcol, values=vcol, aggfunc="mean")
control = wide.get("control", pd.Series(index=wide.index, dtype=float)).to_numpy()
experiment = wide.get("experiment", pd.Series(index=wide.index, dtype=float)).to_numpy()
return control, experiment, True
else:
# Unpaired
control = df.loc[df[gcol] == "control", vcol].to_numpy(dtype=float)
experiment = df.loc[df[gcol] == "experiment", vcol].to_numpy(dtype=float)
return control, experiment, False
raise ValueError("CSV must contain either (control, experiment) columns or (group, value) [and optionally subject_id].")
def preprocess(control: np.ndarray, experiment: np.ndarray, paired: bool, zmax: float = 3.0):
"""
- Remove NaN/inf.
- Robust winsorization by MAD.
- If paired, keep rows where both are valid.
"""
control = np.asarray(control, dtype=float)
experiment = np.asarray(experiment, dtype=float)
if paired:
mask = np.isfinite(control) & np.isfinite(experiment)
control, experiment = control[mask], experiment[mask]
control = winsorize_by_mad(control, zmax=zmax)
experiment = winsorize_by_mad(experiment, zmax=zmax)
else:
control = control[np.isfinite(control)]
experiment = experiment[np.isfinite(experiment)]
control = winsorize_by_mad(control, zmax=zmax)
experiment = winsorize_by_mad(experiment, zmax=zmax)
return control, experiment
def kde_or_none(x: np.ndarray, grid: np.ndarray):
if SCIPY_AVAILABLE and x.size > 1:
try:
kde = gaussian_kde(x)
return kde(grid)
except Exception:
return None
return None
def plot_distributions(control: np.ndarray, experiment: np.ndarray, outdir: Path):
plt.figure(figsize=(8, 5), dpi=120)
bins = max(10, int(np.sqrt(control.size + experiment.size)))
# Histograms
plt.hist(control, bins=bins, density=True, alpha=0.4, color="#1f77b4", label=f"Control (n={control.size})")
plt.hist(experiment, bins=bins, density=True, alpha=0.4, color="#ff7f0e", label=f"Experiment (n={experiment.size})")
# KDE curves if available
xmin = min(np.min(control), np.min(experiment))
xmax = max(np.max(control), np.max(experiment))
grid = np.linspace(xmin, xmax, 200)
d_control = kde_or_none(control, grid)
d_experiment = kde_or_none(experiment, grid)
if d_control is not None:
plt.plot(grid, d_control, color="#1f77b4", lw=2)
if d_experiment is not None:
plt.plot(grid, d_experiment, color="#ff7f0e", lw=2)
# Means
plt.axvline(control.mean(), color="#1f77b4", ls="--", lw=1)
plt.axvline(experiment.mean(), color="#ff7f0e", ls="--", lw=1)
plt.title("Distribution: Control vs Experiment")
plt.xlabel("Value")
plt.ylabel("Density")
plt.legend()
plt.tight_layout()
fp = outdir / "01_distributions.png"
plt.savefig(fp, bbox_inches="tight")
plt.close()
return fp
def plot_error_intervals(control: np.ndarray, experiment: np.ndarray, paired: bool, outdir: Path):
m_c, ci_c_low, ci_c_high, se_c = bootstrap_mean_ci(control)
m_e, ci_e_low, ci_e_high, se_e = bootstrap_mean_ci(experiment)
# Effect size: mean difference
if paired and control.size == experiment.size:
diff = experiment - control
md, md_low, md_high, _ = bootstrap_mean_ci(diff)
effect_text = f"Mean diff (paired): {md:.2f} [{md_low:.2f}, {md_high:.2f}]"
else:
md = m_e - m_c
# Unpaired bootstrap of difference
rng = np.random.default_rng(123)
n_boot = 5000
boots = []
for _ in range(n_boot):
b_c = rng.choice(control, control.size, replace=True).mean()
b_e = rng.choice(experiment, experiment.size, replace=True).mean()
boots.append(b_e - b_c)
md_low, md_high = np.quantile(boots, [0.025, 0.975])
effect_text = f"Mean diff (unpaired): {md:.2f} [{md_low:.2f}, {md_high:.2f}]"
plt.figure(figsize=(7, 5), dpi=120)
means = [m_c, m_e]
ci_lows = [ci_c_low, ci_e_low]
ci_highs = [ci_c_high, ci_e_high]
x = np.arange(2)
colors = ["#1f77b4", "#ff7f0e"]
plt.errorbar(x, means, yerr=[np.array(means) - np.array(ci_lows), np.array(ci_highs) - np.array(means)],
fmt="o", capsize=5, color="black", ecolor="black")
plt.scatter(x, means, c=colors, s=80, zorder=3)
plt.xticks(x, ["Control", "Experiment"])
plt.ylabel("Mean ± 95% CI")
plt.title("Group Means with 95% Confidence Intervals")
plt.text(0.5, 0.05, effect_text, transform=plt.gca().transAxes, ha="center")
plt.tight_layout()
fp = outdir / "02_error_intervals.png"
plt.savefig(fp, bbox_inches="tight")
plt.close()
return fp
def plot_scatter_correlation(control: np.ndarray, experiment: np.ndarray, paired: bool, outdir: Path):
if not paired or control.size != experiment.size:
# Skip if not paired
return None
x = control
y = experiment
n = x.size
# Pearson r
x_c = x - x.mean()
y_c = y - y.mean()
r = (x_c @ y_c) / (np.sqrt((x_c**2).sum()) * np.sqrt((y_c**2).sum()))
r_low, r_high = fisher_r_ci(r, n)
# Regression line and bootstrap band
beta = np.polyfit(x, y, deg=1) # y = beta[0]*x + beta[1]
x_line = np.linspace(x.min(), x.max(), 200)
y_line = beta[0] * x_line + beta[1]
# Bootstrap regression uncertainty
rng = np.random.default_rng(123)
n_boot = 2000
y_boot = np.empty((n_boot, x_line.size))
idx = np.arange(n)
for i in range(n_boot):
b = rng.choice(idx, size=n, replace=True)
xb, yb = x[b], y[b]
bcoef = np.polyfit(xb, yb, deg=1)
y_boot[i] = bcoef[0] * x_line + bcoef[1]
band_low = np.quantile(y_boot, 0.025, axis=0)
band_high = np.quantile(y_boot, 0.975, axis=0)
plt.figure(figsize=(7, 6), dpi=120)
plt.scatter(x, y, alpha=0.7, s=40, color="#2ca02c", edgecolor="white", linewidth=0.5)
plt.plot(x_line, y_line, color="black", lw=2, label="OLS fit")
plt.fill_between(x_line, band_low, band_high, color="gray", alpha=0.2, label="95% bootstrap band")
plt.xlabel("Control")
plt.ylabel("Experiment")
plt.title("Scatter and Correlation (Paired)")
plt.legend()
plt.text(0.05, 0.02, f"Pearson r = {r:.3f} [{r_low:.3f}, {r_high:.3f}] (95% CI)", transform=plt.gca().transAxes)
plt.tight_layout()
fp = outdir / "03_scatter_correlation.png"
plt.savefig(fp, bbox_inches="tight")
plt.close()
return fp
def main():
parser = argparse.ArgumentParser(description="Visualize experimental vs control observations.")
parser.add_argument("--input", type=str, default=None, help="CSV file path. See header format in script description.")
parser.add_argument("--outdir", type=str, default="figures", help="Output directory for figures.")
parser.add_argument("--zmax", type=float, default=3.0, help="MAD-based winsorization threshold.")
args = parser.parse_args()
outdir = Path(args.outdir)
outdir.mkdir(parents=True, exist_ok=True)
csv_path = Path(args.input) if args.input else None
try:
control, experiment, paired = load_data(csv_path)
except Exception as e:
print(f"Error loading data: {e}", file=sys.stderr)
sys.exit(1)
control, experiment = preprocess(control, experiment, paired=paired, zmax=args.zmax)
dist_fp = plot_distributions(control, experiment, outdir)
err_fp = plot_error_intervals(control, experiment, paired, outdir)
scat_fp = plot_scatter_correlation(control, experiment, paired, outdir)
print("Saved figures:")
print(f"- {dist_fp}")
print(f"- {err_fp}")
if scat_fp is not None:
print(f"- {scat_fp}")
else:
print("- Scatter correlation skipped (requires paired data).")
if name == "main":
main()
Notes on methodology:
- Preprocessing: NaN/inf values are removed. Robust winsorization uses MAD to cap extreme outliers without discarding data, which stabilizes estimates and plots.
- Estimation: Means and 95% confidence intervals use bootstrap percentiles, which are reliable for non-normal distributions. For correlations, Fisher’s z-transform provides an analytical 95% CI for r.
- Visualization:
- Distribution plots show histograms (and KDE when SciPy is available) with mean indicators.
- Error interval plot displays mean ± 95% CI per group and annotates the mean difference with its CI (paired or unpaired, based on data).
- Scatter correlation plot requires paired data and includes an OLS regression line with a bootstrap-based 95% uncertainty band.
