MEM1 数据读取与预处理

从数据表中提取四个关键变量：

• Distance (y)：垂体和翼上颌裂之间的距离（mm）
• Age (t)：测量时的年龄（固定在 8、10、12、14 岁）
• Individu：个体标识符，共 27 名受试者（其中 16 名男生，11 名女生）
• Sex：性别信息

数据集包含 108 条观测记录，每个个体最多 4 次测量。

MEM2 线性模型 LM1 的估计与拟合

在 LM1 模型中，所有个体被假设共享相同的截距与斜率：

基于TQ1的理论推导，最大似然估计得到的斜率和截距公式为：

残差方差的MLE估计

注意: Matlab 的 fitlm 默认使用无偏估计(如下公式)，和理论MLE估计结果有细微不同

计算贝叶斯信息准则：

表 MEM2-1 给出了 LM1 的参数估计结果。

表 MEM2-1：LM1 模型参数估计结果

图 MEM2-1 展示了 LM1 的拟合结果，蓝色点为实际观测值，红色直线为模型拟合线，其方程为 y=16.76+0.66×Age。

图 MEM2-1：LM1 全体线性回归拟合


LM1 拟合捕捉到了整体随年龄线性上升的趋势，但不能捕捉到性别差异和个体的变异性，数值表格中，LM1的残差方差相对较大，说明模型对数据的拟合效果并不特别理想

MEM3 - 分性别的LM1模型拟合分析

从图中可见，在各个年龄，男孩的上颌裂距离普遍大于女孩，且随着年龄增长，男孩观测值更加集中，可见过于简化的LM1模型无法充分捕捉数据中的性别差异。

MEM4 - LM2和LM3模型：考虑性别效应

LM2模型（共同截距，不同斜率）：

表 MEM4-1 给出了 LM2 的估计结果。

表 MEM4-1：LM2 模型参数估计结果

图 MEM4-1：LM2 模型拟合结果（共同截距，不同斜率）

图 MEM4-1 展示了 LM2 的拟合效果。蓝色为男孩数据与拟合曲线，粉色为女孩数据与拟合曲线。

LM2 的结果表明，男孩的斜率估计值（\hat{\beta}_{1b}=0.7477）高于女孩（\hat{\beta}_{1g}=0.5329），代表男孩的垂体与翼上颌裂间距的平均增长速度更快。

LM3模型（不同截距，不同斜率）：

表 MEM4-2：LM3 模型参数估计结果

图 MEM4-2：LM3 模型拟合结果

LM3 的结果表明，男孩的斜率（\hat{\beta}_{1b}=0.7844）显著高于女孩（\hat{\beta}_{1g}=0.4795），而女孩的截距（\hat{\beta}_{0g}=17.3727）则高于男孩（\hat{\beta}_{0b}=16.3406）。

MEM5 - 模型比较与个体拟合评估

为了系统评估 LM1、LM2 和 LM3 的拟合性能，我们比较了它们在残差平方和（RSS）、残差方差估计（σ²）、贝叶斯信息准则（BIC）以及均方根误差（RMSE）等指标上的表现。

表 MEM5-1：三个模型的综合比较

从整体拟合指标来看，LM2 和 LM3 相较 LM1 均表现显著更好，残差方差与 RSS 明显更小。但是尽管LM2的BIC更小，但是其方差和残差平方和比LM3略大。但是由于LM2比LM3使用更少的参数，从计算复杂度来说更好，所以我最喜欢LM2

特定个体的模型拟合比较

除了群体水平的指标，我们进一步考察了若干典型个体（Boy5、Boy7、Girl3、Girl4）的拟合效果。

图 MEM5-1：特定个体的模型拟合比较（LM1 vs. LM2 vs. LM3）

结果显示，尽管 LM2 和 LM3 在群体水平表现最佳，但在某些个体上，其拟合效果不如简单的LM1 qui 给出了更接近观测值的拟合。

这一现象表明，群体最优的固定效应模型并不能充分解释个体的差异。随着年龄增长，不同个体的变化趋势并不完全一致，因此在后续分析中引入 个体层面的随机效应，即混合效应模型（LMM）

MEM6 - 线性混合效应模型(LMM)的EM算法实现

LMM模型形式为：

其中 \beta_{0} 和 \beta_{1} 为固定效应，(\beta_{0r,i}, \beta_{1r,i}) 为个体特异性的随机效应，反映个体之间的异质性。

参数估计结果

采用 EM 算法对 LMM 进行估计，结果见表 MEM6-1。

表 MEM6-1：LMM 参数估计结果

LMM的固定效应估计值与LM1完全相同，这是因为LMM就是将LM1扩展为包含随机效应的模型

LMM 的固定效应估计值与 LM1 相同（\hat{\beta}_0=16.76, \hat{\beta}_1=0.66），这是因为 LMM 就是在 LM1 的基础上引入了随机效应。截距方差较大（\hat{\sigma}^2_{0r}= 4.8141）表明个体的起始值有很大的差异，而个体间的生长速度差异(斜率)很小(\hat{\sigma}^2_{1r}=0.0462)。截距与斜率之间存在显著负相关（\hat{\rho}_{01}=-0.5815），意味着起始值较高的个体往往生长速率较低，而起始值较低的个体则更可能具有较高的生长速率。

