775. 缺失值填补方法对比¶

MCAR（Missing Completely At Random，完全随机缺失）：
  缺失与任何已知/未知因素无关
  例如：样本处理时随机打翻了几个管子
  → 最好处理，直接删除或简单填补

MAR（Missing At Random，随机缺失）：
  缺失与已观测的变量有关，但与缺失值本身无关
  例如：低表达的基因在低测序深度样本中更容易缺失
  → 可用多重填补（MICE）

MNAR（Missing Not At Random，非随机缺失，最常见）：
  缺失与缺失值本身有关
  例如：蛋白质丰度太低 → 低于质谱检测限 → 缺失
  → 这是组学数据中最常见的！"低了所以测不到"
  → 需要特殊处理（左截断填补）

组学数据中缺失值的来源：
  RNA-seq: 低表达基因(count=0)——这不是缺失，是真的不表达
  蛋白组: 30-50%缺失率（低丰度蛋白检测不到）
  代谢组: 10-30%缺失率（低于检测限）
  16S: 稀有菌种在某些样本中reads=0

策略一：删除
  删除缺失比例高的特征（基因/蛋白）
  规则：缺失>50%的特征直接删
  优点：简单 | 缺点：丢失信息

策略二：简单填补
  用固定值替代（0、均值、中位数、最小值）
  优点：快 | 缺点：引入偏差

策略三：统计填补
  KNN、MICE、SVD等基于数据结构填补
  优点：利用数据关系 | 缺点：计算慢

策略四：左截断填补（组学专用）
  假设缺失值来自低于检测限的正态分布左尾
  适合MNAR机制的蛋白组/代谢组数据

# ===== 简单填补方法 =====
import pandas as pd  # 导入pandas
import numpy as np  # 导入numpy

# 读取含缺失值的数据
data = pd.read_csv("proteomics_data.csv", index_col=0)  # 蛋白组数据
print(f"数据维度: {data.shape}")  # 打印维度
print(f"缺失率: {data.isna().sum().sum() / data.size * 100:.1f}%")  # 缺失比例

# ===== Step 1: 过滤高缺失特征 =====
missing_ratio = data.isna().mean(axis=1)  # 每个特征的缺失比例
data_filtered = data[missing_ratio < 0.5]  # 保留缺失<50%的特征
print(f"过滤后: {data_filtered.shape} (删除了{data.shape[0]-data_filtered.shape[0]}个特征)")

# ===== 方法一：零值填补 =====
data_zero = data_filtered.fillna(0)  # 用0填充
# 适用：RNA-seq中0确实代表不表达
# 不适用：蛋白组（0不代表蛋白不存在，只是检测不到）

# ===== 方法二：均值填补 =====
data_mean = data_filtered.apply(
    lambda row: row.fillna(row.mean()), axis=1  # 每行用行均值填
)
# 优点：保持特征均值不变
# 缺点：降低方差，夸大显著性

# ===== 方法三：中位数填补 =====
data_median = data_filtered.apply(
    lambda row: row.fillna(row.median()), axis=1  # 每行用行中位数填
)
# 比均值更稳健（对异常值不敏感）

# ===== 方法四：最小值填补（蛋白组常用）=====
def min_value_impute(data, factor=0.5):
    """用最小值的一半填充（模拟低于检测限）"""
    imputed = data.copy()  # 复制数据
    for col in imputed.columns:  # 遍历样本
        min_val = imputed[col].min()  # 该样本的最小值
        imputed[col] = imputed[col].fillna(min_val * factor)  # 用最小值的factor倍
    return imputed

data_minval = min_value_impute(data_filtered)  # 最小值/2填补
# 适用：MNAR机制（缺失因为低于检测限）

# ===== KNN填补 =====
from sklearn.impute import KNNImputer  # 导入KNN填补器

# KNN原理：用最相似的K个样本的值来填补
knn_imputer = KNNImputer(
    n_neighbors=5,        # 使用5个最近邻
    weights="distance",   # 距离加权（近的样本权重大）
    metric="nan_euclidean" # 忽略NaN的欧氏距离
)

