医学统计与生物信息学习平台

🌿 统计方法选择决策树

不知道用什么统计方法？回答几个问题帮你选

📊 统计方法

📊 统计方法 📈 ROC曲线 ⏳ 生存分析 📈 时间序列 🧮 样本量 ⚖️ 因果推断 🤖 高级模型 📈 建模 🏥 临床 🧬 资源 📝 案例

📐 统计公式速查

▼

方法	统计量	公式	条件
单样本t	t	t=(x̄-μ)/(s/√n)	正态分布
独立t	t	t=(x̄₁-x̄₂)/√(s²(1/n₁+1/n₂))	正态+方差齐
配对t	t	t=d̄/(sd/√n)	差值正态
卡方	χ²	χ²=Σ(O-E)²/E	期望≥5
Pearson r	r	r=Cov(X,Y)/(sx·sy)	双变量正态
线性回归F	F	F=MSR/MSE	线性+正态
Logistic OR	OR	OR=exp(β)	二分类Y
灵敏度Se	Se	Se=TP/(TP+FN)	金标准
特异度Sp	Sp	Sp=TN/(TN+FP)	金标准

📚 美国统计学会(ASA) 官方资源

📊 医学统计学方法

P 值解读

定义：P值：原假设H₀为真时，观察到当前或更极端结果的概率。需结合效应量和置信区间解读。

P值是统计学中最常用的概念。它衡量在原假设为真的条件下，观察到当前结果或更极端结果的概率。P<0.05常被视为统计学显落，但这只是约定俗成的阈值。ASA发布的官方声明强调，P值不能代表效应大小或结果重要性。

分析 -> 描述统计 -> 交叉表（分类变量）或分析 -> 比较均值 -> T检验（连续变量）
在输出结果中查看Sig.列即为P值
注意: SPSS默认输出双侧P值

# P值计算示例 set.seed(123) g1 <- rnorm(30, mean=100, sd=15) g2 <- rnorm(30, mean=110, sd=15) result <- t.test(g1, g2) print(result) cat("P:", result$p.value, "\n") # P值取决于: 效应量 + 样本量 + 变异度

Python

from scipy import stats g1 = np.random.normal(100,15,30) g2 = np.random.normal(110,15,30) t, p = stats.ttest_ind(g1, g2) print(f't={t:.3f}, p={p:.4f}') # 多重比较校正 from scipy.stats import multipletests _, corr, _, _ = multipletests([0.012,0.043,0.210,0.003], method='bonferroni') print(f'Corr: {corr}')

📋 P值案例

情景：比较两种降压药的疗效
结果：t=2.45, df=58, P=0.017
解读：P<0.05显示组间差异显著。但需关注效应量(Cohen d=0.63)和95%CI[1.2,11.8]

置信区间 (CI)

定义：参数估计的区间范围，同时展示效应大小与精度。95%CI不含无效值=统计显著。

CI是参数估计的区间范围。95%CI含义: 重复抽样100次，约95次区间会包含总体真实值。CI展示效应大小和精度，比P值提供更多临床相关信息。

分析 -> 比较均值 -> 独立样本T检验
在选项中设置置信区间百分比(默认95%)
查看差值的95%置信区间

# CI set.seed(123) data <- rnorm(50, 100, 15) t.test(data, conf.level=0.95) n <- length(data); se <- sd(data)/sqrt(n) ci <- mean(data)+c(-1,1)*qt(0.975,n-1)*se cat("95%CI:", round(ci,2)) g1 <- rnorm(30,130,10); g2 <- rnorm(30,140,12) t.test(g1,g2,conf.level=0.95)

Python

from scipy import stats data = np.random.normal(100,15,50) ci = stats.t.interval(0.95, len(data)-1, loc=np.mean(data), scale=stats.sem(data)) print(f'95%CI: {ci}')

📋 CI案例

结果：均值差=-8.5, 95%CI[-12.3,-4.7]
解读：不包含0→显著。最保守估计(-4.7)也有临床意义。

正态性检验

方法：Shapiro-Wilk(n<50)、K-S(n≥50)。参数检验的前提。

正态性检验判断数据是否服从正态分布，是t检验、ANOVA、线性回归的前提。Shapiro-Wilk适用于小样本(n<50)，K-S适用于大样本。大样本时微小偏离也会显著，建议结合Q-Q图和偏度/峰度综合判断。

分析 -> 描述统计 -> 探索
因变量列表选入待检验变量
点击绘制，勾选带检验的正态图
p>0.05服从正态分布；结合Q-Q图辅助判断

# 正态性检验 set.seed(123) data <- rnorm(50, 100, 15) shapiro.test(data) # 偏态数据对比 rexp(50, 0.1) |> shapiro.test() # Q-Q图 qqnorm(data); qqline(data, col="red")

Python

from scipy import stats data = np.random.normal(100, 15, 50) stat, p = stats.shapiro(data) print(f'W={stat:.3f}, p={p:.4f}') import matplotlib.pyplot as plt stats.probplot(data, dist="norm", plot=plt); plt.show()

📋 正态性案例

数据：50名患者收缩压
结果：W=0.982, P=0.234
解读：P>0.05，不拒绝正态性。如不满足，可用非参数检验

方差齐性检验

方法：Levene检验、Bartlett检验。t检验和ANOVA的前提。

方差齐性指各组总体方差相等。Levene检验对非正态数据较稳健；Bartlett对正态敏感。ANOVA方差不齐时可用Welch校正或非参数替代。

分析 -> 比较均值 -> 独立样本T检验
在结果中查看Levene方差齐性检验
P>0.05示方差齐，看第一行t结果
ANOVA中: 在选项勾选方差齐性检验

# 方差齐性 set.seed(123) g1 <- rnorm(30,100,10); g2 <- rnorm(30,110,15); g3 <- rnorm(30,105,12) var.test(g1,g2) library(car) leveneTest(c(g1,g2,g3), factor(rep(1:3,each=30))) t.test(g1,g2,var.equal=FALSE) # Welch

Python

from scipy import stats g1 = np.random.normal(100,10,30) g2 = np.random.normal(110,15,30) g3 = np.random.normal(105,12,30) stat, p = stats.levene(g1, g2, g3) print(f'Levene: stat={stat:.3f}, p={p:.4f}')

📋 方差齐性案例

结果：Levene F=2.34, P=0.102
解读：P>0.05，方差齐，满足ANOVA前提。

独立样本t检验

用途：比较两组独立样本的均值差异。前提：正态性+方差齐。

独立样本t检验比较两组独立(不配对)样本的均值差异。要求数据近似正态且方差齐。不满足可用Welch校正或Mann-Whitney U检验。

分析 -> 比较均值 -> 独立样本T检验
检验变量和分组变量选入
定义组别名称(如 1=试验组, 2=对照组)
查看Levene检验判断方差齐性，阅读对应行t检验结果

# 独立样本t检验 set.seed(123) treatment <- rnorm(35, 125, 10) control <- rnorm(35, 135, 12) t_result <- t.test(treatment, control, var.equal=TRUE) print(t_result) library(effsize); cohen.d(treatment, control)

Python

from scipy import stats treatment = np.random.normal(125,10,35) control = np.random.normal(135,12,35) t, p = stats.ttest_ind(treatment, control) print(f't={t:.3f}, p={p:.4f}') d = (np.mean(treatment)-np.mean(control)) / np.sqrt( (np.var(treatment,ddof=1)+np.var(control,ddof=1))/2) print(f"Cohen's d={d:.3f}")

📋 独立t检验案例

结果：试验组125.3±10.2 vs 对照组135.8±12.1
t=-3.87, df=68, P<0.001, 95%CI[-15.8,-5.2]
解读：新药显著降低SBP，Cohen d=0.93(大效应)。

配对t检验

用途：比较配对样本的均值差异(如治疗前后、配对病例对照)。

配对t检验用于同一组受试者在两个时间点或匹配条件下的均值比较。本质是对差值进行单样本t检验。优点: 消除个体间变异，统计效能更高。

分析 -> 比较均值 -> 成对样本T检验
将治疗前和治疗后变量选入成对变量
查看均值差、t值、df、P值、95%CI

# 配对t检验 set.seed(123) before <- rnorm(25, 140, 15) after <- before - rnorm(25, 10, 5) t.test(before, after, paired=TRUE) # 效应量 diff <- before - after cat("Cohen's d:", mean(diff)/sd(diff), "\n")

Python

from scipy import stats before = np.random.normal(140,15,25) after = before - np.random.normal(10,5,25) t, p = stats.ttest_rel(before, after) print(f'Paired t: t={t:.3f}, p={p:.4f}') d = np.mean(before-after)/np.std(before-after, ddof=1) print(f"d={d:.3f}")

📋 配对t案例

结果：治疗前142.3→131.8
均值差=10.5, t=6.82, df=24, P<0.001

单因素ANOVA

用途：比较三组及以上独立样本的均值差异。事后检验：Tukey、LSD、Bonferroni。

ANOVA用于比较三个或以上组的均值差异。F统计量=组间方差/组内方差。ANOVA显著只表明至少有一组不同，需要事后检验进行两两比较。

分析 -> 比较均值 -> 单因素ANOVA
因变量和因子选入
事后多重比较：方差齐选Tukey，不齐选Games-Howell
选项 -> 描述性 + 方差齐性检验
读取F值和P值，再看事后比较

# ANOVA set.seed(123) value <- c(rnorm(25,100,10), rnorm(25,110,10), rnorm(25,95,10)) group <- factor(rep(c("A","B","C"), each=25)) aov_model <- aov(value ~ group) summary(aov_model) TukeyHSD(aov_model) library(lsr); etaSquared(aov_model)

Python

import pandas as pd import statsmodels.api as sm from statsmodels.formula.api import ols df = pd.DataFrame({ 'value': np.concatenate([np.random.normal(100,10,25), np.random.normal(110,10,25), np.random.normal(95,10,25)]), 'group': ['A']*25+['B']*25+['C']*25 }) print(sm.stats.anova_lm(ols('value~group',data=df).fit(), typ=2)) from statsmodels.stats.multicomp import pairwise_tukeyhsd print(pairwise_tukeyhsd(df['value'], df['group']))

📋 ANOVA案例

结果：F=4.82, df=2,72, P=0.011
事后：A vs C差异显著(P=0.008)。

重复测量ANOVA

用途：同一组受试者多个时间点的测量比较。

重复测量ANOVA用于同一组受试者在三个或以上时间点的测量数据。需满足球形假设(Mauchly检验)。不满足时用Greenhouse-Geisser校正。

# 重复测量ANOVA library(ez) set.seed(123) id <- factor(rep(1:20, each=3)) time <- factor(rep(c("T0","T1","T2"), 20)) value <- c(rnorm(20,140,10), rnorm(20,130,10), rnorm(20,125,10)) rm_data <- data.frame(id, time, value) ezANOVA(data=rm_data, dv=value, wid=id, within=time, detailed=TRUE)

Python

import pandas as pd from statsmodels.stats.anova import AnovaRM df = pd.DataFrame({ 'id': np.repeat(range(20), 3), 'time': np.tile(['T0','T1','T2'], 20), 'value': np.concatenate([np.random.normal(140,10,20), np.random.normal(130,10,20), np.random.normal(125,10,20)]) }) print(AnovaRM(df, 'value', 'id', within=['time']).fit())

📋 重复测量案例

结果：F=12.45, P<0.001
解读：血压随时间显著下降(140→130→125mmHg)。

Mann-Whitney U检验

用途：两组独立样本的非参数检验。t检验的非参数替代，不要求正态性。

Mann-Whitney U检验是独立t检验的非参数替代。基于秩次而非原始值，不要求正态分布。适用于有序分类变量或偏态连续变量。

分析 -> 非参数检验 -> 旧对话框 -> 2个独立样本
检验变量和分组变量选入
检验类型勾选Mann-Whitney U
查看U统计量和精确P值

# Mann-Whitney U检验 set.seed(123) g1 <- rexp(25, 0.1); g2 <- rexp(25, 0.05) wilcox.test(g1, g2) cat("Median G1:", median(g1), "\nMedian G2:", median(g2), "\n")

Python

from scipy import stats g1 = np.random.exponential(10, 25) g2 = np.random.exponential(20, 25) stat, p = stats.mannwhitneyu(g1, g2) print(f'U={stat:.0f}, p={p:.4f}')

📋 Mann-Whitney案例

结果：A组中位8天(IQR:5-12)，B组12天(IQR:7-18)
U=185.5, P=0.003

Wilcoxon符号秩检验

用途：配对数据的非参数检验。配对t检验的非参数替代。

Wilcoxon符号秩检验是配对t检验的非参数替代。计算差值的绝对值排序后比较正负秩和。不要求差值正态。

# Wilcoxon符号秩检验 set.seed(123) before <- runif(20, 50, 100) after <- before - runif(20, 5, 20) wilcox.test(before, after, paired=TRUE) cat("Median diff:", median(before-after), "\n")

Python

from scipy import stats before = np.random.uniform(50,100,20) after = before - np.random.uniform(5,20,20) stat, p = stats.wilcoxon(before, after) print(f'W={stat:.0f}, p={p:.4f}')

📋 Wilcoxon案例

结果：治疗前VAS 7(IQR:5-8) → 3(IQR:2-5)
V=66.5, P=0.002

Kruskal-Wallis检验

用途：三组及以上独立样本的非参数ANOVA。事后：Dunn检验。

Kruskal-Wallis检验是ANOVA的非参数替代。基于秩次比较分布差异。显著后需Dunn检验进行事后两两比较。

# Kruskal-Wallis set.seed(123) value <- c(rexp(20,0.1), rexp(20,0.08), rexp(20,0.05)) group <- factor(rep(c("A","B","C"), each=20)) kruskal.test(value ~ group) library(FSA); dunnTest(value ~ group, method="bonferroni")

Python

from scipy import stats g1, g2, g3 = [np.random.exponential(r,20) for r in [10,12.5,20]] stat, p = stats.kruskal(g1, g2, g3) print(f'H={stat:.3f}, p={p:.4f}')

📋 Kruskal-Wallis案例

