纵向数据缺失处理：FIML、TSRE与机器学习方法对比与选择指南

1. 项目概述与核心挑战

在纵向研究的实际操作中，数据缺失几乎是每个研究者都会遇到的“家常便饭”。想象一下，你正在追踪一群学生从入学到毕业的学业表现，计划每学期收集一次数据。但总会有学生因为转学、生病或单纯忘记填写问卷，导致某个时间点的数据“消失”了。这些缺失的数据点，就像拼图里丢失的几块，如果处理不当，最终描绘出的成长轨迹图可能完全失真，导致研究结论南辕北辙。这就是缺失数据处理的核心价值所在——它不仅仅是填补空白，更是为了确保我们从有限且不完整的信息中，能挖掘出最接近真相的规律。

传统上，我们依赖像全信息最大似然估计（FIML）这样的统计方法。它的思路很巧妙：既然数据缺失了，我就不去硬“猜”缺失值是什么，而是直接基于每个个体已有的数据（即“全信息”），通过最大化一个似然函数来估计模型参数。这好比侦探破案，不纠结于丢失的线索本身，而是利用所有现场留下的痕迹，去最合理地还原犯罪过程。FIML在数据符合多元正态分布且缺失机制为随机缺失（MAR）时，表现堪称完美。但现实数据往往“不听话”，存在偏态或异常值（即非正态），此时FIML的估计就可能产生偏差。

于是，两阶段稳健估计（TSRE）应运而生。它像一位经验丰富的质检员，分两步走：第一阶段，先用一种稳健的方法计算数据的均值和协方差矩阵，这个过程会自动给那些远离数据主体的“异常点”赋予较低的权重，从而得到一套“清洗”过的、更能代表主体数据的统计量；第二阶段，再将结构方程模型套用在这套稳健的统计量上进行估计。这种方法对非正态数据展现了良好的韧性。

与此同时，机器学习领域为我们提供了另一套工具箱。比如K近邻（KNN）插补，它的逻辑非常直观：要估计一个缺失值，就在数据集中找到和它最相似的k个“邻居”，然后用这些邻居的值的平均数（连续变量）或众数（分类变量）来填补。另一种强大的方法是随机森林（RF），特别是其单次插补实现missForest，它通过构建大量决策树来预测缺失值，能捕捉变量间复杂的非线性关系。更有甚者，为了考虑填补的不确定性，多重插补方法（如micecart和miceForest）被开发出来，它们不是只生成一个填补后的数据集，而是生成多个，最后将多个数据集的分析结果合并，得到更稳健的推断。

那么，面对纵向数据中的缺失，尤其是在棘手的非随机缺失（MNAR）和非正态分布并存的情况下，我们究竟该信任传统统计方法的深厚理论根基，还是拥抱机器学习方法的灵活与强大？这正是本次模拟研究试图回答的核心问题。我们将通过蒙特卡洛模拟，在可控的条件下，让FIML、TSRE、KNN、missForest、micecart和miceForest这六位“选手”同台竞技，看看在不同样本量、缺失率、缺失机制和数据分布的场景下，谁的表现更胜一筹。

2. 核心方法原理深度拆解

要理解不同方法的优劣，必须先深入其原理。这不仅关乎“怎么用”，更关乎“为什么这么用”以及“什么时候用最合适”。

2.1 传统统计方法：基于模型的优雅与局限

全信息最大似然估计（FIML）的精髓在于“因材施教”。对于纵向数据中的每个个体，无论他/她缺失了哪几个时间点的数据，FIML都会基于该个体实际观测到的数据，为其单独计算一个似然函数。这个似然函数衡量的是，在当前的模型参数下，观测到的数据出现的可能性有多大。最后，将所有个体的似然函数相加（或更准确地说，考虑到独立性，是相乘取对数），得到总似然函数，并通过迭代算法找到使这个总似然函数最大化的那组参数。这就是“最大似然”的含义。

注意：FIML并不直接“猜测”或填补缺失值。它是在模型估计过程中，将缺失机制（如果是MAR）整合进似然函数，从而在参数估计阶段就自动校正了缺失带来的影响。这避免了插补方法可能引入的额外误差。

然而，FIML的“阿喀琉斯之踵”在于其正态性假设。它的数学推导建立在误差项服从多元正态分布的基础上。当数据严重偏态或存在厚尾时，基于正态假设推导出的标准误和检验统计量（如卡方）就会不准，导致假设检验的I类错误率失控。在实际操作中，缺失的存在会让分布检验变得更困难，有时观测数据看起来正常，但总体其实是非正态的，这进一步增加了误用FIML的风险。

