408_跨队列多组学验证¶

import pandas as pd
import numpy as np
from sklearn.preprocessing import StandardScaler
from sklearn.linear_model import LogisticRegression
from sklearn.metrics import roc_auc_score, roc_curve
import matplotlib.pyplot as plt

np.random.seed(42)

# ========== 1. 模拟发现队列（训练数据） ==========
def simulate_cohort(n_samples, n_features, effect_size=1.5, seed=0):
    """模拟一个多组学队列数据"""
    np.random.seed(seed)
    X = np.random.randn(n_samples, n_features)  # 特征矩阵
    y = np.random.binomial(1, 0.4, n_samples)   # 二分类结局（0=对照，1=病例）

    # 前20个特征有真实的预测能力
    for i in range(20):
        X[y == 1, i] += effect_size  # 病例组特征均值更高

    return X, y

# 发现队列（n=200）
X_discovery, y_discovery = simulate_cohort(200, 100, effect_size=1.5, seed=0)
print(f"发现队列: {X_discovery.shape}, 病例数: {y_discovery.sum()}")

# ========== 2. 在发现队列中训练模型 ==========
scaler = StandardScaler()
X_discovery_scaled = scaler.fit_transform(X_discovery)  # 注意：scaler在发现队列上fit

# Lasso逻辑回归（稀疏特征选择）
from sklearn.linear_model import LassoCV
from sklearn.pipeline import Pipeline

model = LogisticRegression(
    penalty="l1",        # Lasso正则化，自动特征选择
    C=0.1,               # 正则化强度（越小越稀疏）
    solver="liblinear",
    random_state=42
)
model.fit(X_discovery_scaled, y_discovery)

# 内部AUC
y_pred_discovery = model.predict_proba(X_discovery_scaled)[:, 1]
auc_discovery = roc_auc_score(y_discovery, y_pred_discovery)
print(f"发现队列内部AUC: {auc_discovery:.3f}")

# ========== 3. 生成独立验证队列（来自不同"医院"） ==========
# 验证队列1：效应量相同，但有轻微批次效应
X_val1, y_val1 = simulate_cohort(150, 100, effect_size=1.3, seed=1)
X_val1 += np.random.randn(*X_val1.shape) * 0.2  # 加入批次噪声

# 验证队列2：不同种群，效应量略小
X_val2, y_val2 = simulate_cohort(120, 100, effect_size=1.0, seed=2)

# ========== 4. 跨队列验证（关键：用发现队列的scaler变换验证集！） ==========
# 错误做法：在验证集上重新fit scaler（数据泄露！）
# 正确做法：用发现队列fit好的scaler转换验证集

X_val1_scaled = scaler.transform(X_val1)  # 用已有scaler，不重新fit！
X_val2_scaled = scaler.transform(X_val2)

# 在验证队列上计算AUC
y_pred_val1 = model.predict_proba(X_val1_scaled)[:, 1]
y_pred_val2 = model.predict_proba(X_val2_scaled)[:, 1]

auc_val1 = roc_auc_score(y_val1, y_pred_val1)
auc_val2 = roc_auc_score(y_val2, y_pred_val2)

print(f"验证队列1 AUC: {auc_val1:.3f}")
print(f"验证队列2 AUC: {auc_val2:.3f}")

# ========== 5. 批次效应校正（ComBat风格） ==========
def simple_combat(data, batch_labels):
    """
    简化版ComBat批次效应校正
    实际使用 combat-seq（R）或 pycombat（Python）
    """
    data_corrected = data.copy()
    unique_batches = np.unique(batch_labels)

    # 计算总体均值和标准差
    global_mean = data.mean(axis=0)
    global_std = data.std(axis=0)

    for batch in unique_batches:
        mask = batch_labels == batch
        # 计算该批次的均值和标准差
        batch_mean = data[mask].mean(axis=0)
        batch_std = data[mask].std(axis=0) + 1e-8  # 防止除零

        # 将该批次数据标准化到总体分布
        data_corrected[mask] = (data[mask] - batch_mean) / batch_std * global_std + global_mean

    return data_corrected

# 模拟合并两个验证队列并校正批次效应
X_combined = np.vstack([X_val1, X_val2])
batch_labels = np.array([0] * 150 + [1] * 120)  # 批次标签
X_corrected = simple_combat(X_combined, batch_labels)

print(f"\n批次效应校正前后对比:")
print(f"校正前 队列1均值: {X_combined[:150, 0].mean():.3f}")
print(f"校正后 队列1均值: {X_corrected[:150, 0].mean():.3f}")

# ========== 6. 绘制多队列ROC曲线对比图 ==========
fig, ax = plt.subplots(figsize=(8, 7))
datasets = [
    (y_discovery, y_pred_discovery, f"Discovery (AUC={auc_discovery:.2f})"),
    (y_val1, y_pred_val1, f"Validation-1 (AUC={auc_val1:.2f})"),
    (y_val2, y_pred_val2, f"Validation-2 (AUC={auc_val2:.2f})")
]
colors_plot = ["#E74C3C", "#2ECC71", "#3498DB"]

for (y_true, y_score, label), color in zip(datasets, colors_plot):
    fpr, tpr, _ = roc_curve(y_true, y_score)
    ax.plot(fpr, tpr, color=color, lw=2, label=label)

ax.plot([0, 1], [0, 1], "k--", lw=1.5, label="Random")
ax.set_xlabel("False Positive Rate")
ax.set_ylabel("True Positive Rate")
ax.set_title("Cross-cohort Validation ROC Curves")
ax.legend(loc="lower right")
plt.tight_layout()
plt.savefig("cross_cohort_roc.png", dpi=150)
print("\nROC曲线已保存")