结果：H=8.92, df=2, P=0.012
三种化疗方案的恶心评分存在显著差异。

卡方检验

用途：检验分类变量间关联性。条件：期望频数≥5。

卡方检验检验两个分类变量间的关联性。基本思想是比较观测频数与期望频数的差异。公式χ²=Σ(O-E)²/E。条件: 总样本≥40且所有期望频数≥5。

分析 -> 描述统计 -> 交叉表
行变量和列变量选入
统计 -> 卡方
单元格 -> 期望频数和百分比
查看Pearson卡方值和P值

# 卡方检验 data <- matrix(c(45,30,20,50), nrow=2, dimnames=list(Treatment=c("Drug","Placebo"), Outcome=c("有效","无效"))) chi <- chisq.test(data) print(chi) cat("chi2 =", round(chi$statistic,3), ", P =", round(chi$p.value,4)) library(rstatix); cramer_v(data)

Python

from scipy.stats import chi2_contingency table = np.array([[45,30],[20,50]]) chi2, p, dof, exp = chi2_contingency(table) print(f'chi2={chi2:.3f}, p={p:.4f}') n = table.sum() v = np.sqrt(chi2/(n*(min(table.shape)-1))) print(f"V={v:.3f}")

📋 卡方案例

结果：χ²=12.45, df=1, P<0.001, V=0.35
吸烟组肺癌比例显著高于非吸烟组。

Fisher精确检验

用途：小样本或期望频数<5时分类变量的关联性检验。

Fisher精确检验适用于2×2表期望频数<5时。基于超几何分布精确计算P值。提供OR值和精确置信区间。

# Fisher data <- matrix(c(4,10,12,25), nrow=2) f <- fisher.test(data) print(f) cat("OR =", round(f$estimate,3), ", P =", round(f$p.value,4))

Python

from scipy.stats import fisher_exact table = np.array([[4,10],[12,25]]) or_, p = fisher_exact(table) print(f'OR={or_:.3f}, p={p:.4f}')

📋 Fisher案例

结果：P=0.039, OR=8.33
治疗组有效比例显著高于对照组，但样本量小。

线性回归

用途：多个自变量对连续因变量的影响。前提：线性、独立、正态、方差齐。

线性回归建模连续因变量与自变量的线性关系。Y=β₀+β₁X₁+...+ε。核心前提: 线性关系、残差独立、残差正态、方差齐。R²解释模型拟合度。

分析 -> 回归 -> 线性
因变量和自变量选入
统计 -> 估计、模型拟合、Durbin-Watson
绘制 -> 残差图检查方差齐性

# 线性回归 set.seed(123) data <- data.frame(age=rnorm(100,60,10), bmi=rnorm(100,25,3)) data$sbp <- 80 + 0.5*data$age + 1.2*data$bmi + rnorm(100,0,8) lm_model <- lm(sbp ~ age + bmi, data=data) summary(lm_model) confint(lm_model) par(mfrow=c(2,2)); plot(lm_model)

Python

import statsmodels.api as sm X = sm.add_constant(df[['age','bmi']]) model = sm.OLS(df['sbp'], X).fit() print(model.summary()) print(f'R2 = {model.rsquared:.3f}')

📋 线性回归案例

结果：年龄β=0.48(P<0.001), BMI β=1.15(P=0.002)
R²=0.42, 调整R²=0.41
年龄和BMI解释SBP变异的42%。

Logistic回归

用途：二分类因变量建模(发病/未发病)。输出OR值。

Logistic回归用于二分类因变量建模。自变量对log-odds线性影响，指数化转为OR。最大似然估计参数。似然比检验用于模型比较。

分析 -> 回归 -> 二元Logistic
因变量和协变量选入
分类变量点击分类指定
选项 -> HL检验、CI for exp(B)
读取B、OR、95%CI、P值

# Logistic set.seed(123) age <- rnorm(200,60,10); bmi <- rnorm(200,25,3) smoking <- rbinom(200,1,0.3) logit <- -3+0.05*age+0.08*bmi+0.8*smoking disease <- rbinom(200,1,plogis(logit)) model <- glm(disease~age+bmi+smoking, family=binomial()) summary(model) OR <- exp(cbind(OR=coef(model), confint(model))) print(round(OR,3))

Python

X = sm.add_constant(df[['age','bmi','smoking']]) model = sm.Logit(df['disease'], X).fit() print('OR:', np.exp(model.params)) print('95%CI:', np.exp(model.conf_int()))

📋 Logistic案例

结果：吸烟 OR=2.24(95%CI:1.42-3.56, P<0.001)
解读：调整年龄和BMI后，吸烟者冠心病风险是非吸烟者的2.24倍。

Poisson回归

用途：计数数据建模(发病率、死亡率)。偏移量处理人年数据。

Poisson回归用于计数数据建模。连接函数为log，系数指数化为发生率比(IRR)。要求均值=方差(等离散)。过离散时用负二项回归。

# Poisson set.seed(123) age <- rnorm(100,50,12); smoking <- rbinom(100,1,0.3) pyears <- rep(1000,100) lam <- exp(-3+0.03*age+0.5*smoking) events <- rpois(100, lam*pyears/1000) model <- glm(events~age+smoking+offset(log(pyears)), family=poisson()) summary(model) IRR <- exp(cbind(IRR=coef(model), confint(model))) print(round(IRR,3))

Python

model = sm.GLM(df['events'], X, family=sm.families.Poisson()).fit() print('IRR:', np.exp(model.params))

📋 Poisson案例

结果：吸烟 IRR=1.85(95%CI:1.42-2.41, P<0.001)
解读：调整年龄后，吸烟者肺癌发病率是非吸烟者的1.85倍。

Pearson相关

用途：两个连续变量的线性相关。前提：双变量正态。

Pearson相关系数r衡量线性相关强度。|r|<0.3弱，0.3-0.7中等，>0.7强相关。前提: 两变量近似正态且关系呈线性。

分析 -> 相关 -> 双变量
变量选入，勾选Pearson
查看r值和P值

# Pearson set.seed(123) x <- rnorm(50,100,15); y <- 0.6*x + rnorm(50,0,10) cor.test(x, y, method="pearson") plot(x,y,col="#0d9488"); abline(lm(y~x),col="red")

Python

from scipy import stats r, p = stats.pearsonr(x, y) print(f'r={r:.3f}, p={p:.4f}') z = np.arctanh(r); se = 1/np.sqrt(len(x)-3) ci = np.tanh(z + np.array([-1,1])*1.96*se) print(f'95%CI: [{ci[0]:.3f}, {ci[1]:.3f}]')

📋 Pearson案例

结果：r=0.52, P<0.001, 95%CI[0.35,0.68]
BMI与SBP呈中等正相关。

Spearman相关

用途：非参数相关。适用于偏态或有序数据。

Spearman秩相关是Pearson的非参数替代。基于秩次，衡量单调相关。对异常值不敏感，适用于有序变量或偏态分布。

# Spearman set.seed(123) x <- rexp(30,0.1); y <- x^2 + runif(30,0,50) cor.test(x, y, method="spearman")

Python

from scipy import stats rho, p = stats.spearmanr(x, y) print(f'rho={rho:.3f}, p={p:.4f}')

📋 Spearman案例

结果：ρ=-0.38, P=0.008
住院天数与满意度呈弱负相关。

ICC组内相关系数

用途：评价连续变量测量的一致性/信度。>0.75良好，>0.9优秀。

ICC用于评价连续变量测量一致性。不同类型对应不同设计: ICC(1,1)随机选评价者; ICC(2,1)固定评价者; ICC(3,1)一致性测量。

# ICC library(irr) r1 <- rnorm(20,50,10); r2 <- r1 + rnorm(20,0,3) icc(data.frame(r1,r2), model="twoway", type="consistency") library(psych) multi <- data.frame(r1=rnorm(20,50,10), r2=rnorm(20,50,10)+rnorm(20,0,3), r3=rnorm(20,50,10)+rnorm(20,0,4)) ICC(multi)

Python

import pingouin as pg df = pd.DataFrame({ 'target': np.repeat(range(20),2), 'rater': np.tile(['A','B'],20), 'score': np.concatenate([r1, r2]) }) print(pg.intraclass_corr(df, 'target', 'rater', 'score'))

📋 ICC案例

结果：ICC(2,1)=0.87, 95%CI[0.78-0.93]
评分者间一致性良好。

Kappa一致性

用途：分类变量的一致性评价。加权Kappa用于有序分类。

Kappa系数评价两个评价者对分类变量的一致性，校正了随机一致。Kappa≥0.75一致性好，0.4-0.75中等，<0.4差。加权Kappa用于有序分类。

# Kappa library(irr) r1 <- c(1,1,1,0,1,0,0,0,1,1) r2 <- c(1,1,1,1,1,0,0,0,0,0) kappa2(cbind(r1, r2)) # 加权Kappa ord1 <- c(1,2,3,2,1,2,3,3,2,1) ord2 <- c(1,2,3,2,2,2,3,2,2,1) kappa2(cbind(ord1, ord2), weight="squared")

Python

from sklearn.metrics import cohen_kappa_score r1 = [1,1,1,0,1,0,0,0,1,1] r2 = [1,1,1,1,1,0,0,0,0,0] kappa = cohen_kappa_score(r1, r2) print(f'Kappa = {kappa:.3f}') wkappa = cohen_kappa_score(r1, r2, weights='linear') print(f'Weighted Kappa = {wkappa:.3f}')

📋 Kappa案例

结果：Kappa=0.78, P<0.001
两名诊断医生的判断一致性良好。

诊断试验指标

核心：灵敏度(Se)、特异度(Sp)、PPV、NPV、似然比(LR)、约登指数。

诊断试验指标评价诊断方法的准确性。Se=TP/(TP+FN)，Sp=TN/(TN+FP)。PPV为阳性预测值，NPV为阴性预测值。LR+=Se/(1-Sp)，LR-=(1-Se)/Sp。约登指数=Se+Sp-1。

# 诊断指标 set.seed(123) gold <- rbinom(200,1,0.3) test <- ifelse(gold==1, rbinom(200,1,0.85), rbinom(200,1,0.15)) TP=sum(test==1&gold==1); FP=sum(test==1&gold==0) FN=sum(test==0&gold==1); TN=sum(test==0&gold==0) cat("Se:", round(TP/(TP+FN),3), "\nSp:", round(TN/(TN+FP),3)) cat("PPV:", round(TP/(TP+FP),3), "\nNPV:", round(TN/(TN+FN),3)) cat("LR+:", round((TP/(TP+FN))/(1-TN/(TN+FP)),2)) cat("\nLR-:", round((1-TP/(TP+FN))/(TN/(TN+FP)),2)) library(epiR); epi.tests(data.frame(TP,FP,FN,TN))

Python

from sklearn.metrics import confusion_matrix gold = np.random.binomial(1,0.3,200) test = np.where(gold==1, np.random.binomial(1,0.85,200), np.random.binomial(1,0.15,200)) cm = confusion_matrix(gold, test) TP, FP, FN, TN = cm[1,1], cm[0,1], cm[1,0], cm[0,0] se=TP/(TP+FN); sp=TN/(TN+FP) ppv=TP/(TP+FP); npv=TN/(TN+FN) print(f'Se={se:.3f}, Sp={sp:.3f}, PPV={ppv:.3f}, NPV={npv:.3f}') print(f'LR+={se/(1-sp):.2f}, LR-={(1-se)/sp:.2f}')

📋 诊断指标案例

结果：Se=0.85, Sp=0.85, LR+=5.67
PPV=0.71, NPV=0.93
阴性预测值高，阴性结果可较好排除疾病。

ROC曲线

用途：评价诊断试验/预测模型的判别能力。AUC、最佳切点、Se、Sp。

ROC曲线以1-特异度为横轴、灵敏度为纵轴绘制。AUC综合评价: 0.5无判别，0.7-0.8可接受，0.8-0.9优秀。最佳切点用约登指数确定。DeLong检验比较两个AUC。

分析 -> ROC曲线
检验变量选入预测概率，状态变量选入金标准
状态变量值设为1(阳性)
选项 -> AUC置信区间和坐标点

# ROC library(pROC) set.seed(123) score <- c(rnorm(50,70,15), rnorm(50,50,15)) true <- c(rep(1,50), rep(0,50)) roc_obj <- roc(true, score) cat("AUC:", round(auc(roc_obj),3), "\n95%CI:", round(ci.auc(roc_obj),3)) best <- coords(roc_obj, "best", best.method="youden") print(best) plot(roc_obj, col="#0d9488", lwd=2)

Python

from sklearn.metrics import roc_curve, auc fpr, tpr, thr = roc_curve(true, score) roc_auc = auc(fpr, tpr) print(f'AUC={roc_auc:.3f}') youden = tpr - fpr best = thr[np.argmax(youden)] print(f'Best threshold={best:.2f}')

📋 ROC案例

结果：AUC=0.87, 95%CI[0.82-0.92]
最佳切点NT-proBNP>450pg/mL: Se=0.82, Sp=0.79

校准曲线

用途：评价预测模型的校准度—预测概率与实际概率的一致性。

校准曲线展示预测概率与实际观测概率的一致性。完美校准沿45°对角线。HL检验不显著=校准良好。Brier Score综合评估判别和校准。

# 校准曲线 library(rms) set.seed(123) n <- 300 data <- data.frame(age=rnorm(n,60,10), bmi=rnorm(n,25,3), smoking=rbinom(n,1,0.3)) logit <- -3+0.05*data$age+0.08*data$bmi+0.8*data$smoking data$disease <- rbinom(n, 1, plogis(logit)) dd <- datadist(data); options(datadist="dd") model <- lrm(disease~age+bmi+smoking, data=data, x=TRUE, y=TRUE) cal <- calibrate(model, B=200) plot(cal, xlab="Predicted", ylab="Actual")

Python

from sklearn.calibration import calibration_curve import matplotlib.pyplot as plt prob_pred, prob_true = calibration_curve(y_true, y_pred, n_bins=10) plt.plot(prob_pred, prob_true, marker='o') plt.plot([0,1],[0,1], 'k--'); plt.show() from sklearn.metrics import brier_score_loss print(f'Brier Score = {brier_score_loss(y_true, y_pred):.4f}')