两阶段稳健估计（TSRE）的设计哲学是“去芜存菁”。它的第一阶段，采用稳健M估计量。你可以把它理解为一个智能过滤器：在计算均值和协方差时，它会根据数据点与主体数据的“距离”动态分配权重。距离中心越远的点（可能是异常值），权重越低；距离越近的点，权重越高。这样计算出的均值和协方差矩阵，对少数异常值不敏感，更能反映数据主体的特征。第二阶段，将标准的SEM拟合函数（如最大似然拟合函数）应用于这套稳健的均值和协方差矩阵，求解模型参数。

实操心得：TSRE的强大之处在于其第一阶段对非正态和异常值的稳健性。但在使用类似rsem这样的R包时，需要留意其默认的稳健估计方法（如Huber权重函数）以及调节参数。对于极端异常值，可能需要调整权重函数或考虑其他稳健方法。

2.2 机器学习方法：数据驱动的灵活与挑战

K近邻（KNN）插补的核心是“物以类聚”。它的成功极度依赖于两个因素：距离度量和k值选择。

• 距离度量：对于混合类型数据（连续变量和分类变量并存），常用高尔距离。它能同时处理不同类型的数据，计算一个综合的距离。如果所有变量都是连续的，欧氏距离就更简单直接。
• k值选择：这是KNN的“艺术”部分。k太小（如k=1），模型会对噪声极度敏感，相当于直接用最近邻的值填充，波动大；k太大，可能会把不相似的“远邻”也包含进来，平滑过度，掩盖了真实的局部结构。通常需要通过交叉验证来选择一个适中的k值。

随机森林（RF）与missForest：RF通过Bootstrap聚合（Bagging）和随机特征选择来构建大量差异化的决策树，然后通过投票（分类）或平均（回归）做出最终预测。missForest将RF用于插补的流程非常系统：

1. 初始化：用均值、中位数等简单方法给所有缺失值一个初始猜测。
2. 循环插补：按缺失比例从低到高排序变量。对每个有缺失的变量，将其作为“因变量”，其他所有变量（包括已用当前估计值填充的）作为“自变量”，训练一个随机森林模型来预测其缺失值。
3. 迭代更新：用新预测的值更新数据集，然后重复步骤2。
4. 收敛判断：比较本次迭代与上次迭代整个数据矩阵的差异（如均方误差）。当差异小于某个阈值或达到最大迭代次数时停止。

注意事项：missForest虽然通过迭代和大量树来提升精度，但其最终输出是单次插补数据集。这意味着它没有考虑插补本身的不确定性，可能导致后续分析的标准误被低估。这是所有单次插补方法的通病。

多重插补（MI）与mice框架：为了克服单次插补的缺陷，Rubin提出了多重插补理论。其核心思想是，既然缺失值不确定，我们就基于某种模型（如回归、决策树）随机抽取多个可能的值来填补，从而生成M个（通常M=20）完整的“副本”数据集。然后对每个数据集分别进行相同的模型分析，最后用Rubin法则将M个结果合并，得到总体参数估计、标准误和置信区间。mice包实现了这一思想，而micecart和miceForest则是用CART树或随机森林作为mice框架中每个变量的插补模型。

关键区别：missForest是单次插补，输出一个数据集；miceForest是多重插补，输出M个数据集，并通过Rubin法则合并结果。后者在理论上能提供更准确的区间估计。

3. 模拟研究设计与实操要点

为了公平、系统地比较这些方法，我们设计了一个基于增长曲线模型（GCM）的蒙特卡洛模拟。GCM是分析纵向数据的利器，它能同时刻画个体的发展轨迹（斜率）和起点差异（截距）。

3.1 数据生成机制

我们模拟一个最简单的无条件线性增长模型，有4个时间点（T=4）。

• 总体模型参数：我们设定潜变量截距的均值为6，斜率的均值为2。这意味着平均起始水平是6，平均每单位时间增长2。截距和斜率的方差都设为1，两者相关系数为0（独立）。测量误差的方差也设为1。
• 数据分布（关键操纵）：我们通过操纵测量误差ei的分布来生成不同形态的数据：
  1. 正态分布：ei ~ N(0,1)，作为基准。
  2. 对数正态分布：ei ~ Log-Normal(0,1)，产生严重的右偏态数据。
  3. t分布（5个自由度）：ei ~ t(5)，产生厚尾分布，极端值出现概率更高。
  4. 正态分布含5%异常值：95%的数据来自N(0,1)，5%的数据来自N(5,1)，模拟数据被少量高离群值污染的情况。
