Resampling Statistics

Beeswarm, histogram, KDE, summary statistics, and big-box NHST resampling for Control (A) vs Treatment (B) T-cell count data. This is the static view of the notebook — click below to run it live in your browser.

▶ Run interactively in browser ↗

Setup — imports & data

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec
import seaborn as sns
from scipy.stats import gaussian_kde
from plot_utils import despine, finish_plot, vline_labeled, MAD


COLOR_A = '#2c99e3'  # blue  — Control (A)
COLOR_B = '#fc6d76'  # pink  — Treatment (B)

grp_A = np.array([53.2, 54.6, 55.1, 55.4, 55.8, 56.2, 55.7, 56.5, 56.6, 56.9, 57.2, 56.7, 57.9, 57.6, 57.8, 58.5, 58.9, 59. , 59.1, 58.7,
                  58.6, 58.8, 59.7, 59.4, 59.8, 60.1, 60.8, 60.5, 61. , 61.7, 61.3, 61.9, 61.5, 62.8, 62.5, 62.6, 62.9, 63.5, 63.8, 65.5,
                  65.1, 64.9, 65.2, 66. , 67.9, 67.8, 68.1, 69.5, 70. ])

grp_B = np.array([57.9, 59.6, 61. , 61.6, 62.4, 62.8, 63.1, 64. , 64.7, 65.2, 65.9, 66.2, 66.4, 66.5, 68.4, 68.9, 68.3, 67.8, 67.7, 68. ,
                  69. , 68.8, 69.5, 70. , 70.1, 70.8, 71.2, 71.3, 72. , 72.1, 71.1, 70.5, 71.1, 71.9, 71.5, 72.5, 73.1, 73.3, 72.8, 73.8,
                  73.5, 74.9, 75.2, 76.2, 75.8])

Figure 6.8 — Beeswarm plots

fig, ax = plt.subplots(figsize=(4, 4))
sns.swarmplot(data=[grp_A, grp_B], ax=ax, palette=[COLOR_A, COLOR_B], size=6)
ax.set_xticks([0, 1])
ax.set_xticklabels(['Control (A)', 'Treatment (B)'])
ax.set_ylim(50, 77)
finish_plot(ax, '', 'T-cell count\n(per unit area)')

Figure 6.9 — Overlapping histograms

bins=np.arange(50, 80, 2)
fig, ax = plt.subplots(figsize=(8, 3.5))
ax.hist(grp_A, bins=bins, color=COLOR_A, alpha=0.5, edgecolor='white', linewidth=1)
ax.hist(grp_B, bins=bins, color=COLOR_B, alpha=0.5, edgecolor='white', linewidth=1)
ax.text(0.97, 0.94, 'Control (A)',   color=COLOR_A, ha='right', va='top',
        transform=ax.transAxes, fontweight='bold')
ax.text(0.97, 0.78, 'Treatment (B)', color=COLOR_B, ha='right', va='top',
        transform=ax.transAxes, fontweight='bold')
ax.set_xlim(45, 85)
finish_plot(ax, 'T-cell count (per unit area)', 'Count')

Figure 6.10 — Histogram + KDE

bins=np.arange(50, 80, 2)
x = np.linspace(44, 88, 300)
fig, ax = plt.subplots(figsize=(8, 3.5))
ax.hist(grp_A, bins=bins, density=True, color=COLOR_A, alpha=0.3, edgecolor='white')
ax.hist(grp_B, bins=bins, density=True, color=COLOR_B, alpha=0.3, edgecolor='white')
ax.plot(x, gaussian_kde(grp_A)(x), color=COLOR_A, lw=2.5)
ax.plot(x, gaussian_kde(grp_B)(x), color=COLOR_B, lw=2.5)
ax.yaxis.set_major_formatter(plt.FuncFormatter(lambda y, _: f'{y*100:.0f}%'))
ax.text(0.97, 0.94, 'Control (A)',   color=COLOR_A, ha='right', va='top',
        transform=ax.transAxes, fontweight='bold')
ax.text(0.97, 0.78, 'Treatment (B)', color=COLOR_B, ha='right', va='top',
        transform=ax.transAxes, fontweight='bold')
ax.set_xlim(45, 85)
finish_plot(ax, 'T-cell count (per unit area)', 'Frequency')