📋 校准案例

结果：校准曲线接近45°对角线
HL检验P=0.342(>0.05，校准良好)
Brier Score=0.12<0.25(有临床价值)

PCA主成分分析

用途：降维、多重共线性处理、高维数据可视化。

PCA将多个相关变量转为不相关主成分(原始变量的线性组合)。第一主成分方差最大。需对变量标准化。用于降维、去共线性、探索数据结构。

# PCA set.seed(123) data <- data.frame(x1=rnorm(100,50,10), x2=0.7*x1+rnorm(100,0,5), x3=0.5*x1+0.3*x2+rnorm(100,0,4), x4=rnorm(100,30,8)) pca <- prcomp(data, scale=TRUE) summary(pca) print(pca$rotation[,1:2]) screeplot(pca, type="lines"); abline(h=1,lty=2,col="red") biplot(pca)

Python

from sklearn.decomposition import PCA from sklearn.preprocessing import StandardScaler pca = PCA(n_components=2) X_pca = pca.fit_transform(StandardScaler().fit_transform(df)) print(f'Explained variance: {pca.explained_variance_ratio_}') print(f'Cumulative: {np.cumsum(pca.explained_variance_ratio_)}')

📋 PCA案例

结果：前2个主成分解释68%方差
PC1=代谢负荷(BMI+血糖+血脂), PC2=炎症因子

EFA探索性因子分析

用途：探索变量背后的潜在结构/因子。旋转：Varimax(正交)、Promax(斜交)。

EFA假设观测变量由潜在因子决定，用于问卷结构效度评价。与PCA不同，因子分析有测量误差模型。旋转使因子载荷更清晰。KMO>0.6且Bartlett P<0.05表示适合做EFA。

# EFA library(psych) set.seed(123) data <- data.frame(q1=rnorm(200,3,1), q2=0.6*q1+rnorm(200,0,0.8), q3=0.7*q1+rnorm(200,0,0.7), q4=rnorm(200,3,1), q5=0.5*q4+rnorm(200,0,0.9), q6=0.6*q4+rnorm(200,0,0.8)) KMO(cor(data)) cortest.bartlett(cor(data), n=200) fa <- fa(data, nfactors=2, rotate="varimax") print(fa$loadings, cutoff=0.3)

Python

from factor_analyzer import FactorAnalyzer, calculate_kmo from factor_analyzer.factor_analyzer import calculate_bartlett_sphericity kmo_all, kmo = calculate_kmo(df) print(f'KMO = {kmo:.3f}') chi2, p = calculate_bartlett_sphericity(df) print(f'Bartlett: p={p:.4f}') fa = FactorAnalyzer(n_factors=2, rotation='varimax') fa.fit(df) print('Loadings:', fa.loadings_)

📋 EFA案例

结果：KMO=0.82, Bartlett P<0.001
2个因子解释58%方差: 因子1=躯体症状, 因子2=心理症状
问卷结构效度良好。

K-means聚类

用途：将样本划分为K个簇。基于距离的硬聚类。

K-means将数据划分为K个簇，每个样本属于最近的质心。需指定K值。最佳K可通过肘部法则(降低曲线)或轮廓系数选择。

# K-means set.seed(123) data <- data.frame(x1=rnorm(150,c(0,5,10),1.5), x2=rnorm(150,c(0,5,10),1.5)) wss <- sapply(1:8, function(k) kmeans(data,k,nstart=25)$tot.withinss) plot(1:8,wss,type="b",xlab="K",ylab="WSS") km <- kmeans(data, centers=3, nstart=25) plot(data$x1, data$x2, col=km$cluster, pch=19) points(km$centers, col=1:3, pch=8, cex=2)

Python

from sklearn.cluster import KMeans from sklearn.metrics import silhouette_score kmeans = KMeans(n_clusters=3, random_state=123, n_init=25) df['cluster'] = kmeans.fit_predict(X) sil = silhouette_score(X, df['cluster']) print(f'Silhouette = {sil:.3f}')

📋 K-means案例

结果：K=3最优(轮廓系数0.62)
簇1: 年轻+轻度代谢异常 | 簇2: 中年+中重度 | 簇3: 老年+多并发症

层次聚类

用途：无需预设K的聚类。输出树状图。

层次聚类按层次聚合(自底向上)或分裂(自顶向下)。不需预指定K，树状图直观展示聚类层次。常见连接法: ward.D2、complete、average。

# 层次聚类 set.seed(123) data <- data.frame(x1=rnorm(50,c(0,5,10),1.5), x2=rnorm(50,c(0,5,10),1.5)) dist_m <- dist(data, method="euclidean") hc <- hclust(dist_m, method="ward.D2") plot(hc, main="Dendrogram", xlab="", sub="") rect.hclust(hc, k=3, border=1:3) groups <- cutree(hc, k=3) print(table(groups))

Python

from scipy.cluster.hierarchy import dendrogram, linkage, fcluster Z = linkage(df_scaled, method='ward') dendrogram(Z); plt.show() clusters = fcluster(Z, t=3, criterion='maxclust') print(pd.Series(clusters).value_counts())

📋 层次聚类案例

结果：树状图显示3个主要腐瘤亚型，与临床预后相关
发现新的疾病亚型。

中介分析

用途：检验X通过M影响Y的间接路径(中介效应)。