# 转置：KNN按行找近邻，所以特征在行
data_T = data_filtered.T  # 样本×特征 → 特征×样本
data_knn = pd.DataFrame(
    knn_imputer.fit_transform(data_T),  # KNN填补
    index=data_T.index,
    columns=data_T.columns
).T  # 转回来

print(f"KNN填补后缺失: {data_knn.isna().sum().sum()}")  # 应该=0

# ===== MICE (Multiple Imputation by Chained Equations) =====
from sklearn.experimental import enable_iterative_imputer  # 启用实验特性
from sklearn.impute import IterativeImputer  # 导入迭代填补器
from sklearn.ensemble import RandomForestRegressor  # 导入随机森林

# MICE原理：
# 对每个有缺失的变量，用其他变量建回归模型预测缺失值
# 迭代多轮直到收敛

# 方法一：线性回归MICE
mice_linear = IterativeImputer(
    max_iter=10,           # 最多迭代10轮
    random_state=42,       # 随机种子
    sample_posterior=True  # 抽样（多重填补）
)

data_mice = pd.DataFrame(
    mice_linear.fit_transform(data_filtered.T),  # 填补
    index=data_filtered.columns,
    columns=data_filtered.index
).T

# 方法二：随机森林MICE（更准但更慢）
mice_rf = IterativeImputer(
    estimator=RandomForestRegressor(
        n_estimators=50,     # 50棵树
        random_state=42      # 随机种子
    ),
    max_iter=5,              # 迭代5轮（RF慢，少迭代）
    random_state=42
)

data_mice_rf = pd.DataFrame(
    mice_rf.fit_transform(data_filtered.T),
    index=data_filtered.columns,
    columns=data_filtered.index
).T

# ===== 左截断填补（蛋白组/代谢组推荐）=====
# 假设：缺失值来自正态分布的左尾（低于检测限）

def left_censored_impute(data, shift=1.8, scale=0.3):
    """
    左截断填补（模拟DEP包的MinProb方法）
    shift: 向左移动多少个标准差（默认1.8）
    scale: 缩小标准差到多少倍（默认0.3）
    """
    imputed = data.copy()  # 复制数据

    for col in imputed.columns:  # 遍历每个样本
        observed = imputed[col].dropna()  # 观测值
        n_missing = imputed[col].isna().sum()  # 缺失数

        if n_missing > 0 and len(observed) > 0:
            mean_obs = observed.mean()  # 观测值均值
            std_obs = observed.std()  # 观测值标准差

            # 从下移+缩窄的正态分布中采样
            fill_mean = mean_obs - shift * std_obs  # 下移均值
            fill_std = std_obs * scale  # 缩小标准差

            np.random.seed(42)  # 固定种子
            fill_values = np.random.normal(
                fill_mean, fill_std, n_missing  # 从正态分布采样
            )

            imputed.loc[imputed[col].isna(), col] = fill_values  # 填充

    return imputed

data_leftcensor = left_censored_impute(data_filtered)
print("左截断填补完成")

# ===== 混合填补策略（推荐）=====
# 思路：MAR缺失用KNN，MNAR缺失用左截断
def hybrid_impute(data, mnar_threshold=0.3):
    """混合填补：根据缺失模式选择方法"""
    imputed = data.copy()

    for feature in imputed.index:  # 遍历特征
        missing_rate = imputed.loc[feature].isna().mean()  # 缺失率

        if missing_rate == 0:
            continue  # 无缺失跳过
        elif missing_rate > mnar_threshold:
            # 高缺失率 → 可能MNAR → 左截断填补
            observed = imputed.loc[feature].dropna()
            n_miss = imputed.loc[feature].isna().sum()
            fill_vals = np.random.normal(
                observed.mean() - 1.8 * observed.std(),
                observed.std() * 0.3,
                n_miss
            )
            imputed.loc[feature, imputed.loc[feature].isna()] = fill_vals
        else:
            # 低缺失率 → 可能MAR → KNN/中位数填补
            imputed.loc[feature] = imputed.loc[feature].fillna(
                imputed.loc[feature].median()
            )

    return imputed