• 缺失机制生成：
  • MAR：让时间点t的观测值是否缺失，依赖于前一个时间点t-1的观测值。例如，如果yi1（第一次测量）的值超过某个阈值c1，则yi2, yi3, yi4缺失。阈值c1根据设定的缺失率（如15%）反推确定。
  • MNAR：生成一个与潜变量斜率bis高度相关（r=0.8）的辅助变量Auxi。如果Auxi超过某个百分位数阈值pt，则对应时间点的yit缺失。由于Auxi不进入增长模型分析，且与未观测的斜率相关，因此缺失机制是非随机的。

3.2 实验因子与软件实现

我们操纵了四个可能影响方法表现的因素，构成一个全面的实验设计：

1. 样本量（N）：100， 200， 500， 1000， 5000。覆盖了小样本到大数据场景。
2. 缺失机制：MAR， MNAR。
3. 缺失率（mr）：0%（基线）， 5%， 15%， 30%。
4. 数据分布：上述四种。

对于每个由（N， 机制， mr， 分布）唯一确定的实验条件，我们使用R语言生成500个模拟数据集。然后，分别用六种方法处理每个数据集：

• FIML：使用lavaan包，在模型设定中直接使用missing = "fiml"选项。
• TSRE：使用rsem包进行稳健估计。
• KNN：使用VIM包中的kNN()函数。关键参数k（邻居数）需要通过预实验优化。
• missForest：使用missForest包中的同名函数。
• micecart：使用mice包，设置method = "cart"。
• miceForest：使用mice包，设置method = "rf"。

对于多重插补方法（micecart， miceForest），我们将插补次数m设为20，这是文献中公认的能达到良好效率与精度平衡的常用值。

3.3 超参数调优实战

机器学习方法的性能对超参数敏感。我们在正式模拟前，针对KNN、CART和RF进行了小规模的预实验来调优。

• 对于KNN：我们尝试了k = 3, 5, 7, 10, sqrt(N)等不同取值。在纵向数据中，由于时间序列的自相关性，k值不宜过大，否则会模糊时间趋势。预实验发现，对于我们的模拟设置，k=5在大多数情况下能在偏差和方差间取得较好平衡。
• 对于CART（micecart）和RF（missForest/miceForest）：主要超参数是树的数量（ntree）和树的最大深度/节点复杂度。我们测试了ntree = 100, 200, 500。结果发现，超过200棵树后，性能提升微乎其微，但计算成本显著增加。因此，在正式模拟中，我们设定ntree = 200。对于树复杂度，我们使用包的默认设置，这通常包含了防止过拟合的剪枝策略。

调优经验：超参数调优没有“一招鲜”。对于实际研究数据，建议在分析前，划出一小部分完整数据，人为制造已知的缺失，然后用不同的超参数设置进行插补，以均方误差（MSE）或平均绝对误差（MAE）作为指标，选择在你自己数据上表现最好的那组参数。这是一个虽繁琐但至关重要的步骤。

4. 性能评估指标与结果解读

我们如何判断一个方法“好”还是“不好”？在模拟研究中，我们关注两个核心维度：准确性和效率。

4.1 评估指标详解

1. 相对偏差（Relative Bias, RB）：衡量参数估计值的平均值与真实参数值之间的系统性偏离程度。RB = (平均估计值 - 真实参数) / 真实参数RB的绝对值越小越好。通常，|RB| < 0.05被认为是可以接受的轻微偏差，|RB| > 0.10则表明存在值得关注的偏差。

2. 均方误差（Mean Squared Error, MSE）：衡量估计值的总体误差，它同时包含了偏差和方差。MSE = 平均[ (估计值 - 真实参数)^2 ] = 方差(估计值) + (偏差)^2MSE越小，说明方法的估计越精确、越稳定。

3. 经验覆盖率（Empirical Coverage Rate）：对于95%的置信区间，我们计算在500次重复模拟中，有多少次置信区间包含了真实的参数值。理想覆盖率是95%。覆盖率显著低于95%说明区间过窄（过于自信），高于95%说明区间过宽（过于保守）。

4.2 主要研究发现与场景化建议

基于模拟结果，我们可以得出一些具有高度实践指导意义的结论：

场景一：数据为MNAR机制（最棘手的情况）