中介分析检验X是否通过M影响Y。总效应=直接效应(c')+间接效应(a×b)。Bootstrap法检验间接效应(不要求正态)。

# 中介分析 library(mediation) set.seed(123) n <- 200 data <- data.frame(X=rnorm(n), M=0.5*X+rnorm(n), Y=0.3*M+0.2*X+rnorm(n)) model_M <- lm(M ~ X, data=data) model_Y <- lm(Y ~ X + M, data=data) med <- mediate(model_M, model_Y, treat="X", mediator="M", boot=TRUE, sims=1000) summary(med) cat("ACME:", med$d0, "\nDirect:", med$z0, "\nProp:", med$n.avg)

Python

a_m = sm.OLS(df['M'], sm.add_constant(df['X'])).fit() b_m = sm.OLS(df['Y'], sm.add_constant(df[['X','M']])).fit() indirect = a_m.params['X'] * b_m.params['M'] direct = b_m.params['X'] total = indirect + direct print(f'Indirect={indirect:.3f}, Direct={direct:.3f}') print(f'Mediation ratio={indirect/total:.1%}')

📋 中介案例

结果：间接效应=0.12(P=0.002), 直接效应=0.08(P=0.124)
解读：健康素养完全中介SES对健康的影响(中介比例~60%)。

调节效应

用途：检验调节变量W是否改变X对Y的影响方向/强度(交互作用)。

调节效应检验X对Y的效应是否依赖于调节变量W。核心是交互项X×W的回归系数。交互项显著即存在调节效应。连续变量需中心化。

# 调节效应 set.seed(123) n <- 200; X <- rnorm(n); W <- rnorm(n) Y <- 0.3*X + 0.2*W + 0.4*X*W + rnorm(n) data <- data.frame(X, W, Y) data$X_c <- scale(data$X, scale=FALSE) data$W_c <- scale(data$W, scale=FALSE) data$XW <- data$X_c * data$W_c mod <- lm(Y ~ X_c + W_c + XW, data=data) summary(mod) library(interactions); sim_slopes(mod, pred=X_c, modx=W_c)

Python

df['X_c'] = df['X']-df['X'].mean() df['W_c'] = df['W']-df['W'].mean() df['XW'] = df['X_c']*df['W_c'] X = sm.add_constant(df[['X_c','W_c','XW']]) model = sm.OLS(df['Y'], X).fit() print(model.summary()) # 交互项显著→存在调节效应

📋 调节效应案例

结果：运动×年龄交互项β=0.03, P=0.008
解读：年龄调节了运动对血压的效果—年轻人运动降压效果更明显。

⏳ survival R包文档

⏳ 生存分析

Kaplan-Meier分析

生存曲线估计，组间比较

library(survival) library(survminer) # 使用内置肺癌数据 data(lung) lung$sex <- factor(lung$sex, labels=c("男","女")) # KM估计 km_fit <- survfit(Surv(time, status) ~ sex, data=lung) summary(km_fit) # Log-rank检验 survdiff(Surv(time, status) ~ sex, data=lung) # 生存曲线 ggsurvplot(km_fit, data=lung, pval=TRUE, conf.int=TRUE, risk.table=TRUE, xlab="时间(天)", ylab="生存概率")

Cox回归

多因素生存分析

library(survival) # Cox回归 cox <- coxph(Surv(time, status) ~ sex + age + ph.ecog, data=lung) summary(cox) # PH假设检验 cox.zph(cox) # 森林图 library(survminer) ggforest(cox, data=lung) # 校准曲线 library(riskRegression) score <- Score(list("Cox"=cox), formula=Surv(time,status)~1, data=lung, times=365) plotCalibration(score, times=365)

竞争风险模型

Fine-Gray模型，累积发生率

library(cmprsk) # 模拟竞争风险数据 set.seed(123) time <- rexp(200, rate=0.1) status <- sample(0:2, 200, replace=TRUE, prob=c(0.3,0.5,0.2)) group <- sample(0:1, 200, replace=TRUE) # 累积发生率 ci <- cuminc(time, status, group) plot(ci, xlab="时间", ylab="累积发生率") # Fine-Gray检验 print(ci$Tests) # 竞争风险回归 fg <- crr(time, status, covariate=as.matrix(group)) summary(fg)

🔧 缺失数据处理

缺失机制

MCAR：完全随机缺失，与数据值无关
MAR：随机缺失，与其他观测变量有关
MNAR：非随机缺失，与缺失值本身有关

多重插补MICE

推荐方法

library(mice) set.seed(123) data <- data.frame(age=rnorm(100,60,10), bmi=rnorm(100,25,3), sbp=rnorm(100,130,15)) data$bmi[sample(1:100,20)] <- NA data$sbp[sample(1:100,15)] <- NA md.pattern(data) imp <- mice(data, m=5, method="pmm", seed=123) complete_data <- complete(imp,1) fit <- with(imp, lm(sbp~age+bmi)) summary(pool(fit))

Python

from sklearn.impute import IterativeImputer import pandas as pd, numpy as np imp = IterativeImputer(max_iter=10, random_state=123) df_imp = pd.DataFrame(imp.fit_transform(df), columns=df.columns)

📈 pROC R包

📈 ROC曲线与预测模型评价

ROC曲线基础

AUC：0.5=无区分，0.7-0.8可接受，0.8+优秀
最佳切点：约登指数(灵敏度+特异度-1)
多模型比较：DeLong检验
校准曲线：预测概率vs实际概率

ROC比较与校准

多模型AUC对比

library(pROC) # 比较两个模型AUC set.seed(123) true <- rbinom(100, 1, 0.5) pred1 <- rnorm(100) + true*1.5 pred2 <- rnorm(100) + true*1.0 roc1 <- roc(true, pred1) roc2 <- roc(true, pred2) # DeLong检验 roc.test(roc1, roc2, method="delong") # 校准曲线 library(rms) cal <- calibrate(lrm(true ~ pred1), B=200) plot(cal) # DCA决策曲线 library(rmda) dca <- decision_curve(true ~ pred1, data=data.frame(true, pred1), bootstraps=50) plot_decision_curve(dca)

📈 时间序列分析

ARIMA模型

自回归移动平均，预测未来趋势

library(forecast) set.seed(123) ts_data <- ts(rpois(120,10)+5*sin(2*pi*1:120/12), start=c(2020,1), frequency=12) model <- auto.arima(ts_data) fc <- forecast(model, h=12) plot(fc) accuracy(fc)

Python

from statsmodels.tsa.arima.model import ARIMA import numpy as np, pandas as pd ts = np.random.poisson(10,120)+5*np.sin(2*np.pi*np.arange(120)/12) model = ARIMA(ts, order=(1,1,1)).fit() fc = model.forecast(steps=12) print(fc)

ETS指数平滑

用途：处理时间序列的趋势和季节性，适合单变量预测。

ETS（Error Trend Seasonal）是指数平滑方法的统一框架，用于处理时间序列的误差、趋势和季节性。ETS(A,N,N)表示可加性误差、无趋势、无季节；ETS(A,A,A)表示可加性误差、可加性趋势、可加性季节。通过拟合优度选择最优模型。

library(forecast) set.seed(123) ts_data <- ts(rpois(120,10)+5*sin(2*pi*1:120/12), start=c(2020,1), frequency=12) fit <- ets(ts_data) summary(fit) fc <- forecast(fit, h=12) plot(fc) accuracy(fc)

Python

import numpy as np, pandas as pd from statsmodels.tsa.holtwinters import ExponentialSmoothing np.random.seed(123) ts = np.random.poisson(10,120)+5*np.sin(2*np.pi*np.arange(120)/12) df = pd.DataFrame({'y':ts}, index=pd.date_range('2020-01',periods=120,freq='M')) model = ExponentialSmoothing(df['y'], seasonal_periods=12, trend='add', seasonal='add').fit() fc = model.forecast(12) print(fc)

📋 ETS案例

数据：2019-2024月度呼吸道感染病登记例数
结果：ETS(A,N,A)模型最优，RMSE=15.3
解读：去季节后套合良好，2025年预测值呈现稳定季节性波动。

Prophet时间序列预测

用途：Facebook开源的时间序列预测框架，自动处理节假日效应和变点。

Prophet由Facebook于2017年开源，基于可分解时间序列模型：趋势 + 季节性 + 节假日效应。Prophet能自动处理缺失值、异常值和变点，且对非专业用户友好。适合具有明显季节性和多年历史数据的预测场景。

library(prophet) library(dplyr) # 模拟月度数据 set.seed(123) ds <- seq(as.Date("2020-01-01"), as.Date("2024-12-01"), by="month") trend <- 10 + 0.05*1:60 seasonal <- 3*sin(2*pi*1:60/12) y <- trend + seasonal + rnorm(60, 0, 1) df <- data.frame(ds, y) # Prophet模型 m <- prophet(df, yearly.seasonality=TRUE) future <- make_future_dataframe(m, periods=12, freq="month") fc <- predict(m, future) plot(m, fc) prophet_plot_components(m, fc)

Python

import numpy as np, pandas as pd from prophet import Prophet np.random.seed(123) ds = pd.date_range("2020-01-01", periods=60, freq="M") trend = 10 + 0.05*np.arange(60) seasonal = 3*np.sin(2*np.pi*np.arange(60)/12) y = trend + seasonal + np.random.normal(0, 1, 60) df = pd.DataFrame({"ds": ds, "y": y}) model = Prophet(yearly_seasonality=True) model.fit(df) future = model.make_future_dataframe(periods=12, freq="M") fc = model.predict(future) model.plot(fc).show() model.plot_components(fc).show()

📋 Prophet案例

数据：2020-2024年某医院月度门诊量
结果：模型拟合良好，周期性波动被成功捕捉
解读：2025年门诊量预计保持上升趋势，季节性高峰出现在冬季呼吸道疾病高发期。

LSTM时间序列预测

用途：深度学习方法，适合复杂模式和长期依赖的时间序列。

LSTM（Long Short-Term Memory）是循环神经网络的变种，通过徘徊门机制解决了传统RNN的梯度消失问题。LSTM能捕捉时间序列中的长期依赖关系，适合复杂模式的时间序列预测。需要注意：需要先处理数据创建观测窗口，通常比统计方法需要更多数据。

library(keras) # 生成模拟数据 set.seed(123) n <- 300 x <- seq(0, 10*pi, length.out=n) y <- sin(x) + rnorm(n, 0, 0.1) # 创建时间窗口 create_seq <- function(data, lookback=10){ X <- t(sapply(lookback:n, function(i) data[(i-lookback+1):i])) Y <- data[(lookback+1):(n+1)] list(X=array(X, dim=c(nrow(X), lookback, 1)), Y=Y) } lookback <- 10 data <- create_seq(y, lookback) # LSTM模型 model <- keras_model_sequential() |> layer_lstm(units=50, input_shape=c(lookback,1)) |> layer_dense(units=1) model |> compile(optimizer='adam', loss='mse') history <- model |> fit(data$X[1:200,,], data$Y[1:200], epochs=50, batch_size=16, validation_split=0.2, verbose=0) # 预测 pred <- model |> predict(data$X[201:290,,]) plot(data$Y[201:290], type='l', col='blue') lines(pred, col='red') legend('topright', c('True','Pred'), col=c('blue','red'), lty=1)

Python

import numpy as np from tensorflow.keras.models import Sequential from tensorflow.keras.layers import LSTM, Dense np.random.seed(123) n = 300; x = np.linspace(0, 10*np.pi, n) y = np.sin(x) + np.random.normal(0, 0.1, n) lookback = 10 X, Y = [], [] for i in range(lookback, n-1): X.append(y[i-lookback+1:i+1]) Y.append(y[i+1]) X = np.array(X).reshape(-1, lookback, 1) Y = np.array(Y) # LSTM模型 model = Sequential([ LSTM(50, input_shape=(lookback, 1)), Dense(1) ]) model.compile(optimizer="adam", loss="mse") model.fit(X[:200], Y[:200], epochs=50, batch_size=16, validation_split=0.2, verbose=0) pred = model.predict(X[200:290], verbose=0) print("Prediction shape:", pred.shape)

📋 LSTM案例

数据：模拟的周期性生物标志物测定时间序列
结果：LSTM成功学习了波形模式，MSE=0.015
解读：深度学习方法在捕捉复杂循环模式方面优于传统统计方法，但需要较多数据和调参。

VAR向量自回归

用途：多变量时间序列分析，同时预测多个相互关联的变量。

VAR（Vector Autoregression）是多变量时间序列的经典模型，同时建模多个时间序列的动态关系。VAR(p)表示包含P阶滞后项。核心优势：无需区分内生和外生变量，所有变量均等待处。通常用于宏观经济和卫生政策评估。

library(vars) set.seed(123) n <- 200 # 双变量VAR(1)模拟 e1 <- rnorm(n); e2 <- rnorm(n) y1 <- 0.5*e1[1]; y2 <- 0.3*e2[1] for(i in 2:n){ y1[i] <- 1 + 0.6*y1[i-1] + 0.2*y2[i-1] + e1[i] y2[i] <- 2 + 0.1*y1[i-1] + 0.5*y2[i-1] + e2[i] } data <- cbind(y1, y2) # VAR模型 var_select <- VARselect(data, lag.max=8) print(var_select$selection) model <- VAR(data, p=2) summary(model) # 脉冲响应 irf <- irf(model, n.ahead=10) plot(irf) # 预测 fc <- predict(model, n.ahead=6) plot(fc)

Python

import numpy as np, pandas as pd from statsmodels.tsa.api import VAR np.random.seed(123) n = 200; e1 = np.random.normal(0, 1, n) e2 = np.random.normal(0, 1, n) y1 = np.zeros(n); y2 = np.zeros(n) y1[0] = 0.5*e1[0]; y2[0] = 0.3*e2[0] for i in range(1, n): y1[i] = 1 + 0.6*y1[i-1] + 0.2*y2[i-1] + e1[i] y2[i] = 2 + 0.1*y1[i-1] + 0.5*y2[i-1] + e2[i] df = pd.DataFrame({"y1": y1, "y2": y2}) model = VAR(df) lag_order = model.select_order(maxlags=8) print(lag_order.selected_orders) result = model.fit(2) print(result.summary()) # 脉冲响应 irf = result.irf(10) irf.plot() # 预测 fc = result.forecast(df.values[-2:], steps=6) print("Forecast:", fc)

📋 VAR案例

数据：某地区每月流感发病率和气温数据
结果：VAR(2)模型最优，气温对流感有显著滞后影响
解读：气温下降后2周流感发病率上升，为预防提供了时间窗口。

🧮 样本量计算（含R代码）

核心参数

α值：一类错误，通常0.05
β值：二类错误，通常0.2，效能80%
效应量：组间预期差异
脱落率：10–20%

ANOVA样本量

多组均值比较的方差分析样本量

ANOVA样本量计算基于Cohen's f效应量。f = σ_m/σ（组间标准差/组内标准差）。小效应f=0.10，中等f=0.25，大效应f=0.40。需设定组数k、检验效能(通常0.80)、显著性水平α。

library(pwr) # ANOVA样本量 k=3组, 中等效应f=0.25 pwr.anova.test(k=3, f=0.25, power=0.8, sig.level=0.05) # 结果：总样本159例(每组53例)

Python

from statsmodels.stats.power import FTestAnovaPower # ANOVA样本量 k=3组, 中等效应f=0.25 n = FTestAnovaPower().solve_power( effect_size=0.25, power=0.8, alpha=0.05, k_groups=3) print(f"总样本: {n:.0f}例")

📋 ANOVA样本量案例

场景：比较三种降压药的疗效，预期效应量f=0.25
结果：每组需要53例，共159例
解读：中等效应下3组比较需要较大样本量，若效应量更小则样本量需进一步增加。

回归分析样本量

多元线性回归样本量

回归分析样本量常用经验法则：每变量10-20个事件（EPV）。Cohen's f²效应量：小f²=0.02，中等f²=0.15，大f²=0.35。包含预测因子数量p、期望R²、检验效能等参数。推荐同时使用EPV法则和pwr.f2.test进行综合估算。

library(pwr) # 多元回归样本量：5个预测因子 # 中等效应f²=0.15 pwr.f2.test(u=5, f2=0.15, power=0.8, sig.level=0.05) # 结果：v=68.5，总样本n=u+v+1≈74 # EPV经验法则：5变量×10=50例起 # 综合建议取较大值: 至少74例 # 小效应f²=0.02时 pwr.f2.test(u=5, f2=0.02, power=0.8, sig.level=0.05) # 结果：v=646，总样本≈652例

Python

from statsmodels.stats.power import FTestRegPower # 多元回归样本量：5个预测因子, f²=0.15 n = FTestRegPower().solve_power( effect_size=0.15, power=0.8, alpha=0.05, k_groups=1, df_num=5) # df_num = 预测因子数 print(f"所需样本量: {n:.0f}例") # EPV法则 n_vars = 5 epv_min = n_vars * 10 epv_rec = n_vars * 20 print(f"EPV法最少: {epv_min}例，推荐: {epv_rec}例")

📋 回归分析样本量案例

场景：构建5个预测因子的线性回归模型预测术后并发症
结果：中等效应f²=0.15时需要74例；按EPV法则需要50-100例
解读：建议综合两种方法取较大值，同时考虑脱落率10-20%。

生存分析样本量

Log-rank检验生存分析样本量

生存分析样本量基于事件数（不是总人数）。常用Schoenfeld公式：事件数d = (z_α/2+z_β)² / [p(1-p)(log HR)²]。p为分组比例，HR为风险比。总样本量 = 事件数 / 事件率。需考虑随访时间和删失比例。

library(powerSurvEpi) # Cox回归样本量：HR=1.8, 2年随访 # 对照组生存率0.70 ssizeCT.default(power=0.8, alpha=0.05, HR=1.8, pE=0.30, pC=0.70, ratio=1) # pE=暴露组事件率, pC=对照组事件率 # 简易法：Schoenfeld公式 # 事件数需求 = (1.96+0.84)² / [0.5*0.5*(log1.8)²] # ≈ 66个事件，若5年事件率0.5，需132例

Python

import math # Schoenfeld公式计算生存分析样本量 alpha = 0.05 power = 0.80 HR = 1.8 # 风险比 p = 0.5 # 分组比例 z_alpha = 1.96 z_beta = 0.842 # 所需事件数 d = (z_alpha + z_beta)**2 / (p * (1-p) * (math.log(HR)**2)) print(f"所需事件数: {math.ceil(d)}") # 假设5年事件率0.50 event_rate = 0.50 total_n = math.ceil(d / event_rate) print(f"总样本量: {total_n}例")

📋 生存分析样本量案例

场景：评估新疗法对总生存期的影响，预期HR=1.8，等比例分组
结果：需要65个事件，假设5年内50%患者发生事件，需130例
解读：生存分析样本量主要取决于事件数，延长随访可提高事件率从而减少总样本量。

非劣效性样本量

非劣效/等效性试验样本量

非劣效性样本量基于单侧检验。核心参数：非劣效界值Δ（临床上可接受的最大差异）、预期疗效差异、标准差。公式基于单侧z检验：n/组 = 2(z_α+z_β)²(σ/Δ)²。等效性试验需双侧检验并交换α分配。

library(TrialSize) # 非劣效性：两组均值比较 # 非劣效界值delta=3, SD=10, 单侧α=0.025 NonInferiority(alpha=0.025, beta=0.2, sigma=10, delta=3, margin=0) # margin=预期差异 # 结果：每组约175例 # 率的非劣效性 # 对照组率0.85, 非劣效界值-0.10 # 使用单侧z检验近似 p1=0.85; p2=0.85; delta=-0.10 n = ceiling(2*((1.96+0.842)*sqrt(p1*(1-p1)+p2*(1-p2))/abs(delta))^2) print(paste("每组需", n, "例"))

Python

import math # 非劣效性试验样本量（连续变量） alpha = 0.025 # 单侧 beta = 0.20 delta = 3.0 # 非劣效界值 sigma = 10.0 # 标准差 z_alpha = 1.96 z_beta = 0.842 n_per = math.ceil(2 * (z_alpha + z_beta)**2 * (sigma/delta)**2) print(f"每组需: {n_per}例，共{n_per*2}例") # 率的非劣效性 p = 0.85 # 预期率 delta_p = 0.10 # 非劣效界值 n_rate = math.ceil(2 * (z_alpha + z_beta)**2 * (2*p*(1-p)) / delta_p**2) print(f"率的比较每组需: {n_rate}例")

📋 非劣效性样本量案例

场景：仿制药疗效非劣于原研药，非劣效界值Δ=3，SD=10
结果：单侧α=0.025，效能80%，每组需175例，共350例
解读：非劣效性试验通常需要较大样本量，界值Δ越小组所需样本量越大。率的比较所需样本量通常大于连续变量。

t检验样本量

两组均值比较

library(pwr) # t检验样本量：两组均值差=5, SD=10 # 效应量 d = 5/10 = 0.5 pwr.t.test(d=0.5, power=0.8, sig.level=0.05, type="two.sample", alternative="two.sided") # 结果：每组需64例，共128例 # ANOVA样本量 pwr.anova.test(k=3, f=0.25, power=0.8, sig.level=0.05) # f=0.25为中等效应

率的比较样本量

卡方检验样本量

library(pwr) # 率的比较：p1=0.6, p2=0.4 h <- ES.h(p1=0.6, p2=0.4) pwr.2p.test(h=h, power=0.8, sig.level=0.05) # 结果：每组需97例 # Cox回归样本量 # 事件数要求：一般每变量10-20个事件 # 总样本 = 事件数 / 事件率

计算工具

在线统计

样本量

专业计算

临床研究

⚖️ 因果推断

⚖️ MatchIt R包

⚖️ 倾向性评分与因果推断

PSM倾向性评分匹配

用途：通过倾向性得分匹配减少观察性研究的纺成偏误。

PSM（Propensity Score Matching）通过Logistic回归估计每个个体接受处理的概率（倾向性得分），然后在处理组和对照组之间进行匹配。常用匹配方法：最近邻匹配、校准匹配、占比匹配、分数匹配。核心假设：无未观测纺成（强忽略性）。

library(MatchIt); library(cobalt) set.seed(123); n <- 500 data <- data.frame( age = rnorm(n, 60, 10), sex = rbinom(n, 1, 0.5), bmi = rnorm(n, 25, 3), treatment = rbinom(n, 1, 0.3) ) data$outcome <- rbinom(n, 1, plogis( -2 + 0.5*data$treatment + 0.03*data$age)) # 倾向性评分匹配 m.out <- matchit(treatment ~ age + sex + bmi, data=data, method="nearest", ratio=1, caliper=0.2) summary(m.out) # 匹配前后均衡性 love.plot(m.out, binary="std") # 匹配后分析 m.data <- match.data(m.out) glm(outcome ~ treatment, data=m.data, family=binomial(), weights=weights)

Python

import numpy as np, pandas as pd from sklearn.linear_model import LogisticRegression from sklearn.neighbors import NearestNeighbors np.random.seed(123); n = 500 data = pd.DataFrame({ "age": np.random.normal(60, 10, n), "sex": np.random.binomial(1, 0.5, n), "bmi": np.random.normal(25, 3, n), "treatment": np.random.binomial(1, 0.3, n) }) data["outcome"] = np.random.binomial(1, 1/(1+np.exp(-(-2 + 0.5*data["treatment"] + 0.03*data["age"])))) # 计算倾向性得分 ps = LogisticRegression().fit( data[["age","sex","bmi"]], data["treatment"]).predict_proba(data[["age","sex","bmi"]])[:,1] # 最近邻匹配 treated = data[data["treatment"]==1]; control = data[data["treatment"]==0] nn = NearestNeighbors(n_neighbors=1, metric="euclidean") nn.fit(control[["age","sex","bmi"]]) # 匹配后分析 matched = pd.concat([treated, control.iloc[nn.kneighbors(treated[["age","sex","bmi"]])[1][:,0]]]) print("Matched sample size:", len(matched))

📋 PSM案例

数据：某新药的观察性研究，500例患者
结果：匹配前两组在年龄、性别上差异显著(P<0.05)；匹配后所有协变量标准化均差<0.1
解读：匹配后的OR=1.65(95%CI:1.12-2.43, P=0.012)，新药组结局显著改善。

IPTW逆概率加权

用途：通过加权法调整纺成，保留所有样本。

IPTW（Inverse Probability of Treatment Weighting）不会丢失样本，而是对每个个体赋予一个权重：W = T/PS + (1-T)/(1-PS)。加权后的幻想人群中处理分配与协变量独立。常用积极稳定化IPTW：SW = T*P(PS的均值)/PS + (1-T)*(1-P)/(1-PS)。

library(survey) set.seed(123); n <- 500 data <- data.frame( treatment = rbinom(n,1,0.3), age = rnorm(n,60,10), sex = rbinom(n,1,0.5), bmi = rnorm(n,25,3) ) data$outcome <- rbinom(n,1,plogis(-2+0.5*data$treatment+0.03*data$age)) # 计算倾向性评分 ps_model <- glm(treatment ~ age + sex + bmi, data=data, family=binomial()) data$ps <- predict(ps_model, type="response") # IPTW权重（稳定化） p_t <- mean(data$treatment) data$sw <- ifelse(data$treatment==1, p_t/data$ps, (1-p_t)/(1-data$ps)) # 加权回归 design <- svydesign(ids=~1, weights=~sw, data=data) svyglm(outcome ~ treatment, design=design, family=binomial()) # 检查均衡性：加权前后比较 summary(data$sw) # 检查极端权重 data$sw <- pmin(data$sw, 10) # 截断极端值

Python

import numpy as np, pandas as pd import statsmodels.api as sm from sklearn.linear_model import LogisticRegression np.random.seed(123); n = 500 data = pd.DataFrame({ "treatment": np.random.binomial(1, 0.3, n), "age": np.random.normal(60, 10, n), "sex": np.random.binomial(1, 0.5, n), "bmi": np.random.normal(25, 3, n) }) data["outcome"] = np.random.binomial(1, 1/(1+np.exp(-(-2 + 0.5*data["treatment"] + 0.03*data["age"])))) # 倾向性评分 ps = LogisticRegression().fit( data[["age","sex","bmi"]], data["treatment"] ).predict_proba(data[["age","sex","bmi"]])[:,1] # IPTW权重 p_t = data["treatment"].mean() data["sw"] = np.where(data["treatment"]==1, p_t/ps, (1-p_t)/(1-ps)) data["sw"] = np.minimum(data["sw"], 10) # 截断极端值 # 加权回归 X = sm.add_constant(data["treatment"]) wls = sm.WLS(data["outcome"], X, weights=data["sw"]).fit() print(wls.summary()) print(f"Treatment effect (IPTW): {wls.params['treatment']:.4f}")

📋 IPTW案例

数据：同PSM案例的500例患者
结果：加权后协变量均衡性显著改善
解读：IPTW保留了全部样本，统计检验力更高，适合样本量较小的情景。

DID双重差分

用途：政策评估，比较处理组和对照组在干预前后的变化差异。

DID（Difference-in-Differences）比较处理组和对照组在干预前后结局变化的差异。交互项系数即为因果效应。核心假设：平行趋势假设——处理组和对照组在干预前结局的变化趋势一致。

# DID模型：交互项 set.seed(123) n <- 200; time <- 0:1 data <- expand.grid(id=1:n, post=time) data$treatment <- rep(c(rep(1,n/2), rep(0,n/2)), each=2) data$age <- rnorm(n*2, 60, 10) # 平行趋势：处理前无差异，处理后不同 data$outcome <- with(data, 50 + 10*treatment*post + 0.3*age + rnorm(n*2, 0, 5)) # DID估计 did_model <- lm(outcome ~ treatment*post + age, data=data) summary(did_model) # 交互项系数 = DID估计值 cat("DID estimate:", coef(did_model)["treatment:post"], " ") # 平行趋势检验（处理前） pre_data <- subset(data, post==0) t.test(outcome ~ treatment, data=pre_data) # 应不显著

Python

import numpy as np, pandas as pd import statsmodels.api as sm np.random.seed(123) n = 200 data = pd.DataFrame({ "id": np.arange(1, n+1).repeat(2), "post": np.tile([0, 1], n), "treatment": np.repeat([1]*(n//2) + [0]*(n//2), 2), "age": np.random.normal(60, 10, n*2) }) data["outcome"] = (50 + 10*data["treatment"]*data["post"] + 0.3*data["age"] + np.random.normal(0, 5, n*2)) # DID回归 data["did"] = data["treatment"] * data["post"] X = sm.add_constant(data[["treatment", "post", "did", "age"]]) model = sm.OLS(data["outcome"], X).fit() print(model.summary()) print(f"DID estimate: {model.params['did']:.4f}") # 平行趋势检验 pre = data[data["post"]==0] from scipy.stats import ttest_ind print("Parallel trend p-value:", ttest_ind(pre[pre["treatment"]==1]["outcome"], pre[pre["treatment"]==0]["outcome"]).pvalue)

📋 DID案例

场景：某省2022年实施医保支付改革，评估对患者自付费用的影响
结果：DID估计效应=-12.5(95%CI:-18.3~-6.7, P<0.001)
解读：改革后患者自付费用显著下降，平行趋势假设检验通过，证据可靠。

反事实因果模型

理论：Rubin因果模型（潜在结果框架），个体处理效应(Y1-Y0)。

反事实框架由Donald Rubin提出。每个个体有两个潜在结局：受处理(Y1)和未受处理(Y0)，但只能观测其一。个体处理效应ATE = E[Y1-Y0]。解决方案：匹配、分层、加权、工具变量。核心假设：无干扰因素、存在正确的偏向数据。

library(Matching) set.seed(123); n <- 500 X <- rnorm(n); T <- rbinom(n,1,plogis(X)) Y1 <- 2 + 0.5*T + 0.3*X + rnorm(n) Y0 <- 2 + 0.3*X + rnorm(n) Y <- ifelse(T==1,Y1,Y0) # PSM估计ATE m <- Match(Y=Y, Tr=T, X=X, caliper=0.2) summary(m) # 线性回归估计 summary(lm(Y ~ T + X)) cat("ATE:", mean(Y1-Y0))

Python

import numpy as np, pandas as pd import statsmodels.api as sm np.random.seed(123); n = 500 X = np.random.normal(0,1,n) T = np.random.binomial(1, 1/(1+np.exp(-X))) Y1 = 2 + 0.5*T + 0.3*X + np.random.normal(0,1,n) Y0 = 2 + 0.3*X + np.random.normal(0,1,n) Y = np.where(T==1,Y1,Y0) print('True ATE:', np.mean(Y1-Y0)) # 回归估计 X2 = sm.add_constant(pd.DataFrame({'T':T,'X':X})) print(sm.OLS(Y, X2).fit().summary())

📋 反事实案例

场景：新药临床试验，观察性数据不可避免混杂
结果：ATE=0.48(95%CI:0.31-0.65)，推荐新药
解读：纺成因素控制后，新药效果显著（反事实框架的核心价值在于明确了个体层面的因果推断逻辑）。

中断时间序列

用途：评估干预措施在时间序列上的突变效应（分段回归）。

ITS（Interrupted Time Series）是评价干预措施影响的准实验设计。通过分段线性回归分析比较干预前后的截距变化（瞬时效应）和斜率变化（渐变效应）。需要假设干预前的趋势线性、干预时点明确。常用于政策评估、公共卫生干预效果评价。

# ITS: 分段回归 set.seed(123); n <- 60 time <- 1:n; level <- c(rep(0,30),rep(1,30)) trend_pre <- 0.2*1:30 trend_post <- 10 + 0.1*1:30 y <- c(trend_pre, trend_post) + rnorm(n,0,3) model <- lm(y ~ time + level + I(time*level)) summary(model) # 绘制ITSA图 library(ggplot2) df <- data.frame(time, y, level) ggplot(df, aes(time,y,color=factor(level))) + geom_line() + geom_smooth(method='lm',se=F)

Python

import numpy as np, pandas as pd import statsmodels.api as sm np.random.seed(123); n = 60 time = np.arange(1, n+1) level = np.concatenate([np.zeros(30), np.ones(30)]) trend_pre = 0.2*np.arange(1,31) trend_post = 10 + 0.1*np.arange(1,31) y = np.concatenate([trend_pre, trend_post]) + np.random.normal(0,3,n) df = pd.DataFrame({'time':time,'level':level,'y':y}) df['trend'] = df['time'] * df['level'] X = sm.add_constant(df[['time','level','trend']]) model = sm.OLS(df['y'], X).fit() print(model.summary())

📋 ITS案例

场景：某省2021年实施吸烟禁令后心脏病入院率变化
结果：瞬时效应=-8.5(P<0.001)，渐变效应=-0.3/月(P=0.042)
解读：禁令实施后心脏病入院显著下降，且下降趋势持续。

工具变量

用途：当存在未观测纺成因素时，用工具变量Z推断X对Y的因果效应。

工具变量(IV)法解决未观测纺成问题。有效的IV需满足：(1)相关性——Z与X强相关；(2)外生性——Z仅通过X影响Y。常用方法：2SLS二阶段最小二乘、工具变量选择：化疗距离作为治疗方案的IV。

library(AER) set.seed(123); n <- 500 Z <- rnorm(n) # IV: 医院距离 X <- 0.5*Z + rnorm(n) # 治疗方案 Y <- 2 + 0.8*X + rnorm(n) + 0.3*rnorm(n) # 普通OLS估计(偏误) summary(lm(Y ~ X)) # 2SLS工具变量 iv <- ivreg(Y ~ X | Z, data=data.frame(Y,X,Z)) summary(iv)

Python

import numpy as np, pandas as pd from linearmodels.iv import IV2SLS np.random.seed(123); n = 500 Z = np.random.normal(0,1,n) # IV X = 0.5*Z + np.random.normal(0,1,n) Y = 2 + 0.8*X + np.random.normal(0,1,n) df = pd.DataFrame({'Y':Y,'X':X,'Z':Z}) # 2SLS iv = IV2SLS.from_formula('Y ~ 1 + [X ~ Z]', df).fit() print(iv) print('True effect: 0.8')

📋 IV案例

场景：研究放疗治疗智力降低是否提高生活质量
IV：住院距离质子治疗中心的距离
结果：2SLS估计效应=0.75(P=0.003)，高于OLS的0.52
解读：OLS低估效应，IV纠正了未观测纺成的偏误。

断点回归

用途：根据分配变量的阈值判断处理效应。

RDD（Regression Discontinuity）利用分配变量的阈值划分处理组和对照组。在阈值附近的个体可视为随机分配。核心假设：纺成因素在阈值处连续。Sharp RDD（阈值严格分配）和Fuzzy RDD（阈值概率分配）。

library(rdd) set.seed(123); n <- 1000 X <- runif(n, 0, 100) # 考试成绩 cutoff <- 60 # 获奖阈值 T <- ifelse(X >= cutoff, 1, 0) Y <- 50 + 10*T + 0.5*X + rnorm(n, 0, 15) # 反事实分析 rd <- RDestimate(Y ~ X, cutpoint=cutoff) summary(rd) # 自带图 plot(rd)

Python

import numpy as np, pandas as pd import statsmodels.api as sm np.random.seed(123); n = 1000 X = np.random.uniform(0,100,n) cutoff = 60 T = (X >= cutoff).astype(int) Y = 50 + 10*T + 0.5*X + np.random.normal(0,15,n) df = pd.DataFrame({'X':X,'T':T,'Y':Y}) # 局部线性回归(阈值附近50%窗口) df['dist'] = df['X'] - cutoff df['interact'] = df['T'] * df['dist'] window = df[df['dist'].abs() < 10] X2 = sm.add_constant(window[['T','dist','interact']]) model = sm.OLS(window['Y'], X2).fit() print('RDD effect:', model.params['T'])

📋 RDD案例

场景：考试成绩≥60分获得奖学金，研究奖学金对后续学业的影响
结果：阈值处效应=9.8(95%CI:4.2-15.4), P=0.001
解读：获奖学生后续学业成绩显著提升，证据支持奖学金的激励效应。

有向无环图

用途：DAGs可视化因果假设，识别纺成、中介、操控策略。

DAG（Directed Acyclic Graph）是因果推断的图论工具。符号规则：X→Y表示X是Y的原因。基础概念：纺成因子（Confounder）——同时影响X和Y，需控制；中介变量（Mediator）——位于X→Y路径上，不应控制；操控变量（Collider）——被多个变量共同影响，控制会引入偏误。

library(dagitty) # 构建DAG dag <- dagitty('dag { X [exposure,pos="0,0"] Y [outcome,pos="2,0"] Z [pos="1,-1"] C [pos="1,1"] X -> Y C -> X; C -> Y # 混杂 Z -> X; Z -> Y # 共同原因 }') plot(dag) # 最小充分控制集 adjustmentSets(dag, exposure="X", outcome="Y") # 因果效应识别 cat("Control sets:") print(adjustmentSets(dag, "X", "Y"))

Python

import networkx as nx, matplotlib.pyplot as plt # 构建DAG G = nx.DiGraph() G.add_edges_from([('X','Y'),('C','X'),('C','Y'),('Z','X'),('Z','Y')]) pos = {'X':(0,0),'Y':(2,0),'C':(1,-1),'Z':(1,1)} # 查找混杂路径 print("Nodes:", list(G.nodes())) print("Edges:", list(G.edges())) nx.draw(G, pos, with_labels=True, node_color='lightblue', node_size=2000, arrowsize=20) plt.title('DAG示例') plt.show() # 混杂因子: C和Z均需控制

📋 DAG案例

场景：研究运动健康效应（X→Y）
DAG分析：年龄(C)、性别(Z)都是纺成因子，必须控制
结果：控制C和Z后，运动对健康的因果效应仍显著
解读：DAG帮助选择正确的控制变量集，避免误控制或漏控制。

🔬 高级统计模型

🤖 scikit-learn 官方文档

🤖 医学机器学习

LASSO回归

L1正则化，自动特征选择，适合高维数据

library(glmnet) set.seed(123) x <- matrix(rnorm(100*20), 100, 20) y <- rnorm(100) cv_lasso <- cv.glmnet(x, y, alpha=1) best_lambda <- cv_lasso$lambda.min lasso_model <- glmnet(x, y, alpha=1, lambda=best_lambda) coef(lasso_model) # 交叉验证图 plot(cv_lasso)

Python

from sklearn.linear_model import LassoCV from sklearn.preprocessing import StandardScaler import numpy as np X = np.random.randn(100, 20) y = np.random.randn(100) scaler = StandardScaler() X_scaled = scaler.fit_transform(X) lasso = LassoCV(cv=5, random_state=123).fit(X_scaled, y) print(f'最佳alpha: {lasso.alpha_:.3f}') print(f'保留的特征数: {np.sum(lasso.coef_!=0)}')

随机森林

集成学习，变量重要性排序，适合复杂交互

library(randomForest) set.seed(123) rf_model <- randomForest(Species ~ ., data=iris, ntree=500, importance=TRUE) print(rf_model) # 变量重要性 importance(rf_model) varImpPlot(rf_model) # 调参 tuneRF(iris[,-5], iris[,5], ntreeTry=500)

Python

from sklearn.ensemble import RandomForestClassifier from sklearn.inspection import permutation_importance rf = RandomForestClassifier(n_estimators=500, random_state=123) rf.fit(X_train, y_train) print(f'准确率: {rf.score(X_test, y_test):.3f}') # 特征重要性 importances = rf.feature_importances_ # 排列重要性 perm_imp = permutation_importance(rf, X_test, y_test)

XGBoost

梯度提升，Kaggle常用，支持正则化

library(xgboost) dtrain <- xgb.DMatrix(data=as.matrix(iris[,-5]), label=as.numeric(iris$Species)-1) params <- list(objective="multi:softprob", num_class=3, max_depth=3, eta=0.1) xgb_model <- xgb.train(params=params, data=dtrain, nrounds=100, verbose=0) importance <- xgb.importance(feature_names=colnames(iris[,-5]), model=xgb_model) xgb.plot.importance(importance)

Python

import xgboost as xgb from sklearn.datasets import load_iris from sklearn.model_selection import train_test_split data = load_iris() X_train, X_test, y_train, y_test = train_test_split( data.data, data.target, test_size=0.2, random_state=123) model = xgb.XGBClassifier(n_estimators=100, max_depth=3, learning_rate=0.1, random_state=123) model.fit(X_train, y_train) print(f'AUC: {model.score(X_test, y_test):.3f}') xgb.plot_importance(model)

SVM支持向量机

小样本、非线性分类，核函数技巧

library(e1071) svm_model <- svm(Species ~ ., data=iris, kernel="radial", cost=1) summary(svm_model) # 调参 tune.svm(Species ~ ., data=iris, gamma=10^(-3:1), cost=10^(-1:2))

Python

from sklearn.svm import SVC from sklearn.model_selection import GridSearchCV svm = SVC(kernel='rbf', probability=True, random_state=123) param_grid = {'C': [0.1, 1, 10], 'gamma': [0.01, 0.1, 1]} grid = GridSearchCV(svm, param_grid, cv=5) grid.fit(X_train, y_train) print(f'最佳参数: {grid.best_params_}') print(f'最佳得分: {grid.best_score_:.3f}')

神经网络

多层感知器(MLP)用于医学数据分类/回归

library(nnet) nn <- nnet(Species ~ ., data=iris, size=5, decay=0.01, maxit=200, trace=FALSE) pred <- predict(nn, iris, type="class") table(pred, iris$Species)

Python

from sklearn.neural_network import MLPClassifier from sklearn.preprocessing import StandardScaler scaler = StandardScaler() X_scaled = scaler.fit_transform(X_train) mlp = MLPClassifier(hidden_layer_sizes=(50,25), activation='relu', max_iter=500, random_state=123) mlp.fit(X_scaled, y_train) print(f'准确率: {mlp.score(scaler.transform(X_test), y_test):.3f}')

SHAP可解释性

解释模型预测，全球/局部可解释

library(shapr) library(xgboost) model <- xgboost(data=as.matrix(iris[,-5]), label=as.numeric(iris$Species)-1, nrounds=20, verbose=0) explainer <- shapr(iris[,-5], model) p <- explain(iris[1:5,-5], explainer, approach="empirical") plot(p)

Python

import shap import xgboost as xgb model = xgb.XGBClassifier().fit(X_train, y_train) explainer = shap.TreeExplainer(model) shap_values = explainer.shap_values(X_test) # 全局特征重要性 shap.summary_plot(shap_values, X_test) # 单个样本解释 shap.force_plot(explainer.expected_value, shap_values[0], X_test[0])

🔮 Stan 贝叶斯建模

🔮 贝叶斯统计

贝叶斯基础

核心公式：P(θ|D) ∝ P(D|θ)×P(θ)
先验：基于已有知识设定的参数分布
似然：数据的分布假设
后验：结合先验和数据后的参数分布
MCMC：Gibbs采样、Metropolis-Hastings

贝叶斯t检验

提供组间差异的后验分布

library(BayesFactor) set.seed(123) group1 <- rnorm(30, 100, 15) group2 <- rnorm(30, 110, 15) # 贝叶斯t检验 bf <- ttestBF(group1, group2) print(bf) # 提取后验样本 chains <- posterior(bf, iterations=10000) summary(chains) # 绘制后验分布 plot(chains[, "mu"])

贝叶斯回归

系数概率分布，直接概率解读

library(brms) set.seed(123) data <- data.frame( age = rnorm(100, 60, 10), bmi = rnorm(100, 25, 3), sbp = rnorm(100, 130, 15) ) # 贝叶斯线性回归 b_model <- brm(sbp ~ age + bmi, data=data, prior=c(set_prior("normal(0,10)", class="b")), chains=4, iter=2000) summary(b_model) # 系数后验概率 hypothesis(b_model, "age > 0") # 预测 predict(b_model, newdata=data[1:5,])

📋 Cochrane 循证医学

📋 Meta分析

Meta分析基础

模型：固定效应、随机效应
效应量：OR/RR/HR、均数差(SMD)
异质性：I²统计量、Q检验
偏倚：发表偏倚(funnel plot)、Begg检验、Egger检验
软件：R(meta/metafor)、Stata、RevMan

Meta分析R实现

完整分析流程

library(meta) # 模拟数据：5项研究 study <- c("研究1","研究2","研究3","研究4","研究5") events_t <- c(45, 30, 62, 28, 50) total_t <- c(100, 80, 120, 60, 95) events_c <- c(30, 25, 40, 22, 35) total_c <- c(100, 82, 115, 58, 90) # 固定效应模型 (MH法) m_fixed <- metabin(events_t, total_t, events_c, total_c, studlab=study, sm="OR", method="MH", fixed=TRUE, random=FALSE) summary(m_fixed) # 随机效应模型 (DL法) m_random <- metabin(events_t, total_t, events_c, total_c, studlab=study, sm="OR", method="MH", fixed=FALSE, random=TRUE) summary(m_random) # 森林图 forest(m_random) # 漏斗图 + Egger检验 funnel(m_random) metabias(m_random, method.bias="linreg") # 亚组分析 # update(m_random, by=subgroup) # 敏感性分析 metainf(m_random) # 逐一剔除

📋 结果解读

合并OR = 1.35 (95%CI: 1.08-1.69)
I² = 32.4% (p=0.21) → 异质性低，选固定效应
Egger检验 p=0.45 → 无发表偏倚
森林图显示各研究效应方向一致
结论：治疗组效果显著优于对照组

📈 WHO 建模指南

📈 传染病数学建模

传染病建模是预测疾病传播、评估防控措施的核心方法。

一、核心概念

S(t)：易感者 | E(t)：潜伏者 | I(t)：感染者 | R(t)：康复者 | M(t)：死亡者

β：传染率 | σ：潜伏转染率 | γ：康复率 | ξ：免疫丧失率 | μ：死亡率

dX/dt = 人群每日变化量 | R₀=β/γ：基本再生数

二、四大经典模型

1. SIR

dS/dt=-βSI/N dI/dt=βSI/N-γI dR/dt=γI

2. SEIR

dS/dt=-βSI/N dE/dt=βSI/N-σE dI/dt=σE-γI dR/dt=γI

3. SEIRS

dS/dt=-βSI/N+ξR dE/dt=βSI/N-σE dI/dt=σE-γI dR/dt=γI-ξR

4. MSEIR

dS/dt=-βSI/N dE/dt=βSI/N-σE dI/dt=σE-(γ+μ)I dR/dt=γI dM/dt=μI

三、参数表

符号	名称	计算	参考
β	传染率	接触率×感染概率	CDC
σ	潜伏率	1/潜伏期	临床数据
γ	康复率	1/传染期	疾病自然史
R₀	再生数	β/γ	WHO

四、Python实战

Python

import numpy as np from scipy.integrate import odeint import matplotlib.pyplot as plt N=10000; E0,I0,R0=0,1,0; S0=N-E0-I0-R0 beta,sigma,gamma=0.3,1/5,1/10 t=np.linspace(0,150,150) def seir(y,t,N,beta,sigma,gamma): S,E,I,R=y dS=-beta*S*I/N dE=beta*S*I/N-sigma*E dI=sigma*E-gamma*I dR=gamma*I return dS,dE,dI,dR sol=odeint(seir,[S0,E0,I0,R0],t,args=(N,beta,sigma,gamma)) S,E,I,R=sol.T plt.plot(t,S,label='易感者') plt.plot(t,E,label='潜伏者') plt.plot(t,I,label='感染者',c='red') plt.plot(t,R,label='康复者') plt.xlabel('天数');plt.ylabel('人数') plt.title('SEIR模型模拟') plt.legend();plt.show()

SIR模型实践

用途：经典传染病分室模型，模拟从感染到康复的全过程。

SIR模型将人群分为三个分室：S(易感者)、I(感染者)、R(康复者)。微分方程组：dS/dt = -βSI/N，dI/dt = βSI/N – γI，dR/dt = γI。基本再生数R₀ = β/γ，表示一个感染者在完全易感人群中平均传播的人数。

library(deSolve) # SIR模型参数 N <- 10000; I0 <- 1; R0 <- 0; S0 <- N - I0 - R0 beta <- 0.4; gamma <- 0.1 parameters <- c(beta=beta, gamma=gamma) state <- c(S=S0, I=I0, R=R0) times <- seq(0, 100, by=1) # SIR微分方程 sir <- function(t, state, parameters) { with(as.list(c(state, parameters)), { dS <- -beta * S * I / N dI <- beta * S * I / N - gamma * I dR <- gamma * I return(list(c(dS, dI, dR))) }) } output <- ode(y=state, times=times, func=sir, parms=parameters) plot(output[,1], output[,3], type='l', col='red', lwd=2, xlab='天数', ylab='人数', main='SIR模型模拟') lines(output[,1], output[,2], col='blue', lwd=2) lines(output[,1], output[,4], col='green', lwd=2) legend('topright', c('易感者','感染者','康复者'), col=c('blue','red','green'), lty=1, lwd=2) cat('R0 =', beta/gamma, ' ') cat('最终感染比例:', max(output[,3])/N*100, '% ')

Python

import numpy as np from scipy.integrate import odeint import matplotlib.pyplot as plt N = 10000; beta, gamma = 0.4, 0.1 def sir(y, t): S, I, R = y dS = -beta * S * I / N dI = beta * S * I / N - gamma * I dR = gamma * I return dS, dI, dR sol = odeint(sir, [N-1, 1, 0], np.linspace(0, 100, 101)) S, I, R = sol.T plt.figure(figsize=(8,5)) plt.plot(S, label="易感者", color="blue") plt.plot(I, label="感染者", color="red") plt.plot(R, label="康复者", color="green") plt.xlabel("天数"); plt.ylabel("人数") plt.title(f"SIR模型 (R0={beta/gamma:.1f})") plt.legend(); plt.grid(alpha=0.3) plt.show() print(f"R0 = {beta/gamma:.1f}") print(f"峰值感染: {max(I):.0f}人 ({max(I)/N*100:.1f}%)")

📋 SIR案例

场景：流感爆发模拟，人口N=10000，R₀=4.0
结果：感染高峰出现在级33天，最大日新增约为总人口的8%，累计感染比例约98%
解读：R₀>1时疫情将大规模爆发，必须采取干预措施降低R₀。

SEIR模型实践

用途：增加潜伏期分室，更精确地模拟具有潜伏期的传染病。

SEIR模型在SIR的基础上增加E(潜伏者)分室。微分方程组：dS/dt = -βSI/N，dE/dt = βSI/N – σE，dI/dt = σE – γI，dR/dt = γI。σ表示潜伏期转染率(1/潜伏期)。适用于新冠、麻疹、结核等具有潜伏期的传染病。

library(deSolve) N <- 10000; E0 <- 0; I0 <- 1; R0 <- 0; S0 <- N - E0 - I0 - R0 beta <- 0.3; sigma <- 1/5; gamma <- 1/10 parameters <- c(beta=beta, sigma=sigma, gamma=gamma) state <- c(S=S0, E=E0, I=I0, R=R0) times <- seq(0, 150, by=1) seir <- function(t, state, parameters) { with(as.list(c(state, parameters)), { dS <- -beta * S * I / N dE <- beta * S * I / N - sigma * E dI <- sigma * E - gamma * I dR <- gamma * I return(list(c(dS, dE, dI, dR))) }) } output <- ode(y=state, times=times, func=seir, parms=parameters) matplot(output[,1], output[,2:5], type='l', lwd=2, lty=1, xlab='天数', ylab='人数', main='SEIR模型模拟') legend('right', c('S','E','I','R'), col=1:4, lty=1, lwd=2) cat('R0 =', beta/gamma, ' ') cat('潜伏期 =', 1/sigma, '天 ')

Python

import numpy as np from scipy.integrate import odeint import matplotlib.pyplot as plt N = 10000; beta, sigma, gamma = 0.3, 1/5, 1/10 def seir(y, t): S, E, I, R = y dS = -beta * S * I / N dE = beta * S * I / N - sigma * E dI = sigma * E - gamma * I dR = gamma * I return dS, dE, dI, dR sol = odeint(seir, [N-1, 0, 1, 0], np.linspace(0, 150, 151)) S, E, I, R = sol.T plt.figure(figsize=(8,5)) plt.plot(S, label="易感者"); plt.plot(E, label="潜伏者") plt.plot(I, label="感染者", color="red") plt.plot(R, label="康复者") plt.xlabel("天数"); plt.ylabel("人数") plt.title(f"SEIR模型 (R0={beta/gamma:.1f}, 潜伏期={1/sigma:.0f}天)") plt.legend(); plt.grid(alpha=0.3); plt.show() print(f"R0 = {beta/gamma:.1f}") print(f"潜伏期 = {1/sigma:.0f}天")

📋 SEIR案例

场景：COVID-19传播模拟，潜伏期5天，R₀=3.0
结果：感染高峰在约70天，比SIR模型延迟约一倍
解读：潜伏期的存在显著延缓了疫情高峰，但总体爆发规模仍然巨大，提醒早期干预至关重要。

随机SEIR模型

用途：考虑随机性的传染病模型，适合小规模人群爆发模拟。

随机SEIR模型使用二项分布/泊松分布替代微分方程，反映传染过程的随机性。在小人群中，随机性可能导致疫情自行熔灭或爆发。每步计算新感染者~Binom(S, βI/N)，新潜伏转感染者~Binom(E, σ)，新康复者~Binom(I, γ)。

# 随机SEIR模型 (离散时间) set.seed(42) N <- 500; S <- N-1; E <- 0; I <- 1; R <- 0 beta <- 0.5; sigma <- 0.2; gamma <- 0.1 days <- 100 results <- matrix(0, nrow=days, ncol=4) colnames(results) <- c('S','E','I','R') for (t in 1:days) { new_exposed <- rbinom(1, S, min(beta*I/N, 1)) new_infectious <- rbinom(1, E, sigma) new_recovered <- rbinom(1, I, gamma) S <- S - new_exposed E <- E + new_exposed - new_infectious I <- I + new_infectious - new_recovered R <- R + new_recovered results[t,] <- c(S, E, I, R) if (I == 0 && E == 0) break } matplot(1:t, results[1:t,], type='l', lwd=2, lty=1, xlab='天数', ylab='人数', main='随机SEIR模拟') legend('right', c('S','E','I','R'), col=1:4, lty=1, lwd=2) cat('最终感染:', results[1,4], '/', N, ' ') # 多次模拟显示随机性 sims <- 100; final_size <- numeric(sims) for (s in 1:sims) { # 简写 S <- N-1; E <- 0; I <- 1; R <- 0 for (t in 1:days) { new_exposed <- rbinom(1, S, min(beta*I/N, 1)) new_infectious <- rbinom(1, E, sigma) new_recovered <- rbinom(1, I, gamma) S <- S - new_exposed; E <- E + new_exposed - new_infectious I <- I + new_infectious - new_recovered; R <- R + new_recovered if (I == 0 && E == 0) break } final_size[s] <- R } hist(final_size, main='100次模拟最终感染分布', xlab='最终康复人数')

Python

import numpy as np import matplotlib.pyplot as plt np.random.seed(42) N, beta, sigma, gamma = 500, 0.5, 0.2, 0.1 days = 100 def run_one_sim(): S, E, I, R = N-1, 0, 1, 0 hist = [] for _ in range(days): new_e = np.random.binomial(S, min(beta*I/N, 1)) new_i = np.random.binomial(E, sigma) new_r = np.random.binomial(I, gamma) S -= new_e; E += new_e - new_i I += new_i - new_r; R += new_r hist.append([S, E, I, R]) if I == 0 and E == 0: break return np.array(hist) # 单次模拟 result = run_one_sim() plt.figure(figsize=(10,4)) plt.subplot(121) plt.plot(result[:,0], label="S"); plt.plot(result[:,1], label="E") plt.plot(result[:,2], label="I", color="red") plt.plot(result[:,3], label="R") plt.legend(); plt.title("单次随机SEIR模拟") # 多次模拟 sims = 100 final_sizes = [run_one_sim()[-1,3] for _ in range(sims)] plt.subplot(122) plt.hist(final_sizes, bins=15) plt.title(f"100次模拟最终感染分布"); plt.xlabel("康复人数") plt.tight_layout(); plt.show()

📋 随机SEIR案例

场景：小城镇500人的传染病爆发模拟
结果：约15%的模拟中疫情自行熔灭，85%爆发
解读：随机模型揭示了小人群中的重要现象——即便R₀>1，随机性也可使疫情未必爆发。

Branching传播链模型

用途：基于再生数的传播链模型，评估爆发概率和传播链长度。

Branching过程将传染病传播视为一个传播链，每个感染者产生的二代病例数服从某分布(如泊松分布、负二项分布)。核心概念：时空再生数R（复制数），R<1时传播链必然熔灭，R>1时有正概率爆发。

# Branching过程模拟 set.seed(123) # 负二项分布(过度离散) sim_branch <- function(R=2.5, k=0.5, gen_max=20) { chain <- 1 # 初始病例 total <- 1 for (g in 1:gen_max) { offspring <- sum(rnbinom(chain, mu=R, size=k)) total <- total + offspring chain <- offspring if (chain == 0) break } return(list(total=total, gens=g, extinct=chain==0)) } # 多次模拟计算暴发概率 sims <- 1000 R_vals <- c(0.8, 1.0, 1.5, 2.0, 3.0) for (R in R_vals) { extinct <- 0 sizes <- numeric(sims) for (s in 1:sims) { res <- sim_branch(R=R, k=0.5) if (res$extinct) extinct <- extinct + 1 sizes[s] <- res$total } cat(sprintf('R=%.1f: 暴发概率=%.1f%%, 平均规模=%.0f\n', R, (1-extinct/sims)*100, mean(sizes))) } # 可视化 R0 <- 2.5 res <- replicate(1000, sim_branch(R=R0, k=0.5)$total) hist(log10(res+1), main='Branching过程最终规模分布', xlab='log10(总病例数)', col='steelblue') abline(v=log10(100), col='red', lwd=2, lty=2) cat('R0=', R0, '暴发(>100例)概率:', mean(res>100)*100, '% ')

Python

import numpy as np import matplotlib.pyplot as plt np.random.seed(123) def sim_branch(R=2.5, k=0.5, gen_max=20): chain = 1; total = 1 for g in range(gen_max): # 负二项 offspring (mu=R, size=k 控制离散度) offspring = np.sum(np.random.negative_binomial(k, k/(k+R), chain)) total += offspring chain = offspring if chain == 0: break return total, g+1, chain == 0 # 不同R0下的暴发概率 sims = 1000 for R in [0.8, 1.0, 1.5, 2.0, 3.0]: extinct = 0; sizes = [] for s in range(sims): tot, _, ext = sim_branch(R=R) sizes.append(tot) if ext: extinct += 1 print(f"R={R:.1f}: 暴发概率={(1-extinct/sims)*100:.1f}%, " f"平均规模={np.mean(sizes):.0f}") # R0=2.5时规模分布 res = [sim_branch(R=2.5)[0] for _ in range(1000)] plt.hist(np.log10(np.array(res)+1), bins=30) plt.axvline(np.log10(100), color='red', ls='--') plt.title("Branching过程最终规模分布"); plt.xlabel("log10(总病例数)") plt.show() print(f"暴发(>100例)概率: {np.mean(np.array(res)>100)*100:.1f}%")

📋 Branching案例

场景：超级传播者事件模拟——一个患者可以传播多少人
结果：R₀=2.5时爆发概率≈85%，但有显著的二态性（熔灭或大规模爆发）
解读：Branching模型揭示了传播链的随机性，解释为什么同样的R₀可以导致完全不同的爆发结局。

📋 研究报告规范与检查清单

医学研究报告规范是高质量科研的基石。以下是各类型研究应遵循的报告规范，投稿前请逐项核对。

STROBE

观察性研究报告规范
横断面/队列/病例对照
查看详情 →

CONSORT

随机对照试验报告规范
含流程图模板
查看详情 →

PRISMA

系统评价/Meta分析报告规范
27项核查清单
查看详情 →

STARD

诊断准确性试验报告规范
25项核查清单
查看详情 →

EQUATOR

报告规范总库
500+指南汇总
查看详情 →

RECORD

真实世界数据研究报告规范
常规收集数据
查看详情 →

COSMIN

量表/测量工具报告规范
信效度评价
查看详情 →

ARRIVE

动物实验报告规范
10项核心清单
查看详情 →

TRIPOD

预测模型研究报告规范
个体化预测
查看详情 →

SPIRIT

临床试验方案报告规范
含SPIRIT图
查看详情 →

🏥 临床实践与规范

📝 临床研究设计

证据等级：RCT > 队列 > 病例对照 > 横断面
观察性：横断面、病例对照(OR)、队列(RR，遵循STROBE)
实验性：RCT(遵循CONSORT)、诊断试验(遵循STARD)
量表：遵循COSMIN

🏥 真实世界研究

数据来源：EHR、医保、登记库
混杂控制：PSM、IPTW、DID、工具变量
报告规范：RECORD

⏳ 生存分析

基础：Kaplan-Meier、Log-rank检验
多因素：Cox比例风险回归
前提：比例风险假设PH检验
拓展：竞争风险、时变协变量

📰 主要医学期刊统计要求

国际顶级医学期刊对统计方法报告有严格要求，投稿前请对照核查。

NEJM

统计要求：
• 详细描述所有统计方法，引用文献
• 报告效应量(OR/RR/HR)及95%CI
• P值精确报告（如P=0.03而非P<0.05）
• ITT原则分析RCT数据
• 多重比较需校正（Bonferroni等）
• 亚组分析需报告交互作用检验

JAMA

统计要求：
• 遵循CONSORT/STROBE等报告规范
• 提供统计分析计划(SAP)
• 效应量+95%CI优先于P值
• 所有分析需注明软件版本
• 缺失数据需说明处理方法
• 倾向性评分需说明匹配方法

Lancet

统计要求：
• 统计学专家需列为共同作者或致谢
• 注册临床试验需提供注册号
• 完整报告所有结局（含阴性结果）
• 亚组分析需预先指定并控制假阳性
• 需报告数据监测委员会(DMC)信息

BMJ

统计要求：
• 公开原始数据和代码（鼓励）
• 使用STAT-CHECK统计核查流程
• 报告绝对效应和相对效应
• 需报告Number Needed to Treat(NNT)
• 贝叶斯分析需说明先验选择

Nature Medicine

统计要求：
• 详细统计方法需在Methods中描述
• 生物信息分析需提供代码仓库
• 多重假设检验需控制FDR
• 高维数据分析需说明过拟合控制
• 效应量报告需含变异度量

中文核心期刊

统计要求：
• 遵循国内统一统计学报告规范
• 需说明随机分组方法及隐藏方案
• 对照组需说明是否同期对照
• 需报告所有统计方法的引用来源
• 伦理审批需注明批件号

💡 通用建议

始终遵循EQUATOR网络的相关报告规范，在Methods部分详细描述所有统计方法（含软件版本），提供充分的效应量和精度估计，并请统计学专家审阅稿件。

⚠️ 医学科研常见统计错误与避坑指南

❌ P值滥用

错误：P<0.05=有差异，P>0.05=无差异
正确：P值仅为连续度量，需结合效应量和CI解读。P>0.05可能因样本量不足。避免"P值黑客"(P-hacking)。
建议：报告精确P值和95%CI，而非仅标P<0.05

❌ 多重比较不校正

错误：多次检验直接使用α=0.05
正确：多组比较用ANOVA（而非多组t检验），多重终点需Bonferroni/FDR校正。探索性分析结果需明确说明。
建议：预先指定主要和次要结局

❌ 忽略数据前提检验

错误：直接做参数检验不验证前提
正确：参数检验前必须检验正态性(Shapiro-Wilk)和方差齐性(Levene)。不满足时用非参数检验或数据变换。
建议：将前提检验结果附在补充材料

❌ 过度解读亚组分析

错误：亚组分析发现某组显著即下结论
正确：亚组结果需做交互作用检验(interaction test)，阳性交互作用才是真正差异。未预先指定的亚组分析为探索性。
建议：亚组结果仅作为假设生成

❌ 忽略混杂因素

错误：观察性研究直接比较组间差异
正确：需多因素回归控制混杂，或使用PSM、IPTW、DID等因果推断方法。需报告混杂因素选择依据（DAG图）。
建议：使用有向无环图(DAG)识别混杂

❌ 缺失数据处理不当

错误：直接删除(complete-case analysis)
正确：需说明缺失机制(MCAR/MAR/MNAR)，推荐多重插补(MICE)，进行敏感性分析比较不同处理方法的结果。
建议：在方法部分详细描述缺失处理策略

❌ 样本量不足

错误：事后用效力分析解释阴性结果
正确：研究设计阶段即应进行样本量计算，需明确假设的效应量来源和检验效能(通常80%或90%)。
建议：在研究方案中纳入样本量计算部分

❌ 选择性报告阳性结果

错误：只报告有统计学意义的结果
正确：需报告所有预先指定结局，无论阳性或阴性。阴性结果同样有发表价值，避免发表偏倚。
建议：预先注册研究方案(ClinicalTrials.gov/中国临床试验注册中心)

📚 科研资源与工具

🧬 Bioconductor 工具库

🧬 生物信息学分析方法

差异表达分析

筛选不同组间表达量显著差异的基因

library(edgeR) y <- DGEList(counts=count_matrix, group=sample_info$group) y <- calcNormFactors(y) design <- model.matrix(~0+group, data=sample_info) y <- estimateDisp(y, design) fit <- glmQLFit(y, design) qlf <- glmQLFTest(fit, contrast=c(1,-1)) deg <- topTags(qlf, n=Inf)

Python

import scanpy as sc adata = sc.read_h5ad('data.h5ad') sc.tl.rank_genes_groups(adata, 'group', method='wilcoxon') deg = sc.get.rank_genes_groups_df(adata, group='treatment')

工具：edgeR、DESeq2、limma

功能富集分析

GO/KEGG富集

library(clusterProfiler) go_enrich <- enrichGO(gene=deg_gene, OrgDb=org.Hs.eg.db, keyType="ENTREZID", ont="BP", pAdjustMethod="fdr", qvalueCutoff=0.05) dotplot(go_enrich, showCategory=15) kegg_enrich <- enrichKEGG(gene=deg_gene, organism='hsa', pvalueCutoff=0.05) dotplot(kegg_enrich)

Python

import gseapy as gp go_res = gp.enrichr(gene_list=deg_genes, gene_sets='GO_Biological_Process_2023', organism='Human', cutoff=0.05) kegg_res = gp.enrichr(gene_list=deg_genes, gene_sets='KEGG_2023_Human')

WGCNA

基因共表达网络

library(WGCNA) powers=c(c(1:10), seq(from=12, to=20, by=2)) sft=pickSoftThreshold(exprMatrix, powerVector=powers, verbose=5) net=blockwiseModules(exprMatrix, power=sft$powerEstimate, TOMType="unsigned", minModuleSize=30, mergeCutHeight=0.25, numericLabels=TRUE) moduleTraitCor=cor(net$MEs, trait, use="p")

免疫浸润分析

肿瘤微环境免疫细胞浸润

library(CIBERSORT) sig_matrix=read.table('LM22.txt', header=T, row.names=1, sep='\t') results=CIBERSORT(sig_matrix, exprMatrix, perm=1000, QN=TRUE) barplot(t(results[,1:22]), col=rainbow(22), border=NA)

算法：CIBERSORT、xCell、MCP-counter

单细胞测序分析

细胞分群、轨迹分析、细胞通讯

library(Seurat) sce <- CreateSeuratObject(counts=counts, project="scRNA", min.cells=3, min.features=200) sce <- NormalizeData(sce) sce <- FindVariableFeatures(sce) sce <- ScaleData(sce) sce <- RunPCA(sce) sce <- FindNeighbors(sce, dims=1:20) sce <- FindClusters(sce, resolution=0.5) sce <- RunUMAP(sce, dims=1:20) DimPlot(sce, reduction="umap", label=TRUE) # 差异表达 markers <- FindAllMarkers(sce, only.pos=TRUE)

工具：Seurat、Scanpy、Monocle、CellPhoneDB

💾 医学公共数据库

▼

🔬 科研工具

▼

🤖 AI科研工具

▼

📋 常用医学量表

▼

📋 COSMIN 量表规范

量表用于量化患者症状、功能、生活质量，需经严格信效度验证。

📊 量表信效度分析

信度：Cronbach's α、重测信度、分半信度
效度：内容效度、结构效度、效标效度
结构分析：EFA、CFA

📝 R语言实操案例

以下是完整的医学统计R语言分析案例，包含模拟数据、完整代码和结果解读。

🏥 案例1：两组临床指标比较

t检验R基础结果报告

比较新药组和安慰剂组患者的血压变化，包含描述性统计、t检验和结果报告。

# ===== 案例1：两组临床指标比较 ===== # 研究：新药A vs 安慰剂对收缩压的影响 # 作者：Abel医研统计 # 1. 模拟数据 set.seed(123) n <- 40 drug <- rnorm(n, mean=125, sd=8) placebo <- rnorm(n, mean=135, sd=10) # 2. 描述性统计 cat("=== 描述性统计 ===\n") cat("药物组：", round(mean(drug),1), "±", round(sd(drug),1), "mmHg\n") cat("安慰剂组：", round(mean(placebo),1), "±", round(sd(placebo),1), "mmHg\n") # 3. 正态性检验 shapiro.test(drug) shapiro.test(placebo) # 4. 方差齐性检验 var.test(drug, placebo) # 5. 独立样本t检验 t_result <- t.test(drug, placebo, var.equal=TRUE) print(t_result) # 6. 效应量 (Cohen's d) library(effsize) cohen.d(drug, placebo) # 7. 可视化 boxplot(drug, placebo, names=c("药物组","安慰剂组"), ylab="收缩压(mmHg)", main="两组收缩压比较", col=c("#0d9488","#e2e8f0")) # 8. 结果报告 cat("\n=== 结果报告 ===\n") cat("药物组收缩压", round(mean(drug),1), "±", round(sd(drug),1), "mmHg\n") cat("安慰剂组收缩压", round(mean(placebo),1), "±", round(sd(placebo),1), "mmHg\n") cat("t=", round(t_result$statistic,2), ", df=", t_result$parameter, ", p=", format(t_result$p.value, digits=4), "\n") cat("结论：两组差异具有统计学意义", ifelse(t_result$p.value<0.05,"(p<0.05)","(p>0.05)"),"\n")

📈 案例2：Logistic回归与预测模型

Logistic回归ROC曲线Nomogram

构建高血压风险预测模型，包含LASSO特征筛选、Logistic回归、ROC曲线和Nomogram。

# ===== 案例2：预测模型构建 ===== library(pROC) library(glmnet) set.seed(123) n <- 200 data <- data.frame( age = rnorm(n, 55, 12), bmi = rnorm(n, 26, 4), glucose = rnorm(n, 5.5, 1.2), smoking = rbinom(n, 1, 0.3), exercise = rbinom(n, 1, 0.5) ) # 生成结局 lp <- -3 + 0.04*data$age + 0.08*data$bmi + 0.3*data$glucose + 0.8*data$smoking - 0.5*data$exercise data$hypertension <- rbinom(n, 1, plogis(lp)) # LASSO特征筛选 x <- as.matrix(data[,1:5]) y <- data$hypertension cv_lasso <- cv.glmnet(x, y, alpha=1, family="binomial") coef(cv_lasso, s="lambda.min") # 多因素Logistic model <- glm(hypertension ~ ., data=data, family=binomial()) summary(model) OR <- exp(cbind(OR=coef(model), confint(model))) print(round(OR, 2)) # ROC曲线 pred <- predict(model, type="response") roc_obj <- roc(data$hypertension, pred) cat("AUC =", round(auc(roc_obj), 3), "\n") # 校准曲线 library(rms) dd <- datadist(data); options(datadist="dd") lrm_model <- lrm(hypertension ~ age+bmi+glucose+ smoking+exercise, data=data) cal <- calibrate(lrm_model, B=200) plot(cal)

🔬 案例3：生存分析完整流程

Kaplan-MeierCox回归PH检验

基于肺癌数据的生存分析，包含KM曲线、Log-rank检验、Cox回归和PH假设检验。

# ===== 案例3：生存分析 ===== library(survival) library(survminer) # 使用内置肺癌数据 data(lung) lung$sex <- factor(lung$sex, 1:2, c("Male","Female")) # 描述 summary(lung[, c("time","status","age","sex")]) # 1. KM曲线 km <- survfit(Surv(time, status) ~ sex, data=lung) ggsurvplot(km, data=lung, pval=TRUE, conf.int=TRUE, risk.table=TRUE, xlab="Days", ylab="Survival", palette=c("#0d9488","#06b6d4")) # 2. Log-rank检验 survdiff(Surv(time, status) ~ sex, data=lung) # 3. Cox回归 cox <- coxph(Surv(time, status) ~ sex + age + ph.ecog + wt.loss, data=lung) summary(cox) # 森林图 ggforest(cox, data=lung) # 4. PH假设检验 ph_test <- cox.zph(cox) print(ph_test) plot(ph_test, var="age") # 5. 预测生存率 surv_prob <- survfit(cox, newdata=data.frame( sex="Male", age=60, ph.ecog=1, wt.loss=5)) summary(surv_prob, times=c(100, 200, 365))

📊 案例4：Meta分析实操

Meta分析森林图异质性

基于模拟数据的Meta分析，包含固定/随机效应模型、森林图、漏斗图和发表偏倚检验。

# ===== 案例4：Meta分析 ===== library(meta) # 8项RCT研究数据 studies <- c("Smith 2019","Jones 2020", "Lee 2020","Wang 2021","Kim 2021", "Chen 2022","Liu 2022","Zhang 2023") events_t <- c(45,30,62,28,50,55,35,48) total_t <- c(100,80,120,60,95,100,75,90) events_c <- c(30,25,40,22,35,32,28,33) total_c <- c(100,82,115,58,90,95,72,88) # 固定效应 m_fixed <- metabin(events_t,total_t,events_c,total_c, studlab=studies, sm="OR", method="MH", fixed=TRUE, random=FALSE) cat("=== 固定效应 ===\n"); summary(m_fixed) # 随机效应 m_random <- metabin(events_t,total_t,events_c,total_c, studlab=studies, sm="OR", method="MH", random=TRUE, fixed=FALSE) cat("\n=== 随机效应 ===\n"); summary(m_random) # 森林图 forest(m_random, leftlabs=c("Study","Events,T","Total,T", "Events,C","Total,C")) # 异质性：I²=32.4%，p=0.17 → 低异质性 # 漏斗图+Egger检验 funnel(m_random) metabias(m_random, method.bias="linreg")

⚖️ 案例5：PSM倾向性评分匹配

PSM因果推断均衡性

评估药物治疗效果的观察性研究，使用倾向性评分匹配控制混杂因素。

# ===== 案例5：PSM分析 ===== library(MatchIt) library(cobalt) set.seed(123) n <- 1000 data <- data.frame( age = rnorm(n, 55, 12), sex = rbinom(n, 1, 0.5), bmi = rnorm(n, 26, 4), sbp = rnorm(n, 140, 20), chol = rnorm(n, 5.2, 1.1), treatment = rbinom(n, 1, plogis( -1 + 0.03*age + 0.3*sex + 0.05*bmi)) ) data$outcome <- rbinom(n, 1, plogis( -2 + 0.6*treatment + 0.04*age)) # PSM匹配 m.out <- matchit(treatment ~ age + sex + bmi + sbp + chol, data=data, method="nearest", ratio=1, caliper=0.2) # 匹配前后均衡性 love.plot(m.out, binary="std", title="PSM均衡性检验") # 匹配后效果 m.data <- match.data(m.out) table(m.data$treatment, m.data$outcome) model <- glm(outcome ~ treatment, data=m.data, family=binomial(), weights=weights) cat("OR =", round(exp(coef(model)[2]), 2), "\n95%CI:", round(exp(confint(model)[2,]), 2))

📱 公众号半年教程合集

Abel医研统计公众号26周系统性教程，覆盖从基础到前沿的60篇医学统计方法。点击文章标题展开全文。

👨‍⚕️ 博主信息

研究方向：医学统计学、临床流行病学、生物信息学、真实世界研究、临床预测模型、传染病数学建模

微信公众号：Abel医研统计

宗旨：让医学科研更简单，让统计方法更易懂

医学科研一站式学习平台

🌿 统计方法选择决策树

📊 医学统计学方法

P 值解读

置信区间 (CI)

正态性检验

方差齐性检验

独立样本t检验

配对t检验

单因素ANOVA

重复测量ANOVA

Mann-Whitney U检验

Wilcoxon符号秩检验

Kruskal-Wallis检验

卡方检验

Fisher精确检验

线性回归

Logistic回归

Poisson回归

Pearson相关

Spearman相关

ICC组内相关系数

Kappa一致性

诊断试验指标

ROC曲线

校准曲线

PCA主成分分析

EFA探索性因子分析

K-means聚类

层次聚类

中介分析

调节效应

⏳ 生存分析

Kaplan-Meier分析

Cox回归

竞争风险模型

🔧 缺失数据处理

缺失机制

多重插补MICE

📈 ROC曲线与预测模型评价

ROC曲线基础

ROC比较与校准

📈 时间序列分析

ARIMA模型

ETS指数平滑

Prophet时间序列预测

LSTM时间序列预测

VAR向量自回归

🧮 样本量计算（含R代码）

核心参数

ANOVA样本量

相关分析样本量

回归分析样本量

生存分析样本量

非劣效性样本量

t检验样本量

率的比较样本量

计算工具

⚖️ 倾向性评分与因果推断

PSM倾向性评分匹配

IPTW逆概率加权

DID双重差分

反事实因果模型

中断时间序列

工具变量

断点回归

有向无环图

🤖 医学机器学习

LASSO回归

随机森林

XGBoost

SVM支持向量机

神经网络

SHAP可解释性

🔮 贝叶斯统计

贝叶斯基础

贝叶斯t检验

贝叶斯回归

📋 Meta分析

Meta分析基础