# ===== 评估不同填补方法的效果 =====
import matplotlib.pyplot as plt  # 导入matplotlib
from sklearn.decomposition import PCA  # 导入PCA
from sklearn.metrics import mean_squared_error  # 导入MSE

# ===== 方法一：人为制造缺失来评估 =====
def evaluate_imputation(complete_data, method_func, missing_rate=0.1, seed=42):
    """
    在完整数据上人为制造缺失，评估填补准确性
    """
    np.random.seed(seed)  # 设置种子

    # 制造缺失
    mask = np.random.random(complete_data.shape) < missing_rate  # 随机掩码
    masked_data = complete_data.copy()  # 复制
    masked_data.values[mask] = np.nan  # 设为缺失

    # 填补
    imputed_data = method_func(masked_data)  # 执行填补

    # 计算误差（只在人为制造的缺失位置）
    rmse = np.sqrt(mean_squared_error(
        complete_data.values[mask],  # 真实值
        imputed_data.values[mask]   # 填补值
    ))

    # 计算相关性
    corr = np.corrcoef(
        complete_data.values[mask].flatten(),
        imputed_data.values[mask].flatten()
    )[0, 1]

    return {"RMSE": round(rmse, 4), "Correlation": round(corr, 4)}

# ===== 方法二：PCA对比可视化 =====
def visualize_imputation_pca(original, methods_dict, output="imputation_pca.png"):
    """用PCA可视化不同填补方法的效果"""
    fig, axes = plt.subplots(1, len(methods_dict), figsize=(5*len(methods_dict), 4))

    # 原始数据PCA
    pca = PCA(n_components=2)  # 2维PCA
    pca.fit(original.T.dropna())  # 在完整数据上fit

    for i, (name, imputed) in enumerate(methods_dict.items()):
        ax = axes[i] if len(methods_dict) > 1 else axes
        scores = pca.transform(imputed.T)  # 投影
        ax.scatter(scores[:, 0], scores[:, 1], alpha=0.7)  # 散点图
        ax.set_title(f"{name}")  # 方法名
        ax.set_xlabel("PC1")  # x轴
        ax.set_ylabel("PC2")  # y轴

    plt.tight_layout()
    plt.savefig(output, dpi=300)  # 保存
    plt.close()

# 对比所有方法
print("=== 填补方法对比 ===")
methods = {
    "Zero": lambda d: d.fillna(0),
    "Median": lambda d: d.apply(lambda r: r.fillna(r.median()), axis=1),
    "KNN": lambda d: pd.DataFrame(
        KNNImputer(n_neighbors=5).fit_transform(d.T),
        index=d.columns, columns=d.index
    ).T,
    "LeftCensor": left_censored_impute
}

for name, func in methods.items():
    result = evaluate_imputation(data_filtered.dropna(), func)
    print(f"{name:12s}: RMSE={result['RMSE']:.4f}, Corr={result['Correlation']:.4f}")

# ===== 缺失值填补速查 =====

# 检查缺失
data.isna().sum()                    # 每列缺失数
data.isna().mean()                   # 每列缺失比例

# 简单填补
data.fillna(0)                       # 零填补
data.fillna(data.median())           # 中位数填补
data.fillna(data.min() * 0.5)        # 最小值/2填补

# KNN填补
from sklearn.impute import KNNImputer
imputer = KNNImputer(n_neighbors=5)
data_filled = imputer.fit_transform(data)

# MICE填补
from sklearn.impute import IterativeImputer
imputer = IterativeImputer(max_iter=10, random_state=42)
data_filled = imputer.fit_transform(data)

# 方法选择
# RNA-seq零值 → 不填补（真实信号）
# 蛋白组MNAR → 左截断填补(MinProb)
# 低缺失率MAR → KNN或MICE
# 混合缺失 → 混合策略(高缺失左截断+低缺失KNN)
# 统计检验 → 多重填补(MICE, M=5-10)