06 · 生物统计与实验数据分析

import numpy as np
from scipy import stats

data = np.array([3.2, 5.8, 1.5, 7.1, 2.4, 4.6, 3.9, 6.2])

data.mean()          # 均值
np.median(data)      # 中位数
data.std(ddof=1)     # 样本标准差（ddof=1 是科研标准！）
data.var(ddof=1)     # 样本方差

def mean_ci(data, confidence=0.95):
    n = len(data)
    mean = np.mean(data)
    se = stats.sem(data)                     # 标准误 = std / sqrt(n)
    t_crit = stats.t.ppf((1 + confidence) / 2, df=n - 1)
    margin = t_crit * se
    return mean - margin, mean + margin

print(mean_ci(data))   # (2.83, 5.57)

import matplotlib.pyplot as plt

# 生成正态分布数据并画直方图
data = np.random.normal(loc=5, scale=1, size=1000)
plt.hist(data, bins=30, density=True, alpha=0.7, color="teal")
# 叠加理论正态曲线
x = np.linspace(2, 8, 200)
plt.plot(x, stats.norm.pdf(x, 5, 1), color="orange", linewidth=2)
plt.show()

# 右偏数据示例（对数正态分布）
data_skewed = np.random.lognormal(mean=0, sigma=0.5, size=1000)
print(f"均值: {data_skewed.mean():.2f}, 中位数: {np.median(data_skewed):.2f}")
# 均值 > 中位数 → 右偏

data = np.random.normal(5, 1, 30)

stat, p = stats.shapiro(data)
print(f"Shapiro-Wilk: W = {stat:.4f}, p = {p:.4f}")

if p > 0.05:
    print("数据符合正态分布（P > 0.05）")
else:
    print("数据不符合正态分布（P < 0.05）")

stat, p = stats.kstest(data, "norm", args=(data.mean(), data.std(ddof=1)))
print(f"K-S: D = {stat:.4f}, p = {p:.4f}")

import scipy.stats as stats
import matplotlib.pyplot as plt

stats.probplot(data, dist="norm", plot=plt)
plt.title("Q-Q 图")
plt.show()
# 点越贴近对角线，越接近正态分布

group1 = np.random.normal(5, 1, 20)
group2 = np.random.normal(6, 1.5, 20)

stat, p = stats.levene(group1, group2)
print(f"Levene: W = {stat:.4f}, p = {p:.4f}")
# P > 0.05 → 方差齐

stat, p = stats.bartlett(group1, group2)

control = np.random.normal(5, 1, 12)
model   = np.random.normal(7, 1, 12)

# 默认 equal_var=True（方差齐）；方差不齐时 equal_var=False
t_stat, p_val = stats.ttest_ind(control, model, equal_var=True)
print(f"t = {t_stat:.4f}, p = {p_val:.4f}")

# 同时报告均值±SD
print(f"对照组: {control.mean():.2f} ± {control.std(ddof=1):.2f}")
print(f"模型组: {model.mean():.2f} ± {model.std(ddof=1):.2f}")

before = np.random.normal(6, 1, 12)
after  = before + np.random.normal(-1.5, 0.5, 12)  # 处理后下降

t_stat, p_val = stats.ttest_rel(before, after)
print(f"配对 t = {t_stat:.4f}, p = {p_val:.4f}")

group_a = np.random.normal(5, 1, 10)
group_b = np.random.normal(6, 1, 10)
group_c = np.random.normal(7, 1, 10)

f_stat, p_val = stats.f_oneway(group_a, group_b, group_c)
print(f"F = {f_stat:.4f}, p = {p_val:.4f}")

import statsmodels.api as sm
from statsmodels.formula.api import ols
import pandas as pd

# 构造数据：药物（drug）× 剂量（dose）
df = pd.DataFrame({
    "drug": np.repeat(["A", "B"], 20),
    "dose": np.tile(np.repeat(["low", "high"], 10), 2),
    "value": np.random.normal(5, 1, 40)
})

