医学预测模型怎么做：LASSO、随机森林与 XGBoost、列线图与外部验证

library(pmsampsize)   # 1.1.3
# 二分类结局：患病率 0.253，参考已发表模型的 C 统计量 0.80（尽量用优化度校正后的值），30 个候选参数
pmsampsize(type = "b", cstatistic = 0.80, parameters = 30, prevalence = 0.253)
#> Criteria 1 (shrinkage >= 0.9)          1136
#> Criteria 2 (R2 optimism <= 0.05)        774
#> Criteria 3 (intercept +/- 0.05)         291
#> Final 1136, 288 events, EPP = 9.58

pmsampsize(type = "b", cstatistic = 0.70, parameters = 30, prevalence = 0.253)
#> Final 2732, 692 events, EPP = 23.04

# 生存结局：Cox-Snell R2 0.10，15 个参数，事件率 0.08/人年，预测 3 年，平均随访 4 年
pmsampsize(type = "s", csrsquared = 0.10, parameters = 15, rate = 0.08,
           timepoint = 3, meanfup = 4)
#> Final 1274, 408 events, EPP = 27.18

library(glmnet)                                     # 本页实测版本 5.0
X <- model.matrix(y ~ ., data = dat)[, -1]          # 因子自动展开为哑变量；X 必须是数值矩阵
y <- dat$y                                          # 0/1

set.seed(2026)
foldid <- sample(rep(1:10, length.out = nrow(X)))   # 固定折分并保存，复现和补充材料都用它
cv <- cv.glmnet(X, y, family = "binomial", foldid = foldid)
# 默认 alpha = 1（LASSO）、standardize = TRUE（内部标准化，系数仍按原尺度返回）、type.measure = "deviance"
plot(cv)                                            # 左竖线 lambda.min，右竖线 lambda.1se
coef(cv, s = "lambda.min")                          # 一定写 s；coef(cv) 不写 s 时返回的是 lambda.1se
coef(cv, s = "lambda.1se")

# 变量集稳定性：换 100 个种子各做一次交叉验证，统计每个变量的入选次数
sel <- sapply(1:100, function(s) {
  set.seed(s)
  b <- coef(cv.glmnet(X, y, family = "binomial", nfolds = 10), s = "lambda.1se")[-1, 1]
  b != 0
})
sort(rowSums(sel), decreasing = TRUE)               # 入选频率
length(unique(apply(sel, 2, paste, collapse = ""))) # 不同变量集的个数

# 降低折分随机性：同一 lambda 序列下平均 100 次折分的交叉验证误差曲线
lam <- glmnet(X, y, family = "binomial")$lambda
cvm <- sapply(1:100, function(s) {
  set.seed(s)
  cv.glmnet(X, y, family = "binomial", lambda = lam,
            foldid = sample(rep(1:10, length.out = nrow(X))))$cvm
})
lam[which.min(rowMeans(cvm))]                       # 平均曲线上的 lambda.min

import numpy as np
from sklearn.model_selection import StratifiedKFold, GridSearchCV, cross_validate
from sklearn.pipeline import Pipeline
from sklearn.impute import SimpleImputer
from sklearn.preprocessing import StandardScaler
from sklearn.linear_model import LogisticRegression
from sklearn.ensemble import RandomForestClassifier
from xgboost import XGBClassifier

inner = StratifiedKFold(n_splits=5, shuffle=True, random_state=1)   # 内层：调参
outer = StratifiedKFold(n_splits=5, shuffle=True, random_state=2)   # 外层：估计性能
models = {
    # scikit-learn 1.8 起 penalty 参数弃用，L1 写作 l1_ratio=1；1.8 之前写 penalty="l1"
    "lasso_lr": (Pipeline([("imp", SimpleImputer(strategy="median")),
                           ("sc", StandardScaler()),
                           ("clf", LogisticRegression(l1_ratio=1, solver="saga", max_iter=5000))]),
                 {"clf__C": np.logspace(-3, 1, 20)}),
    "rf": (Pipeline([("imp", SimpleImputer(strategy="median")),
                     ("clf", RandomForestClassifier(n_estimators=500, random_state=0))]),
           {"clf__max_features": ["sqrt", 0.33, 0.5], "clf__min_samples_leaf": [1, 5, 10, 20]}),
    "xgb": (Pipeline([("clf", XGBClassifier(n_estimators=300, eval_metric="logloss"))]),
            {"clf__learning_rate": [0.03, 0.1], "clf__max_depth": [2, 3, 4],
             "clf__min_child_weight": [1, 5]}),
}
for name, (pipe, grid) in models.items():
    search = GridSearchCV(pipe, grid, cv=inner, scoring="neg_log_loss")  # 用对数损失调参，兼顾校准
    res = cross_validate(search, X, y, cv=outer, scoring=["roc_auc", "neg_brier_score"])
    print(name, res["test_roc_auc"].mean().round(3), res["test_roc_auc"].std().round(3),
          (-res["test_neg_brier_score"]).mean().round(4))

library(glmnet); library(pROC)
auc_of <- function(p, y) as.numeric(auc(roc(y, p, direction = "<", quiet = TRUE)))  # 显式 direction
fit_lasso <- function(X, y) cv.glmnet(X, y, family = "binomial", nfolds = 10)
p_of <- function(cv, X) predict(cv, X, s = "lambda.1se", type = "response")[, 1]