• 发现：FIML在大多数MNAR条件下表现出了最小的相对偏差和最理想的覆盖率，即使是在非正态数据下（如t分布），其偏差也控制在可接受范围内。
• 原理分析：这似乎有悖于FIML的正态假设弱点。但关键在于，在MNAR下，缺失机制与未观测值相关。FIML作为一种基于模型的方法，其似然函数在特定模型设定下（虽然我们模拟的是无条件模型），某种程度上内在地对缺失机制进行了建模（尽管是默认的）。而大多数插补方法（无论是KNN还是随机森林）在MNAR下，仅基于观测数据来预测缺失值，无法校正由缺失机制本身引起的偏差。
• 实操建议：当你有强烈理由怀疑数据是MNAR时（例如，在追踪研究中，因病情恶化而失访），应优先考虑FIML。即使数据非正态，其带来的偏差可能也小于MNAR机制未被校正带来的偏差。TSRE在MNAR下的表现不稳定，有时偏差较大。

场景二：数据为MAR机制，且呈非正态分布

• 发现：TSRE是此场景下的“冠军”。它在处理含异常值或厚尾分布的数据时，偏差和MSE均显著低于FIML。
• 原理分析：MAR机制下，缺失只与观测变量有关，问题相对“单纯”。此时，FIML受非正态性的影响被放大，而TSRE第一阶段稳健估计的优势得以充分发挥，有效抵抗了异常值和偏态的影响。
• 实操建议：对于社会行为科学中常见的偏态数据（如反应时、问卷得分），若缺失机制可近似为MAR，TSRE是比FIML更可靠的选择。在R中，rsem包提供了便捷的TSRE实现。

场景三：大样本（N ≥ 1000）、低缺失率、极端偏态分布

• 发现：missForest在这种特定组合下展现了优势，其估计精度有时能与TSRE媲美甚至略优。
• 原理分析：随机森林需要足够的数据来学习复杂的模式。大样本提供了丰富的“学习材料”，低缺失率意味着模型构建时有更多可用的观测信息，而极端偏态分布下，树模型捕捉非线性关系的能力可能优于基于线性假设的稳健方法。但一旦样本变小或缺失率升高，missForest的性能会迅速下降。
• 实操建议：仅在拥有海量数据且缺失很少，同时数据分布极其复杂、非线性关系强烈时，可以考虑尝试missForest。对于大多数中小样本的纵向研究，它并非首选。

场景四：多重插补 vs. 单次插补

• 发现：在本次模拟的纵向增长模型背景下，micecart和miceForest的表现并未显著优于其对应的单次插补思想或传统方法。有时，由于插补模型与后续分析模型（GCM）不匹配，甚至会产生更大的误差。
• 原理分析：多重插补的理论优势在于考虑不确定性，但这建立在“插补模型正确或至少是兼容的”前提下。在纵向数据中，时间序列的自相关结构和增长轨迹是复杂的关系。用默认的、变量间相互预测的mice框架，可能无法很好地保持这种时间动态结构，导致插补值扭曲了真实的个体变化轨迹。
• 实操建议：对于纵向数据，不要盲目使用通用的多重插补包。如果要用，必须使用专门为纵向数据设计的插补方法（如jomo,pan包），或者确保你的插补模型明确包含了时间结构和随机效应。否则，简单的FIML或TSRE可能更安全。

4.3 各方法优缺点速查与选择指南

5. 实战案例与避坑指南

让我们用一个简化例子，演示如何在R中实际应用FIML和TSRE分析一个含有缺失的纵向数据集。假设我们有一个名为long_data的数据框，包含ID， time， score列，格式为长格式。

