一、赛题背景与核心理解
2026 年江西省研究生数学建模竞赛赛题 3 的题目是“电子健康记录数据补全及其优化算法”。这是一道非常典型的医学数据建模题，核心不是简单地把缺失值填上，而是要围绕电子健康记录数据的缺失机制、数据稀疏性、矩阵补全方法、主成分分析和多算法对比，建立一套完整、可靠、可解释的 EHR 缺失补全与临床指标挖掘方案。
电子健康记录，即 EHR，通常包含患者在住院或随访过程中的生命体征、实验室检查、人口学信息、病情标签等。相比问卷调查数据，EHR 更贴近真实临床过程，但它也有一个非常突出的特点：缺失非常严重。很多指标不是每个小时都测，有些实验室指标只有病情变化时才测，有些患者住院时间短，有些患者记录周期长，这会导致同一个数据集中不同患者、不同变量、不同时间点的缺失程度差异巨大。
本题提供的数据来自类似 CINC2019 的临床时间序列数据。每位患者对应一个 .psv 文件，文件中的每一行代表一次时间观测，每一列代表一个生命体征或临床变量。题目要求选取标准化子集开展建模，重点完成：
缺失率统计、变量筛选与归一化；
简单矩阵缺失值补全；
基于 NMF 的低秩矩阵补全；
基于补全数据的 PCA 主成分分析；
根据 PCA 结果筛选关键生命体征指标；
在不同缺失比例下比较 SVD-PCA、NMF 及其他算法的补全精度；
综合误差、稳定性和复杂度筛选最优 EHR 数据补全模型。
这道题的本质是：
在高缺失、强异质、不同量纲的临床 EHR 数据中，构建可靠的数据补全模型，并验证补全结果是否能提升后续 PCA 分析和关键生理指标挖掘的有效性。
高分论文不能只写“用均值填补、NMF 填补、PCA 分析”。真正好的方案必须回答以下问题：
为什么 EHR 数据会缺失？
缺失是随机缺失还是非随机缺失？
哪些变量缺失率过高，应当剔除？
归一化应在什么时候做？
NMF 为什么适合非负生命体征数据？
分解秩 r 如何选择？
如何避免 NMF 过拟合？
补全后 PCA 的主成分是否更稳定？
补全精度如何在不同缺失率下验证？
最终哪个算法更适合 EHR 数据？
二、整体建模主线
本题建议采用“五层递进式建模框架”：
第一层：数据读取与缺失模式分析。
从 .psv 文件读取患者 EHR 矩阵，统计每个变量、每个患者、每个时间点的缺失率，识别高缺失变量和可用变量。
第二层：数据标准化与低维补全。
对筛选出的有效变量进行 Min-Max 归一化或 Z-score 标准化，标记 NaN 缺失位置，并对题目给出的小矩阵采用均值填补、KNN 填补、SVD 填补或 NMF 填补进行示范。
第三层：NMF 低秩矩阵补全。
将清洗后的 EHR 数据构造成观测矩阵 (A)，使用带缺失掩码的 NMF 模型 (A \\approx WH)，通过迭代更新非负矩阵 (W) 与 (H)，完成缺失点填补。
第四层：PCA 主成分分析与临床规律挖掘。
对 NMF 补全后的完整数据进行 Z-score 标准化，计算协方差矩阵，提取第一、第二主成分，计算方差贡献率和指标载荷，筛选关键生命体征指标。
第五层：多算法、多缺失率补全精度对比。
以补全后的完整数据作为基准，在 10%、20%、30%、40%、50% 五档缺失率下随机掩码，分别用 SVD-PCA、NMF、KNN、MICE、SoftImpute 等方法补全，计算 RSE，并从精度、稳定性和复杂度三方面选择最优模型。
最终论文主线可以概括为：
缺失识别 → 变量筛选 → 矩阵补全 → PCA 降维 → 关键指标筛选 → 多算法鲁棒比较 → 最优补全模型推荐。
三、任务一：数据预处理、缺失率计算、归一化与低维矩阵补全
3.1 数据读取与字段理解
每位患者的 EHR 数据是一个 .psv 文件，字段通常包括：
HR：心率
O2Sat：血氧饱和度
Temp：体温
SBP：收缩压
MAP：平均动脉压
DBP：舒张压
Resp：呼吸频率
EtCO2：呼气末二氧化碳
BaseExcess：碱剩余
HCO3：碳酸氢根
FiO2：吸入氧浓度
pH：酸碱度
PaCO2：动脉二氧化碳分压
SaO2：动脉血氧饱和度
AST、BUN、Creatinine、Glucose、Lactate 等实验室指标
Age、Gender、Unit1、Unit2、HospAdmTime、ICULOS 等基础特征
SepsisLabel：脓毒症标签
题目要求对 p000001.psv 的前 40 个生命特征以及缺失信息进行统计，并选择前 36 个特征变量，找出其中缺失率小于 90% 的前 11 个特征标量。这里的“前 36 个特征变量”一般可以理解为剔除相对稳定的基础变量或标签后的动态生命体征变量，重点关注随时间变化的临床指标。
3.2 缺失率计算
对某一特征 (j)，缺失率定义为：
[
MR_j=\\frac{\\sum_{i=1}^{n}I(x_{ij}\\text{ is NaN})}{n}
]
其中 (n) 是该患者的观测时间点数量，(I) 是指示函数。
若 (MR_j < 0.9)，说明该变量至少有 10% 的有效观测，可以保留进入后续建模。
对每个患者也可以计算整体缺失率：
[
MR_i=\\frac{\\sum_{j=1}^{p}I(x_{ij}\\text{ is NaN})}{p}
]
该指标可用于筛选观测质量较高的患者样本。
3.3 为什么要剔除高缺失变量？
EHR 中某些实验室指标缺失率极高，例如 TroponinI、Fibrinogen、Bilirubin_direct 等指标可能并不是所有患者都检测。若直接保留，会产生两个问题：
第一，补全结果主要由模型猜测得到，可信度很低；
第二，这些变量会对 PCA 和距离计算产生不稳定影响。
因此，缺失率过高的变量不应盲目补全。合理做法是：
缺失率低于 30%：优先保留；
缺失率 30%–70%：视临床意义和模型需要选择保留；
缺失率 70%–90%：谨慎使用；
缺失率超过 90%：一般剔除。
题目明确要求筛选缺失率小于 90% 的变量，这是为了确保后续补全具有基本数据支撑。
3.4 归一化方法
生命体征变量量纲不同。比如心率单位是次/分钟，体温单位是摄氏度，血压单位是 mmHg，pH 范围约 7 左右。若不归一化，数值范围大的变量会在 NMF、SVD、PCA 中占据不合理权重。
常用归一化方法有两种。
Min-Max 归一化
[
x'=\\frac{x-x_{\\min}}{x_{\\max}-x_{\\min}}
]
优点是将数据压缩到 ([0,1])，适合 NMF，因为 NMF 要求非负输入。
Z-score 标准化
[
z=\\frac{x-\\mu}{\\sigma}
]
优点是适合 PCA，因为 PCA 基于方差结构，Z-score 可以消除量纲影响。
建议流程是：
NMF 补全前使用 Min-Max 归一化，保证非负性；
PCA 分析前使用 Z-score 标准化，保证协方差分析公平。
3.5 任务一代码：读取 psv 文件、统计缺失率、筛选变量
import os
import zipfile
import numpy as np
import pandas as pd
from io import TextIOWrapper
ZIP_PATH = "附件2_data_cinc2019.zip"
def read_psv_from_zip(zip_path, inner_path):
with zipfile.ZipFile(zip_path) as z:
with z.open(inner_path) as f:
return pd.read_csv(TextIOWrapper(f), sep="|")
def list_training_setA_files(zip_path):
with zipfile.ZipFile(zip_path) as z:
names = [
n for n in z.namelist()
if n.startswith("cinc2019/training/training_setA/")
and n.endswith(".psv")
and "__MACOSX" not in n
]
return sorted(names)
files_A = list_training_setA_files(ZIP_PATH)
print("training_setA 文件数量:", len(files_A))
print("前5个文件:", files_A[:5])
df1 = read_psv_from_zip(ZIP_PATH, "cinc2019/training/training_setA/p000001.psv")
print(df1.head())
print(df1.shape)
# 前 40 个生命特征，不含 SepsisLabel
feature_cols_40 = df1.columns[:40]
missing_rate = df1[feature_cols_40].isna().mean().sort_values()
missing_table = pd.DataFrame({
"feature": missing_rate.index,
"missing_rate": missing_rate.values
})
print(missing_table)
# 前 36 个特征中缺失率小于 90% 的前 11 个
feature_cols_36 = df1.columns[:36]
mr36 = df1[feature_cols_36].isna().mean()
selected_11 = mr36[mr36 < 0.90].sort_values().head(11).index.tolist()
print("缺失率小于90%的前11个变量:")
print(selected_11)
3.6 选取 50 份 36×41 规格患者数据矩阵
题目要求从 training_setA 中筛选 50 份 36×41 规格患者临床数据矩阵。这里的 36 行表示选择前 36 个时间点，41 列表示完整字段，包括 40 个特征和 SepsisLabel。
若某患者记录长度不足 36，则不选；若长度超过 36，则取前 36 行或按任务要求取标准化片段。为了保证比较公平，可以统一取每位患者的前 36 条记录。
def load_standardized_patients(zip_path, n_patients=50, n_rows=36, n_cols=41):
files = list_training_setA_files(zip_path)
matrices = []
names = []
for path in files:
df = read_psv_from_zip(zip_path, path)
if df.shape[0] >= n_rows and df.shape[1] >= n_cols:
mat = df.iloc[:n_rows, :n_cols].copy()
matrices.append(mat)
names.append(os.path.basename(path))
if len(matrices) >= n_patients:
break
return names, matrices
names50, mats50 = load_standardized_patients(ZIP_PATH, n_patients=50)
print("选取患者数:", len(mats50))
print("第一个矩阵形状:", mats50[0].shape)
3.7 标记 NaN 缺失位置
缺失掩码矩阵 (M) 定义为：
[
M_{ij}=
\\begin{cases}
1, & A_{ij}\\text{ 已观测}\\
0, & A_{ij}\\text{ 缺失}
\\end{cases}
]
后续 NMF 和误差计算都要依赖这个掩码。
def nan_mask(matrix_df):
A = matrix_df.values.astype(float)
observed_mask = ~np.isnan(A)
missing_mask = np.isnan(A)
return A, observed_mask, missing_mask
A0, obs0, miss0 = nan_mask(mats50[0])
print("缺失数量:", miss0.sum())
print("观测数量:", obs0.sum())
3.8 小矩阵缺失补全示例
题目给出一个含 NaN 的小矩阵，要求选取两种以上方法完成补全。这里可展示三种：均值填补、KNN 填补、低秩 SVD 填补。
方法一：列均值填补
特点：简单、快速、稳定，但忽略变量间相关性。
方法二：KNN 填补
特点：利用相似样本信息，适合局部结构明显的数据，但对距离度量敏感。
方法三：低秩 SVD 填补
特点：利用矩阵低秩结构，适合变量之间存在潜在相关结构的数据。
from sklearn.impute import SimpleImputer, KNNImputer
from sklearn.experimental import enable_iterative_imputer
from sklearn.impute import IterativeImputer
A_small = np.array([
[5, 3, 8, 1],
[9, 7, 4, np.nan],
[6, 5, 2, 9],
[1, 8, 3, 7],
[4, 6, 2, 5],
[8, 1, 6, np.nan],
[np.nan, np.nan, np.nan, np.nan],
[np.nan, np.nan, np.nan, np.nan],
[np.nan, np.nan, np.nan, np.nan]
], dtype=float)
# 1. 均值填补
mean_imputer = SimpleImputer(strategy="mean")
A_mean = mean_imputer.fit_transform(A_small)
# 2. KNN 填补
knn_imputer = KNNImputer(n_neighbors=2)
A_knn = knn_imputer.fit_transform(A_small)
# 3. MICE/迭代回归填补
iter_imputer = IterativeImputer(random_state=2026, max_iter=20)
A_iter = iter_imputer.fit_transform(A_small)
print("均值填补:")
print(np.round(A_mean, 3))
print("KNN填补:")
print(np.round(A_knn, 3))
print("迭代回归填补:")
print(np.round(A_iter, 3))
论文中要说明：低维小矩阵示例只是展示不同补全思想，真正高维 EHR 数据应采用能利用低秩结构和变量相关性的模型，例如 NMF、SVD-PCA 或 MICE。
四、任务二：NMF 缺失填充与 PCA 主成分分析建模
4.1 为什么选择 NMF？
NMF，即非负矩阵分解，假设非负数据矩阵 (A) 可以近似分解为：
[
A \\approx WH
]
其中：
[
W \\ge 0,\\quad H \\ge 0
]
(W) 表示样本在潜在生理模式上的权重，(H) 表示每个潜在模式在不同生命体征上的贡献。
NMF 适合本题有三个原因：
第一，生命体征经 Min-Max 归一化后为非负数据，符合 NMF 非负约束。
第二，EHR 变量之间存在潜在低秩结构，例如血压、心率、呼吸、血氧等可能共同反映循环和呼吸状态。
第三，NMF 分解具有一定可解释性，每个潜在因子可以看作某类生理状态模式。
4.2 带缺失掩码的 NMF 目标函数
由于原始矩阵存在缺失，不能直接对 NaN 计算误差。应只在观测位置上拟合：
[
\\min_{W,H\\ge0}\\frac{1}{2}|M\\odot(A-WH)|_F^2+\\lambda(|W|_F^2+|H|_F^2)
]
其中 (M) 是观测掩码矩阵，(\\odot) 表示逐元素乘积。
这个模型的含义是：只要求 (WH) 在原本有观测值的位置尽量接近真实值，而缺失位置由低秩结构自动估计。
4.3 分解秩 r 的选择
分解秩 (r) 决定潜在生理模式数量。若 (r) 太小，模型欠拟合；若 (r) 太大，模型可能过拟合噪声。
可以用以下方法选择 (r)：
交叉验证：随机遮盖一部分已知值，比较不同 (r) 的 RSE；
重构误差拐点法：观察误差随 (r) 增加的下降曲线；
临床可解释性：通常选择 3–10 个潜在生理模式更易解释；
稳定性分析：多次初始化后结果波动小的 (r) 更可靠。
建议论文中选用交叉验证法确定 (r)，并给出误差曲线。
4.4 NMF 缺失补全代码
import numpy as np
def masked_nmf(A, mask, rank=5, max_iter=500, lr=0.005, reg=1e-4, random_state=2026):
"""
带缺失掩码的 NMF。
A: 已归一化矩阵，缺失位置可先填0
mask: 观测位置为1，缺失位置为0
"""
rng = np.random.default_rng(random_state)
n, p = A.shape
W = rng.random((n, rank))
H = rng.random((rank, p))
losses = []
A_filled = np.nan_to_num(A, nan=0.0)
M = mask.astype(float)
for it in range(max_iter):
WH = W @ H
E = M * (WH – A_filled)
grad_W = E @ H.T + reg * W
grad_H = W.T @ E + reg * H
W -= lr * grad_W
H -= lr * grad_H
W = np.maximum(W, 1e-8)
H = np.maximum(H, 1e-8)
if it % 20 == 0:
loss = 0.5 * np.sum(E ** 2) + 0.5 * reg * (np.sum(W ** 2) + np.sum(H ** 2))
losses.append(loss)
A_hat = W @ H
A_completed = A_filled.copy()
A_completed[~mask] = A_hat[~mask]
return A_completed, W, H, losses
4.5 构造 50 份患者合并矩阵
题目可以有两种理解方式：
第一种，把每个患者 36×41 矩阵展平为一个样本向量，形成 50×1476 的矩阵；
第二种，把 50 个患者的所有时间点拼接，形成 1800×41 的矩阵。
更推荐第二种，因为 PCA 和 NMF 面向生命体征变量，保留列含义更清楚。
def concat_patient_matrices(mats, use_cols=None):
frames = []
for df in mats:
if use_cols is not None:
frames.append(df[use_cols])
else:
frames.append(df)
return pd.concat(frames, axis=0, ignore_index=True)
# 只使用前40个生命体征，不含标签
feature_cols = mats50[0].columns[:40].tolist()
X_raw = concat_patient_matrices(mats50, use_cols=feature_cols)
print(X_raw.shape)
4.6 Min-Max 归一化并保留缺失位置
def minmax_scale_with_nan(df):
X = df.values.astype(float)
X_scaled = X.copy()
mins = np.nanmin(X_scaled, axis=0)
maxs = np.nanmax(X_scaled, axis=0)
ranges = maxs – mins
ranges[ranges == 0] = 1.0
X_scaled = (X_scaled – mins) / ranges
return X_scaled, mins, ranges
X_scaled, mins, ranges = minmax_scale_with_nan(X_raw)
mask = ~np.isnan(X_scaled)
X_nmf, W, H, losses = masked_nmf(
X_scaled,
mask,
rank=6,
max_iter=600,
lr=0.003,
reg=1e-4
)
print("补全后是否仍有 NaN:", np.isnan(X_nmf).sum())
4.7 PCA 主成分分析
PCA 分析要使用 Z-score 标准化：
[
Z_{ij}=\\frac{x_{ij}-\\mu_j}{\\sigma_j}
]
协方差矩阵：
[
C=\\frac{1}{n-1}Z^TZ
]
特征分解：
[
Cv_k=\\lambda_k v_k
]
第 (k) 个主成分方差贡献率：
[
CR_k=\\frac{\\lambda_k}{\\sum_j\\lambda_j}
]
样本在第 (k) 个主成分上的得分：
[
Score_k=Zv_k
]
4.8 PCA 代码
from sklearn.preprocessing import StandardScaler
from sklearn.decomposition import PCA
def run_pca(X_completed, feature_names, n_components=2):
scaler = StandardScaler()
Z = scaler.fit_transform(X_completed)
pca = PCA(n_components=n_components)
scores = pca.fit_transform(Z)
loadings = pd.DataFrame(
pca.components_.T,
index=feature_names,
columns=[f"PC{i+1}" for i in range(n_components)]
)
explained = pca.explained_variance_ratio_
score_df = pd.DataFrame(
scores,
columns=[f"PC{i+1}_score" for i in range(n_components)]
)
return pca, score_df, loadings, explained
pca, score_df, loadings, explained = run_pca(X_nmf, feature_cols, n_components=2)
print("第一主成分贡献率:", explained[0])
print("第二主成分贡献率:", explained[1])
print("累计贡献率:", explained[:2].sum())
print("载荷最高的指标:")
print(loadings.abs().sort_values("PC1", ascending=False).head(10))
五、任务三：基于 PCA 结果的临床数据规律挖掘
5.1 主成分载荷的解释
PCA 的关键不是只算出方差贡献率，而是解释主成分代表什么临床含义。主成分载荷越大，说明该生命体征对该主成分贡献越大。
例如：
若 PC1 在 HR、Resp、Temp、WBC、Lactate 上载荷较高，可能代表炎症/感染应激状态；
若 PC2 在 SBP、MAP、DBP 上载荷较高，可能代表循环血流动力学状态；
若某主成分在 O2Sat、SaO2、PaCO2、FiO2 上载荷较高，可能代表呼吸氧合状态；
若 BUN、Creatinine 载荷较高，可能代表肾功能状态。
本题虽然背景提到肿瘤筛查，但数据标签中包含脓毒症标签，因此论文写作时可以更稳妥地表述为“患者高危生理状态预警指标”，不要强行说某指标一定是肿瘤特异指标。
5.2 关键预警指标筛选方法
可以用以下规则筛选关键指标：
在 PC1 或 PC2 中绝对载荷排名前 10；
在两类主成分中均具有较高载荷；
方差贡献率较高主成分中的高载荷变量优先；
缺失率不能过高；
具有明确临床意义。
综合评分可定义为：
[
S_j=CR_1|l_{j1}|+CR_2|l_{j2}|-\\alpha MR_j
]
其中 (CR_1,CR_2) 是第一、第二主成分贡献率，(l_{j1},l_{j2}) 是载荷系数，(MR_j) 是缺失率。该指标兼顾主成分贡献和数据可用性。
5.3 关键指标筛选代码
def select_key_indicators(loadings, explained, missing_rate, top_k=10, alpha=0.2):
mr = missing_rate.reindex(loadings.index).fillna(0)
score = (
explained[0] * loadings["PC1"].abs()
+ explained[1] * loadings["PC2"].abs()
– alpha * mr
)
result = pd.DataFrame({
"PC1_loading": loadings["PC1"],
"PC2_loading": loadings["PC2"],
"missing_rate": mr,
"importance_score": score
}).sort_values("importance_score", ascending=False)
return result.head(top_k), result
missing_rate_all = X_raw.isna().mean()
top_indicators, indicator_table = select_key_indicators(
loadings,
explained,
missing_rate_all,
top_k=12
)
print(top_indicators)
5.4 补全前后 PCA 对比
题目要求分别对 NMF 补全前残缺原始数据、补全后完整数据开展同等条件 PCA，对比主成分组成和方差贡献率变化。
由于 PCA 不能直接处理 NaN，因此补全前可以用简单均值填补作为“原始残缺数据的基础处理”基准，然后与 NMF 补全结果比较。
对比指标包括：
第一、第二主成分贡献率变化；
累计贡献率变化；
载荷排序稳定性；
关键指标是否更符合临床解释；
主成分得分散点是否更连续、更少异常点。
from sklearn.impute import SimpleImputer
from scipy.stats import spearmanr
# 补全前：简单均值填补作为基准
mean_imputer = SimpleImputer(strategy="mean")
X_mean = mean_imputer.fit_transform(X_scaled)
pca_mean, score_mean, load_mean, exp_mean = run_pca(X_mean, feature_cols, n_components=2)
pca_nmf, score_nmf, load_nmf, exp_nmf = run_pca(X_nmf, feature_cols, n_components=2)
print("均值基准 PCA 贡献率:", exp_mean)
print("NMF补全 PCA 贡献率:", exp_nmf)
# 载荷排序相关性
rank_corr_pc1 = spearmanr(load_mean["PC1"].abs(), load_nmf["PC1"].abs()).correlation
rank_corr_pc2 = spearmanr(load_mean["PC2"].abs(), load_nmf["PC2"].abs()).correlation
print("PC1载荷排序相关:", rank_corr_pc1)
print("PC2载荷排序相关:", rank_corr_pc2)
论文中可以这样解释：
如果 NMF 补全后第一、第二主成分贡献率更高，说明补全增强了数据中的主要结构；
如果载荷排序更稳定，说明关键指标识别更可靠；
如果主成分散点图中异常离散点减少，说明补全改善了数据质量。
六、任务四：多算法、多缺失率补全精度对比与最优模型筛选
6.1 缺失率模拟实验设计
题目要求设置 10%、20%、30%、40%、50% 五档缺失比例，在基准完整数据集上随机掩码生成缺失样本。
实验流程如下：
将 NMF 初步补全或较完整的数据作为基准完整数据；
在观测数据中随机选择一定比例位置设为 NaN；
用不同算法进行补全；
只在被人工遮盖的位置计算误差；
重复多次随机掩码，统计平均 RSE 和标准差。
6.2 RSE 指标
相对平方误差：
[
RSE=\\frac{\\sum_{(i,j)\\in\\Omega}(x_{ij}-\\hat{x}{ij})^2}{\\sum{(i,j)\\in\\Omega}x_{ij}^2}
]
也可以使用开方形式：
[
RSE=\\sqrt{\\frac{\\sum_{(i,j)\\in\\Omega}(x_{ij}-\\hat{x}{ij})^2}{\\sum{(i,j)\\in\\Omega}x_{ij}^2}}
]
论文中应明确使用哪一种。推荐使用带根号形式，便于和 RMSE 类似解释。
6.3 多算法补全方法
建议比较以下方法。
方法一：均值填补
优点：计算最快，基准方法。
缺点：低估方差，破坏变量相关性。
方法二：KNN 填补
优点：利用相似样本局部结构。
缺点：缺失率高时邻居不稳定，计算成本较高。
方法三：SVD-PCA 低秩补全
优点：利用全局低秩结构，适合强相关变量。
缺点：可能产生负值，对非负生命体征需后处理。
方法四：NMF 补全
优点：非负约束强，结果更符合生命体征非负特性；具有潜在生理模式解释。
缺点：对秩 r 和初始化敏感，迭代成本较高。
方法五：MICE 迭代回归填补
优点：变量间关系表达灵活。
缺点：高维高缺失时容易不稳定，计算较慢。
6.4 SVD-PCA 补全代码
def svd_impute(X_missing, rank=5, max_iter=100, tol=1e-5):
X = X_missing.copy()
mask = ~np.isnan(X_missing)
# 初始用列均值填补
col_mean = np.nanmean(X, axis=0)
inds = np.where(np.isnan(X))
X[inds] = np.take(col_mean, inds[1])
last_X = X.copy()
for it in range(max_iter):
U, S, Vt = np.linalg.svd(X, full_matrices=False)
S_rank = np.diag(S[:rank])
X_hat = U[:, :rank] @ S_rank @ Vt[:rank, :]
# 只更新缺失位置
X[~mask] = X_hat[~mask]
diff = np.linalg.norm(X – last_X) / (np.linalg.norm(last_X) + 1e-8)
if diff < tol:
break
last_X = X.copy()
# 非负裁剪
X = np.clip(X, 0, 1)
return X
6.5 多算法统一评估代码
from sklearn.impute import SimpleImputer, KNNImputer
from sklearn.experimental import enable_iterative_imputer
from sklearn.impute import IterativeImputer
def create_random_mask(X_full, missing_rate=0.2, random_state=2026):
rng = np.random.default_rng(random_state)
X_missing = X_full.copy()
n, p = X_full.shape
total = n * p
num_missing = int(total * missing_rate)
all_indices = np.arange(total)
selected = rng.choice(all_indices, size=num_missing, replace=False)
rows = selected // p
cols = selected % p
eval_mask = np.zeros_like(X_full, dtype=bool)
eval_mask[rows, cols] = True
X_missing[eval_mask] = np.nan
return X_missing, eval_mask
def rse_score(X_true, X_pred, eval_mask):
numerator = np.sum((X_true[eval_mask] – X_pred[eval_mask]) ** 2)
denominator = np.sum(X_true[eval_mask] ** 2) + 1e-8
return np.sqrt(numerator / denominator)
def impute_by_method(X_missing, method, rank=6):
if method == "mean":
return SimpleImputer(strategy="mean").fit_transform(X_missing)
if method == "knn":
return KNNImputer(n_neighbors=5).fit_transform(X_missing)
if method == "mice":
return IterativeImputer(max_iter=20, random_state=2026).fit_transform(X_missing)
if method == "svd":
return svd_impute(X_missing, rank=rank)
if method == "nmf":
mask = ~np.isnan(X_missing)
X_imp, _, _, _ = masked_nmf(
X_missing,
mask,
rank=rank,
max_iter=400,
lr=0.003,
reg=1e-4
)
return X_imp
raise ValueError("未知方法")
def evaluate_imputation_methods(X_full, missing_rates=(0.1,0.2,0.3,0.4,0.5),
methods=("mean","knn","mice","svd","nmf"),
repeats=5):
records = []
for rate in missing_rates:
for rep in range(repeats):
X_missing, eval_mask = create_random_mask(
X_full,
missing_rate=rate,
random_state=2026 + rep
)
for method in methods:
try:
X_pred = impute_by_method(X_missing, method, rank=6)
X_pred = np.clip(X_pred, 0, 1)
rse = rse_score(X_full, X_pred, eval_mask)
records.append({
"missing_rate": rate,
"repeat": rep,
"method": method,
"RSE": rse
})
except Exception as e:
records.append({
"missing_rate": rate,
"repeat": rep,
"method": method,
"RSE": np.nan
})
return pd.DataFrame(records)
# 以 X_nmf 作为近似完整基准数据进行模拟遮盖
result_df = evaluate_imputation_methods(X_nmf, repeats=3)
summary = result_df.groupby(["missing_rate", "method"])["RSE"].agg(["mean", "std"]).reset_index()
print(summary)
6.6 稳定性与复杂度评价
仅看 RSE 不够。最优模型要综合考虑：
平均 RSE；
RSE 标准差；
缺失率升高时误差增长速度；
计算时间；
临床解释性；
是否满足非负约束。
可定义综合得分：
[
Score_m=\\alpha \\overline{RSE}_m+\\beta SD(RSE_m)+\\gamma T_m+\\delta C_m
]
分数越低越好。其中 (T_m) 是归一化计算时间，(C_m) 是模型复杂度惩罚。
import time
def timed_impute(X_missing, method):
start = time.time()
X_pred = impute_by_method(X_missing, method)
elapsed = time.time() – start
return X_pred, elapsed
def model_selection_score(summary, time_table=None):
# 示例：仅用RSE均值和标准差
df = summary.copy()
df["score"] = df["mean"] + 0.5 * df["std"].fillna(0)
return df.sort_values(["missing_rate", "score"])
论文中可以给出结论：
低缺失率下，KNN 或 SVD-PCA 可能表现较好；
中高缺失率下，NMF 由于非负低秩约束更稳定；
极高缺失率下，所有算法误差都会明显增大，但 NMF 和 SVD-PCA 的退化更可控；
若考虑临床解释性，NMF 更适合作为 EHR 数据补全主模型。
七、论文写作结构建议
摘要
摘要应突出：EHR 高缺失问题、缺失率筛选、NMF 补全、PCA 关键指标挖掘、多算法 RSE 对比和最优模型推荐。
可写成：
本文针对电子健康记录数据高缺失、高维度和多量纲特征共存的问题，建立了基于缺失模式分析、非负矩阵分解补全和主成分分析的 EHR 数据修复与指标挖掘框架。首先对 training_setA 中患者 .psv 文件进行读取与缺失率统计，筛选缺失率小于 90% 的有效生命体征变量，并构建 50 份 36×41 标准化临床矩阵。随后采用 Min-Max 归一化与缺失掩码记录方法，构建带非负约束的 NMF 低秩补全模型，实现 EHR 缺失值填充。在补全数据基础上，利用 Z-score 标准化和 PCA 主成分分析提取主要生理状态模式，并基于主成分载荷和方差贡献率筛选关键预警指标。最后，在 10% 至 50% 五档人工缺失率下，比较均值填补、KNN、SVD-PCA、NMF 和 MICE 等方法的 RSE、稳定性与计算复杂度。实验表明，NMF 在中高缺失率下具有较优稳定性和临床解释性，适合作为本类 EHR 稀疏数据的核心补全模型。
第一章：问题重述与数据特点
说明 EHR 数据来源、变量类型、患者时序矩阵结构、缺失严重性和任务目标。
重点强调：
EHR 缺失不是普通随机缺失；
临床数据量纲差异大；
补全结果会影响后续 PCA 和风险指标筛选；
必须同时考虑补全精度和临床可解释性。
第二章：数据预处理与缺失分析
包括：
读取 p000001.psv；
统计前 40 个生命体征缺失率；
在前 36 个特征中筛选缺失率小于 90% 的前 11 个变量；
选取 50 份 36×41 标准矩阵；
归一化；
构造缺失掩码。
第三章：低维矩阵补全示例
对题目给出的小矩阵，采用均值填补、KNN 填补和迭代回归填补。说明不同方法特点，为后续高维补全做铺垫。
第四章：NMF 缺失补全模型
写清楚：
为什么 NMF 适合 EHR；
NMF 模型公式；
缺失掩码目标函数；
分解秩选择；
迭代优化方法；
补全结果生成完整数据集。
第五章：PCA 主成分分析
包括：
Z-score 标准化；
协方差矩阵；
特征值分解；
第一、第二主成分得分；
方差贡献率；
主成分载荷解释。
第六章：基于 PCA 的关键指标挖掘
根据载荷系数和贡献率筛选关键生命体征指标。
建议结合医学含义讨论 HR、O2Sat、Temp、SBP、MAP、Resp、Lactate、WBC、Creatinine 等变量可能对应的循环、呼吸、炎症、代谢、肾功能等状态。
第七章：补全前后 PCA 对比
对比均值填补基准和 NMF 补全后 PCA：
主成分贡献率是否提高；
载荷结构是否更稳定；
关键指标是否更集中；
补全是否提升数据质量。
第八章：多算法、多缺失率精度对比
设置 10%、20%、30%、40%、50% 缺失率，比较：
均值填补；
KNN；
MICE；
SVD-PCA；
NMF。
评价指标：
RSE；
RSE 标准差；
计算时间；
模型可解释性；
非负性约束。
第九章：结论与建议
总结：
高缺失变量不宜盲目补全；
NMF 能利用非负低秩结构，适合中高缺失 EHR；
PCA 可帮助筛选患者生理状态关键指标；
多缺失率模拟验证比单次补全更可信；
最优补全模型应综合精度、稳定性和临床解释性。
八、最终博客总结
赛题 3 的高分关键不是简单把 NaN 填上，而是建立一套完整的 EHR 缺失数据分析体系。首先，要从缺失机制出发，识别哪些变量值得补全、哪些变量应当剔除；其次，要根据 NMF 的非负低秩结构，构建适合生命体征数据的补全模型；再次，要通过 PCA 判断补全后的数据是否增强了主要生理结构；最后，要通过多算法、多缺失率的 RSE 实验验证模型的鲁棒性。
一句话概括：
本题赢在“缺失模式识别清楚、NMF 补全逻辑严谨、PCA 临床解释充分、多算法对比完整”。
如果要发布博客，可以在文末说明：
完整资源包含：论文全文、数据读取代码、缺失率统计结果、NMF 补全代码、PCA 主成分分析图表、多缺失率 RSE 对比表、算法稳定性分析和可复现实验说明。本文仅展示核心思路与部分代码。