cv0 <- fit_lasso(X, y)
apparent <- auc_of(p_of(cv0, X), y)                 # 表观 AUC：同一数据建模又评估
set.seed(2026)
opt <- replicate(200, {
  i <- sample(nrow(X), replace = TRUE)
  cvb <- fit_lasso(X[i, ], y[i])                    # 变量筛选与 lambda 选择在每个 bootstrap 样本内重做
  auc_of(p_of(cvb, X[i, ]), y[i]) - auc_of(p_of(cvb, X), y)
})
c(apparent = apparent, optimism = mean(opt), corrected = apparent - mean(opt))

from imblearn.pipeline import Pipeline          # 用 imblearn 的 Pipeline，sklearn 的 Pipeline 不接受采样器
from imblearn.over_sampling import SMOTE
from sklearn.preprocessing import StandardScaler
from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import StratifiedKFold, cross_val_score

pipe = Pipeline([("sc", StandardScaler()),
                 ("smote", SMOTE(random_state=0)),   # 只作用于每个训练折，测试折保持原始分布
                 ("clf", LogisticRegression(max_iter=1000))])
cv = StratifiedKFold(n_splits=5, shuffle=True, random_state=0)
print(cross_val_score(pipe, X, y, cv=cv, scoring="roc_auc").mean())
# 重采样后的预测概率整体偏高，要报告风险时需在原始分布的数据上重新校准（TRIPOD+AI 第 13 条）

library(pROC); library(rms); library(dcurves)
# p：外部验证集（或外层交叉验证、bootstrap 校正后）的预测概率；y：0/1 结局
r <- roc(y, p, direction = "<", quiet = TRUE)    # 验证集上显式写 direction
ci.auc(r)                                        # 默认 DeLong 法 95% CI
# roc.test(r1, r2)                               # 两个非嵌套模型在同一人群上的配对 DeLong 检验

val.prob(p, y)                                   # 输出 C、Brier、校准截距、校准斜率、Emax 并画平滑校准曲线

df <- data.frame(y = y, model = p)
dca(y ~ model, data = df, thresholds = seq(0, 0.5, by = 0.01)) |>
  plot(smooth = TRUE)                            # 同时画出 Treat All 与 Treat None

library(rms)                                     # 本页核对版本 8.1.1
dd <- datadist(dat); options(datadist = "dd")    # 必须在拟合前设置，且用同一个数据框
fit <- lrm(y ~ rcs(age, 4) + sex + scoma + meanbp + pafi + alb + bili + crea,
           data = dat, x = TRUE, y = TRUE)       # x、y = TRUE 供 validate 与 calibrate 使用
nom <- nomogram(fit, fun = plogis, funlabel = "Risk of in-hospital death",
                fun.at = c(0.05, 0.1, 0.2, 0.4, 0.6, 0.8))
plot(nom)

set.seed(2026)
validate(fit, B = 200)       # Dxy、斜率等的优化度校正；C = Dxy / 2 + 0.5；B 默认只有 40
cal <- calibrate(fit, B = 200)
plot(cal)                    # Apparent、Bias-corrected 两条线；报告 Bias-corrected

library(survival); library(glmnet); library(rms); library(timeROC)
d <- na.omit(subset(pbc, !is.na(trt),
     select = c(time, status, age, albumin, bili, copper, protime, ast, edema, stage)))
d$death <- as.integer(d$status == 2)             # 2 = 死亡；1 = 肝移植，按删失处理
d$status <- NULL

# LASSO-Cox：y 直接用 Surv 对象（glmnet 4.1 起支持）
X <- model.matrix(~ . - time - death, d)[, -1]
set.seed(2026)
cvc <- cv.glmnet(X, Surv(d$time, d$death), family = "cox", type.measure = "C")
coef(cvc, s = "lambda.1se")

# Cox 列线图：surv = TRUE 且 time.inc 等于之后 calibrate 的 u
dd <- datadist(d); options(datadist = "dd")
fit <- cph(Surv(time, death) ~ log(bili) + albumin + age + protime + edema,
           data = d, x = TRUE, y = TRUE, surv = TRUE, time.inc = 1826)
cox.zph(fit)                                     # 比例风险假定
sv <- Survival(fit)
plot(nomogram(fit, fun = list(function(lp) sv(1826, lp)), funlabel = "5-year survival"))
set.seed(2026)
cal <- calibrate(fit, u = 1826, cmethod = "KM", m = 50, B = 200)   # 默认 cmethod = "hare" 需要 polspline 包
plot(cal)

# 时间依赖 ROC 与 C-index
lp <- predict(fit)                               # 线性预测值，越大风险越高
tr <- timeROC(T = d$time, delta = d$death, marker = lp, cause = 1,
              times = c(730, 1826), iid = TRUE)
tr$AUC; confint(tr)
concordance(Surv(time, death) ~ lp, data = d, reverse = TRUE)                   # Harrell C
concordance(Surv(time, death) ~ lp, data = d, reverse = TRUE, timewt = "n/G2")  # Uno C

import shap                                      # 0.52
explainer = shap.TreeExplainer(xgb_model)        # xgb_model：已拟合的 XGBClassifier
sv = explainer(X_test)                           # 在测试集或外层折上计算，不用训练集
shap.plots.beeswarm(sv)                          # 二分类时横轴单位是 log-odds，不是概率
shap.plots.scatter(sv[:, "age"], color=sv)       # 依赖图：看非线性与交互