# 加载必要的包 library(lavaan) # 用于FIML library(rsem) # 用于TSRE library(ggplot2) # 用于绘图 # 假设数据已加载为 long_data # 1. 将数据转换为宽格式，以便SEM软件使用 wide_data <- reshape(long_data, timevar = "time", idvar = "ID", direction = "wide") # 2. 使用FIML拟合线性增长曲线模型 model_fiml <- ' # 定义潜变量 i =~ 1*score.1 + 1*score.2 + 1*score.3 + 1*score.4 # 截距因子 s =~ 0*score.1 + 1*score.2 + 2*score.3 + 3*score.4 # 斜率因子（线性增长） # 估计截距和斜率的均值与方差 i ~~ i s ~~ s i ~~ s # 测量误差方差（可设定为相等） score.1 ~~ res*score.1 score.2 ~~ res*score.2 score.3 ~~ res*score.3 score.4 ~~ res*score.4 ' fit_fiml <- growth(model_fiml, data = wide_data, missing = "fiml") summary(fit_fiml, fit.measures = TRUE) # 3. 使用TSRE拟合相同模型 # rsem通常需要直接指定变量名和模型 # 注意：rsem的语法可能与lavaan略有不同，请参考其文档 # 这里仅为示意流程 # fit_tsre <- rsem(model_fiml, data = wide_data) # summary(fit_tsre) # 4. 模型比较与诊断 # 检查FIML的正态性假设（在缺失数据下需谨慎） lavInspect(fit_fiml, "sampstat") # 获取样本矩 # 可以尝试对观测数据部分进行正态性检验（如MVN包），但解释力有限。 # 绘制个体增长轨迹（仅使用观测数据） ggplot(long_data, aes(x = time, y = score, group = ID)) + geom_line(alpha = 0.3) + stat_smooth(aes(group = 1), method = "lm", se = TRUE, color = "red", size = 1.2) + labs(title = "个体增长轨迹与总体趋势（含缺失）", x = "测量时间点", y = "得分")

常见问题与排查技巧实录

1. 问题：使用FIML时，模型无法收敛或出现“非正定矩阵”错误。

   • 可能原因：样本量过小、缺失率过高导致信息不足；模型设定错误（如过于复杂）；数据存在极端共线性。
   • 排查：
     • 检查样本量，尝试用更简单的模型。
     • 检查变量间的相关性矩阵，看是否有|r| > 0.9的变量。
     • 尝试用listwise删除法先分析完整案例，看模型是否设定有误。
     • 对于增长模型，确保时间编码正确，因子载荷矩阵设定无误。

2. 问题：TSRE估计结果与FIML差异巨大。

   • 可能原因：数据中存在大量或极强的异常值，TSRE的稳健估计将其权重降得很低。
   • 排查：
     • 绘制数据的箱线图或散点图，直观检查异常值。
     • 分别用FIML和TSRE拟合后，比较个案的影响统计量（如Cook‘s D），看哪些个案被TSRE严重降权。
     • 考虑异常值的实质意义：是数据录入错误，还是真实的极端个体？根据研究目的决定是保留（用TSRE）还是剔除（需谨慎并报告）。

3. 问题：使用机器学习插补（如missForest）后，后续的GCM分析结果变得不合理（如斜率方差为负）。

   • 可能原因：插补过程破坏了变量间的协方差结构，特别是时间序列的自相关。单次插补忽略了不确定性，导致“确定性”过强。
   • 排查：
     • 比较插补前后数据的相关矩阵和描述统计量。
     • 强烈建议：对于纵向分析，优先使用直接处理缺失的方法（FIML， TSRE）或专门为纵向设计的多重插补。如果必须使用通用插补，考虑在插补模型中纳入时间变量作为预测因子，并尝试使用能保持时间结构的算法。

4. 问题：如何判断数据是MAR还是MNAR？

   • 这是一个根本性难题，无法仅从数据本身得到确证。
   • 实践策略：
     • 理论驱动：基于研究设计和领域知识进行判断。例如，在药物临床试验中，因副作用退出更可能是MNAR；在邮寄调查中，忘记填写更可能是MCAR/MAR。
     • 敏感性分析：这是处理MNAR怀疑的黄金标准。使用不同的方法（如FIML， 选择模型， 模式混合模型）分析数据，看结论是否一致。如果结论稳健，则对缺失机制的假设依赖较小。
     • 收集辅助信息：尽可能收集与缺失可能相关的变量（即使它不是你研究的核心变量），并将其纳入模型，尝试将MNAR转化为MAR。

在我多年的数据分析经历中，最深的一点体会是：没有一种处理缺失数据的方法是“银弹”。FIML和TSRE等模型方法与你的分析模型紧密结合，理论性质明确，是我在严肃的纵向研究中的首选。机器学习插补方法更像一把灵活的“瑞士军刀”，在数据探索、预处理或预测任务中非常有用，但在用于复杂的纵向模型推断时，需要格外警惕其与后续分析模型的兼容性。最终的选择，必须基于对数据缺失机制的理解、数据分布的认识、样本量的大小，以及最重要的——通过敏感性分析来检验你的结论是否经得起不同处理方式的考验。当FIML和TSRE给出相似结论时，我对结果就更有信心了；当它们分歧很大时，这正是需要深入挖掘、并在论文中坦诚报告的时候。缺失数据处理，终究是一场与不确定性的共舞，我们的目标是让舞步尽可能稳健，而非追求绝对完美。