model = ols("value ~ C(drug) * C(dose)", data=df).fit()
anova_table = sm.stats.anova_lm(model, typ=2)
print(anova_table)
# 看 drug、dose、drug:dose（交互项）各自的 F 和 P 值

from statsmodels.stats.anova import AnovaRM

# df 需为长格式：subject, time, value
df_long = pd.DataFrame({
    "subject": np.repeat(range(10), 4),
    "time": np.tile(["t0", "t24", "t48", "t72"], 10),
    "value": np.random.normal(5, 1, 40)
})

result = AnovaRM(df_long, "value", "subject", within=["time"]).fit()
print(result)

group1 = np.random.normal(5, 1, 12)
group2 = np.random.normal(7, 1.5, 12)

u_stat, p_val = stats.mannwhitneyu(group1, group2, alternative="two-sided")
print(f"U = {u_stat:.4f}, p = {p_val:.4f}")

before = np.random.normal(6, 1, 12)
after  = before + np.random.normal(-1.5, 0.5, 12)

w_stat, p_val = stats.wilcoxon(before, after)
print(f"W = {w_stat:.4f}, p = {p_val:.4f}")

group_a = np.random.normal(5, 1, 10)
group_b = np.random.normal(6, 1.5, 10)
group_c = np.random.normal(7, 2, 10)

h_stat, p_val = stats.kruskal(group_a, group_b, group_c)
print(f"H = {h_stat:.4f}, p = {p_val:.4f}")

# 同一批对象多个时间点
data_matrix = np.random.randn(10, 4)  # 10个对象 × 4个时间点
stat, p_val = stats.friedmanchisquare(*[data_matrix[:, i] for i in range(4)])
print(f"χ² = {stat:.4f}, p = {p_val:.4f}")

# 2×2 列联表
table = np.array([[30, 10],   # 给药组：存活30，死亡10
                  [15, 25]])  # 对照组：存活15，死亡25

chi2, p_val, dof, expected = stats.chi2_contingency(table)
print(f"χ² = {chi2:.4f}, p = {p_val:.4f}")

from statsmodels.stats.multicomp import pairwise_tukeyhsd

# 长格式数据
df = pd.DataFrame({
    "group": np.repeat(["control", "low", "high"], 10),
    "value": np.concatenate([np.random.normal(5, 1, 10),
                             np.random.normal(6, 1, 10),
                             np.random.normal(7, 1, 10)])
})

result = pairwise_tukeyhsd(df["value"], df["group"], alpha=0.05)
print(result)
# 看 reject 列：True 表示该对差异显著

from statsmodels.stats.multitest import multipletests

p_values = [0.01, 0.04, 0.06, 0.2, 0.5]

# Bonferroni
bonf = multipletests(p_values, method="bonferroni")
print(bonf[1])   # 校正后的 P 值

# FDR (Benjamini-Hochberg)
fdr = multipletests(p_values, method="fdr_bh")
print(fdr[1])

x = np.random.normal(0, 1, 50)
y = 0.8 * x + np.random.normal(0, 0.5, 50)

r, p_val = stats.pearsonr(x, y)
print(f"r = {r:.4f}, p = {p_val:.4f}")

r, p_val = stats.spearmanr(x, y)
print(f"ρ = {r:.4f}, p = {p_val:.4f}")

import statsmodels.api as sm

x = np.random.normal(0, 1, 100)
y = 2 + 1.5 * x + np.random.normal(0, 0.5, 100)

X = sm.add_constant(x)          # 加截距
model = sm.OLS(y, X).fit()
print(model.summary())
# 关键输出：coef（系数）、P>|t|（P值）、R-squared（R²）

X_multi = np.random.randn(100, 3)
y = 1 + 0.5*X_multi[:,0] - 0.3*X_multi[:,1] + 0.2*X_multi[:,2] + np.random.randn(100)*0.2

X_multi = sm.add_constant(X_multi)
model = sm.OLS(y, X_multi).fit()
print(model.summary())

import statsmodels.api as sm