表 MEM6-2：LMM和三个传统线性模型比较

LMM 的残差方差显著低于 LM2（1.7162 vs. 4.9154），说明随机效应成功捕捉了大部分个体间的变异性，且LMM整体残差平方和大幅度减小，这表明考虑个体特异性的混合效应模型能够更准确的拟合观测

图 MEM6-1：LMM 与 LM2 在典型个体上的拟合比较

可以看到，在个体上，LMM 的拟合曲线显著更贴近观测点，相应的 RMSE 也明显小于 LM2。例如，Boy5 的 RMSE 在 LM2 下为 2.24，而在 LMM 下仅为 1.08。

第2部分：处理鲁棒性

R1 - 添加Student's t分布噪声

为每个个体i生成独立的自由度：\eta_i \sim \mathcal{U}[2.1, 3] ，较小的自由度产生厚尾分布，增加异常值概率

R2 - 鲁棒性分析结果

在添加 t 分布噪声后，首先对 LM1、LM2 和 LM3 模型重新进行了拟合。图 R2-1 至 R2-4 分别展示了带噪声数据下三类模型的拟合情况。

图 R2-1：LM1 带噪声数据的整体拟合

图 R2-1 展示了 LM1 在带噪声数据下的整体拟合情况。与原始数据相比，可以看到数据点的分布更加分散，线性回归被异常噪声拉扯

图 R2-2：LM1 按性别分组的拟合结果

图 R2-2 分别展示了 LM1 对男孩和女孩的拟合结果。分析同之前，男孩女孩的数据点的分布更加分散

图 R2-3：LM2 带噪声数据的拟合（共同截距，不同斜率）

图 R2-4：LM3 带噪声数据的拟合（不同截距，不同斜率）

表 R2-1：原始数据 vs. t 噪声数据的模型比较

从表中可见，所有模型在噪声数据下残差方差与 BIC 都显著上升，说明传统线性模型对异常值较为敏感。尽管如此，LM2 和 LM3 的性能仍优于 LM1，但整体拟合质量下降明显。

图 R2-5：特定个体的模型拟合比较（t噪声数据）

图 R2-5 给出了 Boy5、Boy7、Girl3 和 Girl4 的个体拟合比较。在这些个体中，t 噪声导致数据点出现强烈偏离。同样可见，即使LM2 LM3在群体上表现较好，但是在个体上可能不如LM1模型，因为个体的变异性没有考虑

t-LMM鲁棒混合效应模型

为了提高对异常值的鲁棒性，实现了基于Student's t分布的线性混合效应模型（t-LMM）。其思想是在随机效应与噪声项均采用 t 分布建模，以增强对厚尾和离群点的适应性。

表 R2-2：LMM 与 t-LMM 参数估计结果比较

图 R2-6：t-LMM 与 LM2 的个体拟合比较（t 噪声数据）

可见，t-LMM 在这些个体上的拟合均优于 LM2，即在存在离群点的情形下，拟合曲线更为稳健。这表明 t-LMM 在个体水平上显著提升了鲁棒性。

表 R2-3：t 噪声数据下的模型比较

可见，基于Student's t分布噪声的LMM模型有对噪声很好的鲁棒性，成功缓解了异常值的影响，在噪声环境下保持了更稳定的参数估计，有更少的残差、方差 和BIC，显著优于 LM1–LM3

R3 强异常值条件下的鲁棒性分析

为了测试模型在极端条件下的表现，向原始数据添加10%的强异常值，总共有11个异常数据点，数据范围从[16.50, 31.50]扩展至[7.86, 42.14]

异常值的添加逻辑不同会导致后续的结果显示不同，如果生成正负两个方向的异常值，且数量相当，即异常值分布相均匀，这会导致异常值的影响会政府正负抵消，这会影响我们对模型鲁棒性的分析，因此我们考虑80%异常值在平均观测之上上，20%在平均观测之下的情况，这样可以更明显的观察模型拟合线被异常值拖拽的情况。

图 R3-1：LM1 在异常值数据下的整体拟合

拟合直线明显被右上方的大异常值“拖拽”，整体曲线被抬升，少数下方异常点则造成拟合方差的扩大，模型对正常数据的解释力显著下降

图 R3-2：LM1 在异常值数据下的性别分组拟合

结论同上，无论男孩还是女孩样本，都由于高异常点导致拟合线明显抬高，模型对正常数据的解释力显著下降

图 R3-3：LM2 在异常值数据下的拟合

图 R3-4：LM3 在异常值数据下的拟合

LM2 LM3的结果图同理，被异常值向上拖拽，截距、斜率均被改变，出现了错误的推断

表R3-1：参数估计的偏差分析