Figure 6.11 — Beeswarm + boxplot + summary statistics

fig = plt.figure(figsize=(5.5, 7))
gs = gridspec.GridSpec(2, 1, height_ratios=[4.5, 2], hspace=0.05)
ax = fig.add_subplot(gs[0])

# Beeswarm
sns.swarmplot(data=[grp_A, grp_B], ax=ax, palette=[COLOR_A, COLOR_B], size=6)

# Boxplots — hide built-in median line, will draw manually below
bp = ax.boxplot([grp_A, grp_B], positions=[0.3, 1.3], widths=0.25,
                patch_artist=True,
                medianprops=dict(linewidth=0),
                whiskerprops=dict(linewidth=1.5),
                capprops=dict(linewidth=1.5),
                flierprops=dict(visible=False))
for patch, color in zip(bp['boxes'], [COLOR_A, COLOR_B]):
    patch.set_facecolor(color)
    patch.set_alpha(0.4)
    patch.set_zorder(2)

# Mean: dashed line inside box; Median: thick solid line, wider than box
for x, grp in zip([0.3, 1.3], [grp_A, grp_B]):
    ax.hlines(np.mean(grp),   x - 0.125, x + 0.125, colors='black', linewidth=2, zorder=4)
    ax.hlines(np.median(grp), x - 0.16, x + 0.16, colors='black', linewidth=4, zorder=5)

ax.set_xticks([])
ax.set_xlim(-0.4, 1.6)
ax.set_ylim(50, 77)
ax.set_ylabel('T-cell count\n(per unit area)')
ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)
ax.spines['bottom'].set_visible(False)
ax.tick_params(bottom=False)

# Summary stats — two sections aligned under each group
ax2 = fig.add_subplot(gs[1])
ax2.axis('off')

rows_lbl = ['Mean', 'MAD', 'Median']
vals_A = [f'{np.mean(grp_A):.0f}', f'{MAD(grp_A):.1f}', f'{np.median(grp_A):.1f}']
vals_B = [f'{np.mean(grp_B):.0f}', f'{MAD(grp_B):.1f}', f'{np.median(grp_B):.1f}']

for vals, bbox, header in [
    (vals_A, [0.04, 0.0, 0.44, 1.0], 'Control (A)'),
    (vals_B, [0.54, 0.0, 0.44, 1.0], 'Treatment (B)')
]:
    tbl = ax2.table(
        cellText=[[r, v] for r, v in zip(rows_lbl, vals)],
        colLabels=['', header],
        cellLoc='center', loc='center',
        bbox=bbox)
    tbl.auto_set_font_size(False)
    tbl.set_fontsize(11)
    for (r, c), cell in tbl.get_celld().items():
        cell.set_linewidth(0)
        if r == 0:
            cell.set_facecolor('white')
            cell.set_text_props(fontweight='bold')
        else:
            cell.set_facecolor('#eeeeee')
            cell.get_text().set_ha('left' if c == 0 else 'right')

plt.show()

Figure 6.12 — NHST: big-box resampling

dobs = np.median(grp_B) - np.median(grp_A)
bigbox = np.concatenate([grp_A, grp_B])
np.random.seed(10)
N = 10_000
ds = (np.median(np.random.choice(bigbox, size=(N, len(grp_B))), axis=1) -
      np.median(np.random.choice(bigbox, size=(N, len(grp_A))), axis=1))

n_extreme = int(np.sum(ds >= dobs) + np.sum(ds <= -dobs))
p_val = n_extreme / N
print(f"\u0394_obs = {dobs:.2f}  |  p = {p_val:.4f}  ({n_extreme} of {N} simulations)")

fig, ax = plt.subplots(figsize=(6, 2))
counts, _, _ = ax.hist(ds, bins=10, facecolor='none', edgecolor='black', linewidth=0.8)
vline_labeled(ax,  dobs, '\u0394_obs')
vline_labeled(ax, -dobs, '\u2212\u0394_obs', ls='--')
ax.text(dobs + 0.2, counts.max() * 0.7, 'p < 0.0001', va='center')
ax.set_xlim(-dobs * 1.6, dobs * 1.6)
finish_plot(ax, '\u0394\u1d62 (T-cells/unit area)', 'Count')