X = np.random.randn(200, 3)
y = (X[:, 0] + X[:, 1] > 0).astype(int)

X = sm.add_constant(X)
model = sm.Logit(y, X).fit()
print(model.summary())
# 关注 OR 值（优势比）和 P 值

from lifelines import KaplanMeierFitter

# 模拟数据：时间 + 是否发生事件（1=事件，0=删失）+ 分组
durations = np.random.exponential(10, 40)
events = np.random.binomial(1, 0.7, 40)
groups = np.repeat(["control", "drug"], 20)

kmf = KaplanMeierFitter()
for g in ["control", "drug"]:
    mask = groups == g
    kmf.fit(durations[mask], event_observed=events[mask], label=g)
    kmf.plot_survival_function()

plt.xlabel("Time (days)")
plt.ylabel("Survival probability")
plt.show()

from lifelines.statistics import logrank_test

mask_c = groups == "control"
mask_d = groups == "drug"
result = logrank_test(durations[mask_c], durations[mask_d],
                      event_observed_A=events[mask_c],
                      event_observed_B=events[mask_d])
print(f"log-rank p = {result.p_value:.4f}")

def cohens_d(group1, group2):
    n1, n2 = len(group1), len(group2)
    s1, s2 = group1.std(ddof=1), group2.std(ddof=1)
    sp = np.sqrt(((n1-1)*s1**2 + (n2-1)*s2**2) / (n1+n2-2))  # 合并标准差
    return (group1.mean() - group2.mean()) / sp

d = cohens_d(np.random.normal(5,1,20), np.random.normal(6,1,20))
print(f"Cohen's d = {d:.3f}")

from statsmodels.stats.power import TTestIndPower

# 独立样本 t 检验样本量估算
analysis = TTestIndPower()
n = analysis.solve_power(effect_size=0.8,   # Cohen's d
                         alpha=0.05,
                         power=0.80,
                         alternative="two-sided")
print(f"每组需要样本量: {np.ceil(n)}")

拿到数据
│
├─ 因变量是连续变量？
│   ├─ 否 → 分类变量 → 卡方检验 / Fisher 精确检验
│   └─ 是 ↓
│
├─ 比较几组？
│   ├─ 两组 → 数据正态？
│   │   ├─ 是 → 方差齐？
│   │   │   ├─ 是 → 独立样本 t 检验
│   │   │   └─ 否 → Welch's t 检验
│   │   └─ 否 → Mann-Whitney U 检验
│   │
│   └─ ≥3 组 → 数据正态？
│       ├─ 是 → 方差齐？
│       │   ├─ 是 → One-way ANOVA → Tukey HSD 事后比较
│       │   └─ 否 → Welch's ANOVA → Games-Howell
│       └─ 否 → Kruskal-Wallis → Dunn 事后比较
│
├─ 配对/重复测量？
│   ├─ 两组配对 → 正态？→ 配对 t 检验 / Wilcoxon 符号秩检验
│   └─ 多时间点 → 重复测量 ANOVA / Friedman 检验
│
├─ 两个连续变量关系 → 相关分析（Pearson / Spearman）
├─ 预测/影响因素 → 回归分析（线性 / Logistic）
└─ 生存时间 → Kaplan-Meier + log-rank

import matplotlib.pyplot as plt
import seaborn as sns

# 两组比较图（箱线图 + 散点 + 显著性标注）
plt.figure(figsize=(5, 5))
sns.boxplot(data=df, x="group", y="value", width=0.4, palette="Set2")
sns.stripplot(data=df, x="group", y="value", color="black", alpha=0.5, size=6)

# 添加显著性标注
x1, x2 = 0, 1
y_max = df["value"].max() + 1
plt.plot([x1, x1, x2, x2], [y_max, y_max+0.3, y_max+0.3, y_max], color="black", lw=1.5)
plt.text((x1+x2)/2, y_max+0.4, "**", ha="center", fontsize=14)

plt.ylim(df["value"].min()-1, y_max+1.5)
plt.show()