同之前分析，传统线性模型参数被异常值大幅度改变，所有模型的 截距 β₀ 被大幅抬高，斜率 β₁ 被显著削弱(LM2 女孩斜率甚至降为 -0.0032，近乎为0)

图 R3-5：特定个体的模型拟合比较（异常值数据）

在个体图上，有和之前相同的结论，首先 总体表现好的 LM2 LM3 不如更简单的LM1，其次，所有LM1-LM3均受到异常值影响，代表传统线性模型（LM1-LM3）对异常值敏感

t-LMM的鲁棒性表现

为应对异常值，采用基于 Student’s t 分布的混合效应模型（t-LMM）。

表R3-2：异常值数据下的性能比较

可以看到，传统的 LM1–LM3 模型均出现了明显的性能下降，而 t-LMM 的表现显著优于其他模型：

图 R3-2：典型个体的 LM2 与 t-LMM 拟合对比（异常值数据）

可见，基于Student's t分布噪声的LMM模型有对噪声很好的鲁棒性，通过引入潜在权重对异常观测降权，成功缓解了异常值的影响，在噪声环境下保持了更稳定的参数估计，有更少的残差、方差 和BIC，拟合曲线更贴近主要观测趋势。

这里有问题

第3部分：处理缺失值

M1 - 缺失数据生成与可视化

模拟10% 完全随机缺失（MCAR）数据的情况，比较不同缺失处理策略对模型拟合的影响

图 M1-1：原始数据与缺失位置（10% MCAR）

在原始 108 条观测中，完全随机的缺失 10%数据（11 条记录），并用标记缺失位置

M2 - 缺失数据处理策略比较

两种处理策略:

• CC（Complete Case）：删除含缺失的数据行
• MI（Mean Imputation）：采用分组均值插补缺失值。

并将 完整观测（Full Observation, FO） 作为基准进行对照。

M2.1 完整案例（CC）策略

这里少个图

表 M2-1：[CC] LM1 模型参数估计结果

图 M2-2：[CC] LM1 性别分组拟合

图 M2-3：[CC] LM2 共同截距，不同斜率

图 M2-4：[CC] LM3 不同截距，不同斜率

表 M2-1: FO与CC策略下的模型比较

从表中可见，CC 策略下的残差方差普遍增大，说明模型拟合精度下降。然而 BIC 数值却降低，这是由于 BIC 的惩罚项与样本量 N 有关，N 减少导致对数似然项减小，从而影响了 BIC 的数值。

M2.2 均值插补（MI）策略

在 MI 策略下，利用 性别 × 年龄分组均值 填补缺失数据，恢复了完整的 108 条观测。结果如下：

图 M2-5：[MI] LM1 均值插补后拟合

图 M2-6：[MI] LM1 性别分组拟合

图 M2-7：[MI] LM2 共同截距，不同斜率

图 M2-8：[MI] LM3 不同截距，不同斜率

表 M2-2: FO与MI策略下的模型比较

与 FO 相比，MI 的残差方差与 BIC 基本接近，甚至略有改善。这表明在 10% MCAR 缺失 的情况下，均值插补能够较好地恢复数据结构，拟合效果可信度较高

M3 - 基于EM算法的LMM处理缺失数据

虽然我们使用CC与MI方法对10% MCAR 缺失数据进行了处理，但是这两种方法均存在局限，因此基于之前提出的线性混合效应模型（LMM），使用EM 算法来进行缺失数据的处理

基于条件期望塔性质，将缺失数据作为潜在变量，在 E 步补齐缺失并估计个体随机效应分布，M 步更新固定效应与方差参数。

参数估计结果

EM 收敛后得到的参数估计如表 M3-1 所示。

表 M3-1：EM算法的LMM参数估计

表 M3-2：不同方法的模型评估指标对比

M3-LMM 的残差平方和与残差方差均最小，相比于完整观测数据下的传统线性模型，或者CC/MI处理策略下的线性模型来讲，其模型拟合结果非常理想，证明了使用高斯算法来处理LMM方法在处理缺失数据方面有非常优秀的表现

相比于MEM6-LMM，我们发现两者结果基本相同，首先由于M3-LMM在缺失数据背景下运行，因此总数据量较MEM6-LMM 小，所以RSS和BIC都偏小，所以这并不意味着比MEM6-LMM更优秀，但是残差方差方面，其结果更好，说明说明算法EM在处理缺失数据时能够更好地捕捉数据的真实变异性

图 M3-1：个体上的模型拟合对比

从个体拟合比较图可以看出，M3-LMM 与 MEM6-LMM 在典型个体（Boy5、Boy7、Girl3、Girl4）上的拟合曲线几乎重合，两者的 RMSE 也非常接近。这说明在个体预测精度上，M3-LMM 的缺失数据处理并未削弱模型的解释力。相比之下，FO、CC、MI的拟合效果明显偏离观测点。

综上所述，M3-LMM 能够在缺失数据场景下保持与完整数据模型相近的拟合性能，验证了 EM 算法在补偿缺失信息和保持模型稳健性方面的有效性。