基于非规则数据的模型统计学习 —— 应用于正畸数据分析
背景介绍
主要针对正畸数据建模，学习参数混合效应模型（MEM），并分析鲁棒性，以应对可能存在的异常值和/或缺失值的情况。
数据来源于 16 名男孩和 11 名女孩的正畸研究，年龄范围为 8 至 14 岁。响应变量是垂体与翼腭裂之间的距离（单位为毫米），在 8、10、12 和 14 岁分别测量了每个男孩和女孩的距离。
一开始，我们尝试用简单的线性回归来描述每个孩子的变化趋势。但在实际操作中发现，有些孩子的变化速度（斜率）明显比其他孩子快，数据的波动更大（尤其男孩），拟合出来的线不足以描述这种波动，因此我们意识到除了整体的共同趋势外，还有个体差异性。由此我们采用混合效应模型（MEM），它不仅能捕捉整体的变化趋势，还能反映个体差异，非常适用于处理异常值或者缺失值的非规则数据。
传统方法
我们先介绍几种传统方法，也就是最基本的线性回归模型。
首先令 y_{ij} \in \mathbb{R}^{I \times J} 表示第 i 个个体在第 j 个年龄 t_j 时的正畸测量值。我们假设正畸测量值和年龄之间的关系是线性的。
TQ1 模型 LM1
LM1 最简单，它假设所有个体在正畸测量和年龄之间有相同的线性关系（斜率和截距相同），其模型表达式为：
y_{ij}=\beta_0+\beta_1t_j+\epsilon_{ij},\quad \epsilon_{ij}\sim\mathcal{N}(0,\sigma^2)
- y_{ij} \in \mathbb{R}^{I \times J} 表示第 i 个个体在第 j 个年龄 t_j 时的正畸测量值。比如说第一个哥们在12岁时候的正畸测量值。
- \beta_0 就是截距，它代表了当年龄 t_j 为 0 时（刚出生的时候），正畸测量值是多少。实际上，如果我们对年龄进行了中心化处理， \beta_0 就表示在参考时间点（比如平均年龄）时的测量值。
- \beta_1 是斜率，表示每长大一岁，正畸测量值增加（减少）多少
- t_j 是年龄变量
- \epsilon_{ij} 是模型中的随机误差项（残差），标准操作，将确定性模型变成概率模型，从而才可以进行统计推断方法
- 我们现在手上一共有 N = I \times J 个测量数据
现在我们要利用这个模型找到最合适的 \beta_0 和 \beta_1 的取值，来让模型更好的解释现在已有的数据，这就是最大似然估计的思想，因此我们假设了误差项 \epsilon_{ij} 服从正态分布，然后准备采用最大似然估计方法求出参数。
由于我们已知在正态误差假定下，最大似然估计与最小二乘估计是等价的，我们实际上就把问题变成最小化平方误差和
S=\sum_{i=1}^{I}\sum_{j=1}^{J}(y_{ij}-\beta_0-\beta_1t_j)^2
证明: 在正态误差假定下最大似然估计和最小二乘估计等价
根据正态误差假定，y_{ij} 在给定参数 \beta_0,\beta_1 时服从正态分布，均值为 \beta_0+\beta_1t_j ，方差为 \sigma^2。因此，每个观测值的概率密度为
p(y_{ij}\mid\beta_0,\beta_1,t_j)
=\frac{1}{\sqrt{2\pi\sigma^2}}
\exp\!\left(-\frac{(y_{ij}-\beta_0-\beta_1t_j)^2}{2\sigma^2}\right)
由于各个观测值独立（视作独立），则所有数据的联合似然函数为
L(\beta_0,\beta_1,\sigma^2)
=\prod_{i=1}^{I}\prod_{j=1}^{J}\frac{1}{\sqrt{2\pi\sigma^2}}
\exp\!\left(-\frac{(y_{ij}-\beta_0-\beta_1t_j)^2}{2\sigma^2}\right)
对似然函数取对数 ——> 对数似然函数
\begin{aligned}
\log L(\beta_0,\beta_1,\sigma^2)
&=\sum_{i=1}^{I}\sum_{j=1}^{J}\left[-\frac{1}{2}\log(2\pi\sigma^2)
-\frac{(y_{ij}-\beta_0-\beta_1t_j)^2}{2\sigma^2}\right] \\
&= -\frac{N}{2}\log(2\pi\sigma^2)
-\frac{1}{2\sigma^2}\sum_{i=1}^{I}\sum_{j=1}^{J}(y_{ij}-\beta_0-\beta_1t_j)^2
\end{aligned}
注意到第一项-\frac{N}{2}\log(2\pi\sigma^2)与\beta_0,\beta_1无关，因此在最大化\log L时，我们只需使第二项 达到最大值
-\frac{1}{2\sigma^2}\sum_{i=1}^{I}\sum_{j=1}^{J}(y_{ij}-\beta_0-\beta_1t_j)^2
因此等价于使下式达到最小值，等价于最小化平方误差和的最小二乘估计方法
S=\sum_{i=1}^{I}\sum_{j=1}^{J}(y_{ij}-\beta_0-\beta_1t_j)^2
证毕。
所以现在我们的目标是最小化平方误差和
S=\sum_{i=1}^{I}\sum_{j=1}^{J}(y_{ij}-\beta_0-\beta_1t_j)^2
- 对于每一个数据点，我们计算预测值
\beta_0+\beta_1t_j
与实际观测值
y_{ij}
之间的差异，这个差异就是
\epsilon_{ij}
- 平方的原因在于防止正负误差互相抵消，并更重视大误差
- 所有平方误差加起来，形成整体指标 S
- 调整参数
\beta_0
和
\beta_1
使得 $S$ 达到最小值，找到使得所有误差平方和最小的那个参数组合
为求解 \beta_0 和 \beta_1，我们对 S 分别对这两个参数求偏导并令其为零
对 \beta_0 求偏导为零：
\frac{\partial S}{\partial \beta_0} = -2 \sum_{i=1}^{I} \sum_{j=1}^{J} \Bigl( y_{ij} - \beta_0 - \beta_1 t_j \Bigr) = 0
整理后得到
\sum_{i=1}^{I} \sum_{j=1}^{J} y_{ij} = N \beta_0 + \beta_1 \sum_{i=1}^{I} \sum_{j=1}^{J} t_j \tag{1}
对 \beta_1 求偏导为零：
\frac{\partial S}{\partial \beta_1} = -2 \sum_{i=1}^{I} \sum_{j=1}^{J} t_j \Bigl( y_{ij} - \beta_0 - \beta_1 t_j \Bigr) = 0
整理后得到
\sum_{i=1}^{I} \sum_{j=1}^{J} t_j y_{ij} = \beta_0 \sum_{i=1}^{I} \sum_{j=1}^{J} t_j + \beta_1 \sum_{i=1}^{I} \sum_{j=1}^{J} t_j^2 \tag{2}
为了简化表达，我们定义总体均值
\bar{y} = \frac{1}{N} \sum_{i=1}^{I} \sum_{j=1}^{J} y_{ij}, \quad \bar{t} = \frac{1}{N} \sum_{i=1}^{I} \sum_{j=1}^{J} t_j
因此，公式（1）可写为:
N\,\bar{y} = N\,\beta_0 + \beta_1 \,N\,\bar{t}
两边同时除以 $N$，得到
\bar{y} = \beta_0 + \beta_1\,\bar{t}
因此截距的估计表达式为
\hat{\beta}_0^{LM1} = \bar{y} - \hat{\beta}_1^{LM1} \bar{t}
原公式（2）为
\sum_{i,j} t_j y_{ij} = \beta_0\sum_{i,j} t_j + \beta_1\sum_{i,j} t_j^2
将 \hat{\beta}_0^{LM1} 代入（2）可得
\sum_{i,j} t_j y_{ij} = \Bigl(\bar{y}-\beta_1\bar{t}\Bigr)\,N\,\bar{t} + \beta_1\sum_{i,j} t_j^2
展开
\sum_{i,j} t_j y_{ij} = N\,\bar{y}\,\bar{t} - \beta_1 N\,\bar{t}^2 + \beta_1\sum_{i,j} t_j^2
把含有
\beta_1
的项移到一边，得到
\sum_{i,j} t_j y_{ij} - N\,\bar{y}\,\bar{t} = \beta_1\Bigl[\sum_{i,j} t_j^2 - N\,\bar{t}^2\Bigr]
解得斜率的估计为
\hat{\beta}_1^{LM1} = \frac{\sum_{i,j} t_j y_{ij} - N\,\bar{y}\,\bar{t}}{\sum_{i,j} t_j^2 - N\,\bar{t}^2}
分子和分母可以写成中心化的形式
\sum_{i,j}(t_j-\bar{t})(y_{ij}-\bar{y}) = \sum_{i,j} t_j y_{ij} - N\,\bar{y}\,\bar{t}
\sum_{i,j}(t_j-\bar{t})^2 = \sum_{i,j} t_j^2 - N\,\bar{t}^2
证明: 中心化形式等价
\begin{aligned}
\sum_{i,j}(t_j-\bar{t})(y_{ij}-\bar{y})
&=\sum_{i,j}\Bigl[t_jy_{ij} - t_j\bar{y} - \bar{t}y_{ij} + \bar{t}\bar{y}\Bigr] \\
&=\sum_{i,j}t_jy_{ij} - \bar{y}\sum_{i,j}t_j - \bar{t}\sum_{i,j}y_{ij} + \sum_{i,j}\bar{t}\bar{y}
\end{aligned}
由于 \bar{y} 和 \bar{t} 是常数，不随 i,j 变化，所以 \sum_{i,j}\bar{t}\bar{y} = N\,\bar{t}\bar{y} 。 与此同时，我们已知 \sum_{i,j}t_j = N\,\bar{t},\quad \sum_{i,j}y_{ij} = N\,\bar{y} ，我们把上述结果代入上式
\begin{aligned}
\sum_{i,j}(t_j-\bar{t})(y_{ij}-\bar{y})
&=\sum_{i,j}t_jy_{ij} - \bar{y}(N\,\bar{t}) - \bar{t}(N\,\bar{y}) + N\,\bar{t}\bar{y} \\
&=\sum_{i,j}t_jy_{ij} - N\,\bar{t}\bar{y} - N\,\bar{t}\bar{y} + N\,\bar{t}\bar{y} \\
&=\sum_{i,j}t_jy_{ij} - N\,\bar{t}\bar{y}
\end{aligned}
第一个公式证毕
\begin{aligned}
\sum_{i,j}(t_j-\bar{t})^2
&=\sum_{i,j}\Bigl(t_j^2 - 2t_j\bar{t} + \bar{t}^2\Bigr)\\[1mm]
&=\sum_{i,j}t_j^2 - 2\bar{t}\sum_{i,j}t_j + \sum_{i,j}\bar{t}^2
\end{aligned}
我们依旧有
\sum_{i,j}t_j = N\,\bar{t}\quad \text{和}\quad \sum_{i,j}\bar{t}^2 = N\,\bar{t}^2
因此
\begin{aligned}
\sum_{i,j}(t_j-\bar{t})^2
&=\sum_{i,j}t_j^2 - 2\bar{t}(N\,\bar{t}) + N\,\bar{t}^2\\[1mm]
&=\sum_{i,j}t_j^2 - 2N\,\bar{t}^2 + N\,\bar{t}^2\\[1mm]
&=\sum_{i,j}t_j^2 - N\,\bar{t}^2
\end{aligned}
第二个公式证毕
中心化就是原始项减去均值，目的是去掉常数项，专注于数据相对于均值的变化情况。
因此，斜率估计可简写为
\hat{\beta}_1^{LM1}=\frac{\sum_{i,j}(t_j-\bar{t})(y_{ij}-\bar{y})}{\sum_{i,j}(t_j-\bar{t})^2}
综上所述，LM1下的最大似然估计（MLE）为
\begin{cases}
\hat{\beta}_0^{LM1} = \bar{y} - \hat{\beta}_1^{LM1}\bar{t}\\
\hat{\beta}_1^{LM1} = \frac{\sum_{i,j}(t_j-\bar{t})(y_{ij}-\bar{y})}{\sum_{i,j}(t_j-\bar{t})^2}
\end{cases}
TQ2 模型 LM2：
LM2 考虑了性别对变化速率的影响（之前提到男生变化速度更快），因此我们假定所有个体共享截距（ \beta_0 ），但是斜率不同（由两个参数描述，分别对应男孩女孩）
设定模型为
y_{ij}=\beta_0+\beta_{1b}t_j\mathbb{I}_{\text{sex=boy}}+\beta_{1g}t_j\mathbb{I}_{\text{sex=girl}}+\epsilon_{ij},\quad \epsilon_{ij}\sim\mathcal{N}(0,\sigma^2)
对于男孩（\mathbb{I}_{\text{sex=boy}} = 1，\mathbb{I}_{\text{sex=girl}} = 0），模型为
y_{ij} = \beta_0 + \beta_{1b} t_j + \epsilon_{ij}
而对于女孩（\mathbb{I}_{\text{sex=girl}} = 1，\mathbb{I}_{\text{sex=boy}} = 0），模型为
y_{ij} = \beta_0 + \beta_{1g} t_j + \epsilon_{ij}
我们构造平方误差和为两部分之和
S = \sum_{\text{boy}} \left( y_{ij} - \beta_0 - \beta_{1b} t_j \right)^2 + \sum_{\text{girl}} \left( y_{ij} - \beta_0 - \beta_{1g} t_j \right)^2
接下来，对每个参数求偏导并令其为零。
对 \beta_0 求偏导：
\frac{\partial S}{\partial \beta_0} = -2 \sum_{\text{boy}} \left( y_{ij} - \beta_0 - \beta_{1b} t_j \right) - 2 \sum_{\text{girl}} \left( y_{ij} - \beta_0 - \beta_{1g} t_j \right) = 0
整理可得
N \beta_0 + \beta_{1b} \sum_{\text{boy}} t_j + \beta_{1g} \sum_{\text{girl}} t_j = \sum_{i,j} y_{ij}
因此，\beta_0 的估计为
\hat{\beta}_0^{LM2} = \frac{1}{N} \left[ \sum_{i,j} y_{ij} - \hat{\beta}_{1b}^{LM2} \sum_{\text{boy}} t_j - \hat{\beta}_{1g}^{LM2} \sum_{\text{girl}} t_j \right]
对男孩部分的 \beta_{1b} 求偏导（仅考虑男孩数据）:
\frac{\partial S}{\partial \beta_{1b}} = -2 \sum_{\text{boy}} t_j \left( y_{ij} - \beta_0 - \beta_{1b} t_j \right) = 0
解得
\hat{\beta}_{1b}^{LM2} = \frac{\sum_{\text{boy}} t_j \left( y_{ij} - \hat{\beta}_0^{LM2} \right)}{\sum_{\text{boy}} t_j^2}
对女孩部分的 \beta_{1g} 求偏导（仅考虑女孩数据）：
\frac{\partial S}{\partial \beta_{1g}} = -2 \sum_{\text{girl}} t_j \left( y_{ij} - \beta_0 - \beta_{1g} t_j \right) = 0
得到
\hat{\beta}_{1g}^{LM2} = \frac{\sum_{\text{girl}} t_j \left( y_{ij} - \hat{\beta}_0^{LM2} \right)}{\sum_{\text{girl}} t_j^2}
因此，LM2 模型下的最大似然估计（MLE）：
\begin{cases}
\hat{\beta}_0^{LM2} = \frac{1}{N} \left[ \sum_{i,j} y_{ij} - \hat{\beta}_{1b}^{LM2} \sum_{\text{boy}} t_j - \hat{\beta}_{1g}^{LM2} \sum_{\text{girl}} t_j \right] \\
\hat{\beta}_{1b}^{LM2} = \frac{\sum_{\text{boy}} t_j \left( y_{ij} - \hat{\beta}_0^{LM2} \right)}{\sum_{\text{boy}} t_j^2} \\
\hat{\beta}_{1g}^{LM2} = \frac{\sum_{\text{girl}} t_j \left( y_{ij} - \hat{\beta}_0^{LM2} \right)}{\sum_{\text{girl}} t_j^2}
\end{cases}
在实际应用中，我们也可以将该问题写成矩阵形式，我们把模型写成
\underline{y} = \underline{X}\,\underline{\beta} + \epsilon
加下划线代表向量，这里 \underline{\beta} 表达式为
\underline{\beta}=\begin{pmatrix}
\beta_0\\[1mm]
\beta_{1b}\\[1mm]
\beta_{1g}
\end{pmatrix}
因此 \underline{X} 矩阵可得为
X = \begin{pmatrix}
1 & t_j\, \mathbb{I}_{\text{sex=boy}} & t_j\, \mathbb{I}_{\text{sex=girl}}
\end{pmatrix}
然后利用普通最小二乘公式
\underline{\beta} = (\underline{X}^\top \underline{X})^{-1} \underline{X}^\top \underline{y}
这样可以一次性求出包含
\beta_0,\ \beta_{1b},\ \beta_{1g}
的参数向量
\underline{\beta}
TQ3 模型 LM3：
LM3 模型更加复杂一些，设定为男孩和女孩完全不同，即假设斜率和截距均不同。
y_{ij}=\beta_{0b}\mathbb{I}_{\text{sex=boy}}+\beta_{0g}\mathbb{I}_{\text{sex=girl}}+\beta_{1b}t_j\mathbb{I}_{\text{sex=boy}}+\beta_{1g}t_j\mathbb{I}_{\text{sex=girl}}+\epsilon_{ij},\quad \epsilon_{ij}\sim\mathcal{N}(0,\sigma^2)
这里，指示函数
\mathbb{I}_{\text{sex=boy}}
和
\mathbb{I}_{\text{sex=girl}}
分别指示个体是否为男孩或女孩。这样做的目的是允许男孩和女孩在截距和斜率上均有不同表现，也就是说，它们的行为完全独立。
y_{ij} = \beta_{0b} + \beta_{1b}\,t_j + \epsilon_{ij}
y_{ij} = \beta_{0g} + \beta_{1g}\,t_j + \epsilon_{ij}
由于误差项服从正态分布，最大似然估计与最小二乘估计是一致的，因此，我们可以对男孩和女孩分别做简单线性回归来获得各自的参数估计。
对于男孩：
设男孩的观测数为
N_b
，并定义平均年龄 和 平均正畸测量值
\bar{t}_b = \frac{1}{N_b} \sum_{\text{boy}} t_j,\quad \bar{y}_b = \frac{1}{N_b} \sum_{\text{boy}} y_{ij}
根据简单线性回归的公式，男孩组的斜率和截距的估计为
\hat{\beta}_{0b}^{LM3} = \bar{y}_b - \hat{\beta}_{1b}^{LM3}\,\bar{t}_b
\hat{\beta}_{1b}^{LM3} = \frac{\sum_{\text{boy}} (t_j - \bar{t}_b)(y_{ij} - \bar{y}_b)}{\sum_{\text{boy}} (t_j - \bar{t}_b)^2}
对于女孩：
设女孩的观测数为
N_g
，并定义
\bar{t}_g = \frac{1}{N_g} \sum_{\text{girl}} t_j,\quad \bar{y}_g = \frac{1}{N_g} \sum_{\text{girl}} y_{ij}
同样，女孩组的斜率和截距的估计为
\hat{\beta}_{0g}^{LM3} = \bar{y}_g - \hat{\beta}_{1g}^{LM3}\,\bar{t}_g
\hat{\beta}_{1g}^{LM3} = \frac{\sum_{\text{girl}} (t_j - \bar{t}_g)(y_{ij} - \bar{y}_g)}{\sum_{\text{girl}} (t_j - \bar{t}_g)^2}
由于模型中男孩和女孩的参数完全独立（即截距与斜率均不同），在估计时两组数据的贡献不会相互干扰，所以两组数据之间相互独立，似然函数就可以分解为男孩和女孩两部分。这样，每一部分的最大似然估计可以单独用线性回归问题求解。
我们之前证明过，在正态分布误差的假设下，最大似然估计与最小二乘估计是一致的，因此，上述利用最小二乘法求得的估计值
\hat{\beta}_{0b}^{LM3}
、
\hat{\beta}_{1b}^{LM3}
、
\hat{\beta}_{0g}^{LM3}
和
\hat{\beta}_{1g}^{LM3}
同时也是最大似然估计。
TQ4 其他 LM 模型的变体
异方差模型：
男孩和女孩的残差变异性不同，可以允许误差项的方差依性别而变化
\epsilon_{ij} \sim
\begin{cases}
\mathcal{N}(0, \sigma_b^2) & \text{男孩} \\
\mathcal{N}(0, \sigma_g^2) & \text{女孩}
\end{cases}
非线性时间模型：
引入年龄的二次项，捕捉非线性增长趋势，不同性别的二次项系数不同。
y_{ij} = \beta_0 + (\beta_{1b} t_j + \beta_{2b} t_j^2) \mathbb{I}_{\text{sex=boy}} + (\beta_{1g} t_j + \beta_{2g} t_j^2) \mathbb{I}_{\text{sex=girl}} + \epsilon_{ij}
混合效应模型:
加入其他协变量，比如说牙齿发育情况等等，用 x_i 表示
y_{ij} = \beta_0 + \beta_{1b} t_j \mathbb{I}_{\text{sex=boy}} + \beta_{1g} t_j \mathbb{I}_{\text{sex=girl}} + \beta_2 x_i + \epsilon_{ij}
第 1 部分：处理局部变异性
接下来我们使用一项来自正畸研究的数据集，包含 16 名男孩和 11 名女孩，，年龄在 8 到 14 岁之间。响应变量是垂体与翼腭裂之间的距离（单位：毫米），在 8、10、12 和 14 岁时分别测量每位男孩和女孩的距离。
TQ5
传统的固定效应模型假设个体在所属性别组内的行为一致，包括 LM3 ，即使男孩女孩的截距和斜率不同，但是不同男孩、不同女孩之间的个体差异并没有考虑。这个情况下我们使用线性混合效应模型 LMM，来在整体趋势下引入个体行为变化，不仅能捕捉性别之间的差异，还可以捕捉同一性别内个体差异。换句话说，我们希望为每个个体的斜率和截距引入随机效应:
y_{ij} = (\beta_0 + \beta_{0r,i}) + (\beta_1 + \beta_{1r,i})t_j + \epsilon_{ij}
其中
- \beta_0 和 \beta_1 是总体的固定效应
- \beta_{0r,i} 表示第 i 个个体相对于总体截距 \beta_0 的随机偏差
- \beta_{1r,i} 表示第 i 个个体相对于总体斜率 \beta_1 的随机偏差
- \epsilon \sim \mathcal{N}(0, \sigma^2 I) 为全局残差
我们假设这些随机效应满足
\beta_{0r,i} \sim \mathcal{N}(0, \sigma_{0r}^2),\quad \beta_{1r,i} \sim \mathcal{N}(0, \sigma_{1r}^2)
每个个体的截距和斜率都是从一个分布中抽取的偏离值，其均值为零，方差反映了个体之间的差异大小。
注意这里是 固定效应 + 随机效应，这是因为个体的偏差在整体效应的基础上叠加不会改变变量之间的线性关系。换句话说，个体的特有偏差可以看作是总体效应的线性、独立的附加部分，这样既保留了总体趋势，又允许个体间存在随机差异。
TQ5.1
首先我们把上述 LMM 公式以向量形式重写成以下形式
\mathbf{y}_i =
\begin{bmatrix}
y_{i1} \\
\vdots \\
y_{iJ}
\end{bmatrix}
=
A
\begin{bmatrix}
\beta_0 \\
\beta_1
\end{bmatrix}
+
A
\begin{bmatrix}
\beta_{0r,i} \\
\beta_{1r,i}
\end{bmatrix}
+
\begin{bmatrix}
\epsilon_{i1} \\
\vdots \\
\epsilon_{iJ}
\end{bmatrix}
显而易见，A 矩阵为如下形式
A=\begin{bmatrix}
1 & t_1\\
1 & t_2\\
\vdots & \vdots\\
1 & t_J
\end{bmatrix}
原式变成:
\begin{bmatrix}
y_{i1} \\
y_{i2} \\
\vdots \\
y_{iJ}
\end{bmatrix}
=
\begin{bmatrix}
1 & t_1 \\
1 & t_2 \\
\vdots & \vdots \\
1 & t_J
\end{bmatrix}
\begin{bmatrix}
\beta_0 \\
\beta_1
\end{bmatrix}
+
\begin{bmatrix}
1 & t_1 \\
1 & t_2 \\
\vdots & \vdots \\
1 & t_J
\end{bmatrix}
\begin{bmatrix}
\beta_{0r,i} \\
\beta_{1r,i}
\end{bmatrix}
+
\begin{bmatrix}
\epsilon_{i1} \\
\epsilon_{i2} \\
\vdots \\
\epsilon_{iJ}
\end{bmatrix}
\underline{\beta}=\begin{bmatrix}
\beta_0\\[1mm]
\beta_1
\end{bmatrix}
\underline{b}_i=\begin{bmatrix}
\beta_{0r,i}\\[1mm]
\beta_{1r,i}
\end{bmatrix}
\underline{\epsilon}_i=\begin{bmatrix}
\epsilon_{i1}\\[1mm]
\epsilon_{i2}\\[1mm]
\vdots\\[1mm]
\epsilon_{iJ}
\end{bmatrix}
\underline{y_i} = A\,\underline{\beta} + A\,\underline{b_i} + \underline{\epsilon_i},\quad \epsilon_i\sim\mathcal{N}(0,\sigma^2I)
TQ5.2
描述: 为了进行最大似然估计，我们希望写出观测数据
\mathbf{y}_i
的似然函数。然而，由于
b_i
是未观察到的，我们需要通过对其积分（边缘化）来得到
\mathbf{y}_i
的边缘概率密度函数（将随机效应积分掉）：
p(\mathbf{y}_i \mid \theta) = \int p(\mathbf{y}_i \mid b_i, \theta)\, p(b_i \mid \theta)\, db_i
得到边缘概率密度后，我们可以构造整个样本的似然函数，从而对模型参数（包括固定效应和方差成分）进行估计。同时，求得
\mathbf{y}_i
的概率密度函数也是后续使用EM算法时构造完整数据对数似然、进行E步和M步的基础。
因此，对每个个体 i，将随机效应积分后可得 \mathbf{y}_i 的边缘分布为
\mathbf{y}_i\sim\mathcal{N}\Bigl(A\,\beta,\; A\,\Sigma_b\,A^T+\sigma^2I\Bigr)
于是，其概率密度函数写为
f\bigl(\mathbf{y}_i\mid\beta,\sigma_{0r}^2,\sigma_{1r}^2,\sigma^2\bigr)
=\frac{1}{(2\pi)^{J/2}\,\Bigl|A\,\Sigma_b\,A^T+\sigma^2I\Bigr|^{1/2}}
\exp\!\Biggl\{-\frac{1}{2}\Bigl(\mathbf{y}_i-A\,\beta\Bigr)^T\Bigl(A\,\Sigma_b\,A^T+\sigma^2I\Bigr)^{-1}\Bigl(\mathbf{y}_i-A\,\beta\Bigr)\Biggr\}
证明: 概率密度函数公式
首先我们给出条件分布:
\mathbf{y}_i \mid b_i \sim \mathcal{N}\Bigl(A\,\beta + A\,b_i,\;\sigma^2 I\Bigr)
随后我们利用正态分布的卷积性质求边缘分布
\underline{y}_i = A\,\underline{\beta} + \Bigl(A\,\underline{b}_i + \underline{\epsilon}_i\Bigr)
A\,\underline{b}_i 是随机效应部分，也服从正态分布
A\,\underline{b}_i \sim \mathcal{N}\Bigl(0,\; A\,\Sigma_b\,A^T\Bigr)
又因为
\underline{\epsilon}_i \sim \mathcal{N}\Bigl(0,\sigma^2 I\Bigr)
且 A\,\underline{b}_i 与 \underline{\epsilon}_i 独立，根据正态分布和的性质，有
A\,\underline{b}_i + \underline{\epsilon}_i \sim \mathcal{N}\Bigl(0,\; A\,\Sigma_b\,A^T + \sigma^2 I\Bigr)
因此，将 A\,\underline{\beta} 加上去，得到
\underline{y}_i \sim \mathcal{N}\Bigl(A\,\underline{\beta},\; A\,\Sigma_b\,A^T + \sigma^2 I\Bigr)
这就是将随机效应 \underline{b}_i 积分掉后得到的 \underline{y}_i 的边缘分布。证毕
然后我们写出 \underline{y}_i 的概率密度函数
多元正态分布的概率密度函数为：对于 J 维向量 \underline{y} \sim \mathcal{N}(\mu, \Sigma) :
f(\underline{y}) = \frac{1}{(2\pi)^{J/2} |\Sigma|^{1/2}} \exp\!\Bigl\{-\frac{1}{2} (\underline{y}-\mu)^T \Sigma^{-1} (\underline{y}-\mu)\Bigr\}
因此，对于个体 i，其概率密度函数为
f\bigl(\underline{y}_i\mid\underline{\beta},\sigma_{0r}^2,\sigma_{1r}^2,\sigma^2\bigr)
=\frac{1}{(2\pi)^{J/2}\,\Bigl|A\,\Sigma_b\,A^T+\sigma^2 I\Bigr|^{1/2}}
\exp\!\Biggl\{-\frac{1}{2}\Bigl(\underline{y}_i-A\,\underline{\beta}\Bigr)^T\Bigl(A\,\Sigma_b\,A^T+\sigma^2 I\Bigr)^{-1}\Bigl(\underline{y}_i-A\,\underline{\beta}\Bigr)\Biggr\}
证毕
由于这个公式是靠积分掉随机效应 \underline{b}_i 完成的，这就使得计算过程非常复杂，因此固定效应 \beta 的最大似然估计很难写出闭式解，所以我们下面要采用 EM 算法策略，我们将随机效应 \underline{b}_i 看作缺失数据，这样完整数据的对数似然函数具有较简单的形式，不需要再进行积分。
TQ5.3：基于 EM 算法求解参数
提出一种期望最大化（EM）算法，以接近 \hat{\beta}_0^{\text{LMM}},\ \hat{\beta}_1^{\text{LMM}},\ \hat{\sigma}^2,\ \hat{\sigma}_{0r}^2,\ \hat{\sigma}_{1r}^2,\ \hat{\rho}_{0r,1r},\ \hat{\sigma}^2 的最大似然估计值。
EM 算法的核心就是，我们没有完整数据，但是我们可以看作我们有完整数据，只不过缺了一部分，这样完整数据的对数似然函数就很简单。但是我们 还是缺一部分数据怎么办呢？我们就要在 E 步中计算缺失数据的条件期望，从而填补这些缺失数据，然后在 M 步中更新参数。
完全数据设定
在我们的线性混合效应模型中，将每个被试 i 的随机效应 \underline{b}_i 看作缺失数据，每个被试 i 的观测数据 \underline{y}_i 与随机效应 \underline{b}_i=\begin{bmatrix} \beta_{0r,i}\\ \beta_{1r,i} \end{bmatrix} 构成了完整数据，因此完整数据为
X_i=(\underline{y}_i,\, \underline{b}_i)
进而利用 EM 算法对参数 \theta=(\beta,\sigma_{0r}^2,\sigma_{1r}^2,\sigma^2) 进行最大似然估计。
完整数据对数似然函数
根据模型
\underline{y}_i=A\,\beta + A\,\underline{b}_i +\epsilon_i,\quad \epsilon_i\sim\mathcal{N}(0,\sigma^2I)
且随机效应
\underline{b}_i\sim\mathcal{N}\Bigl(0,\Sigma_b\Bigr),\quad \Sigma_b=\begin{bmatrix}\sigma_{0r}^2 & 0\\[1mm]0 & \sigma_{1r}^2\end{bmatrix}
每个 ID 的完整数据对数似然为
\begin{aligned}
\log p(\underline{y}_i, \underline{b}_i \mid \theta)
=&\; \log p(\underline{y}_i \mid \beta, \underline{b}_i, \sigma^2)
+\log p(\underline{b}_i\mid \sigma_{0r}^2,\sigma_{1r}^2)\\[1mm]
=&\; -\frac{J}{2}\log(2\pi\sigma^2)
-\frac{1}{2\sigma^2}\bigl\|\underline{y}_i-A\,\beta-A\,\underline{b}_i\bigr\|^2 -\frac{1}{2}\log\Bigl[(2\pi)^2|\Sigma_b|\Bigr]
-\frac{1}{2}\underline{b}_i^T\,\Sigma_b^{-1}\,\underline{b}_i.
\end{aligned}
证明: 每个个体的完整数据对数似然表达式
完整数据为 (\underline{y}_i,\underline{b}_i)，故其联合密度为‘
p(\underline{y}_i, \underline{b}_i \mid \theta)
=p(\underline{y}_i \mid \beta, \underline{b}_i, \sigma^2)\;p(\underline{b}_i \mid \Sigma_b)
（1）条件于 \underline{b}_i 的 \underline{y}_i 的概率密度函数为
p(\underline{y}_i \mid \beta, \underline{b}_i, \sigma^2)
=\frac{1}{(2\pi)^{J/2}(\sigma^2)^{J/2}}
\exp\!\Biggl\{-\frac{1}{2\sigma^2}\Bigl\|\underline{y}_i-A\,\beta-A\,\underline{b}_i\Bigr\|^2\Biggr\}
（2）随机效应 \underline{b}_i 的概率密度函数为
p(\underline{b}_i \mid \Sigma_b)
=\frac{1}{(2\pi)^{p/2}|\Sigma_b|^{1/2}}
\exp\!\Biggl\{-\frac{1}{2}\underline{b}_i^T\,\Sigma_b^{-1}\,\underline{b}_i\Biggr\}
因为 \underline{b}_i 是 2 维向量，其中 p=2，故有
p(\underline{b}_i \mid \Sigma_b)
=\frac{1}{(2\pi)|\Sigma_b|^{1/2}}
\exp\!\Biggl\{-\frac{1}{2}\underline{b}_i^T\,\Sigma_b^{-1}\,\underline{b}_i\Biggr\}
代入上述两个密度表达式到联合密度中
\begin{aligned}
p(\underline{y}_i, \underline{b}_i \mid \theta)
=&\; \frac{1}{(2\pi)^{J/2}(\sigma^2)^{J/2}}
\exp\!\Biggl\{-\frac{1}{2\sigma^2}\Bigl\|\underline{y}_i-A\,\beta-A\,\underline{b}_i\Bigr\|^2\Biggr\}
\times \frac{1}{(2\pi)|\Sigma_b|^{1/2}}
\exp\!\Biggl\{-\frac{1}{2}\underline{b}_i^T\,\Sigma_b^{-1}\,\underline{b}_i\Biggr\}
\end{aligned}
取自然对数得到
\begin{aligned}
\log p(\underline{y}_i, \underline{b}_i \mid \theta)
=&\; -\frac{J}{2}\log(2\pi)-\frac{J}{2}\log(\sigma^2)
-\frac{1}{2\sigma^2}\Bigl\|\underline{y}_i-A\,\beta-A\,\underline{b}_i\Bigr\|^2 -\log(2\pi)-\frac{1}{2}\log|\Sigma_b|
-\frac{1}{2}\underline{b}_i^T\,\Sigma_b^{-1}\,\underline{b}_i
\end{aligned}
合并后得到:
\begin{aligned}
\log p(\underline{y}_i, \underline{b}_i \mid \theta)
=&\; -\frac{J}{2}\log(2\pi\sigma^2)
-\frac{1}{2\sigma^2}\Bigl\|\underline{y}_i-A\,\beta-A\,\underline{b}_i\Bigr\|^2 -\frac{1}{2}\log\Bigl[(2\pi)^2|\Sigma_b|\Bigr]
-\frac{1}{2}\underline{b}_i^T\,\Sigma_b^{-1}\,\underline{b}_i
\end{aligned}
证明完毕
EM 算法的五步骤
(1) 初始化: 设初始参数
\theta^{(0)}=(\beta^{(0)},\sigma_{0r}^{2(0)},\sigma_{1r}^{2(0)},\sigma^{2(0)})
(2) 条件分布计算（E 步）: 在第 $m 次迭代时，利用当前参数 \theta^{(m)}，计算每个个体的随机效应 \underline{b}_i$ 的条件后验分布
p(\underline{b}_i\mid \underline{y}_i,\theta^{(m)})\sim \mathcal{N}\Bigl(\mu_i^{(m)},\,V_i^{(m)}\Bigr)
均值：
\mu_i^{(m)}=\Sigma_b^{(m)}\,A^T\Bigl(A\,\Sigma_b^{(m)}\,A^T+\sigma^{2(m)}I\Bigr)^{-1}\Bigl(\underline{y}_i-A\,\beta^{(m)}\Bigr)
协方差：
V_i^{(m)}=\Sigma_b^{(m)}- \Sigma_b^{(m)}\,A^T\Bigl(A\,\Sigma_b^{(m)}\,A^T+\sigma^{2(m)}I\Bigr)^{-1}A\,\Sigma_b^{(m)}
其中
\Sigma_b^{(m)}=\begin{bmatrix}\sigma_{0r}^{2(m)} & 0\\[1mm]0 & \sigma_{1r}^{2(m)}\end{bmatrix}
证明: 条件后验分布的均值和协方差公式
首先引入多元正态分布的条件分布公式
有一个联合正态分布的随机向量
\begin{pmatrix}
X \\
Y
\end{pmatrix}
\sim \mathcal{N}\Biggl(
\begin{pmatrix}
\mu_X \\
\mu_Y
\end{pmatrix},\,
\begin{pmatrix}
\Sigma_{XX} & \Sigma_{XY} \\
\Sigma_{YX} & \Sigma_{YY}
\end{pmatrix}
\Biggr)
那么根据多元正态分布的性质，条件分布
X\mid Y=y\sim\mathcal{N}\Bigl(\mu_X+\Sigma_{XY}\Sigma_{YY}^{-1}(y-\mu_Y),\; \Sigma_{XX}-\Sigma_{XY}\Sigma_{YY}^{-1}\Sigma_{YX}\Bigr)
在我们的混合效应模型中，我们把随机效应 \underline{b}_i 看作 X，而观测数据 \underline{y}_i 看作 Y。
\underline{b}_i\sim\mathcal{N}\Bigl(0,\Sigma_b\Bigr)
\underline{y}_i \sim \mathcal{N}\Bigl(A\,\underline{\beta},\; A\,\Sigma_b\,A^T + \sigma^2 I\Bigr)
可以推导出联合分布为
\begin{pmatrix}
\underline{b}_i \\
\underline{y}_i
\end{pmatrix}
\sim \mathcal{N}\Biggl(
\begin{pmatrix}
0 \\
A\,\beta
\end{pmatrix},\,
\begin{pmatrix}
\Sigma_b & \Sigma_b\,A^T \\
A\,\Sigma_b & A\,\Sigma_b\,A^T+\sigma^2 I
\end{pmatrix}
\Biggr)
对比一下可以得到:
条件均值
\mu_i = E\bigl(\underline{b}_i \mid \underline{y}_i\bigr)
=\mu_X+\Sigma_{XY}\Sigma_{YY}^{-1}(\underline{y}_i-\mu_Y)
=\Sigma_b\,A^T\Bigl(A\,\Sigma_b\,A^T+\sigma^2 I\Bigr)^{-1}\Bigl(\underline{y}_i-A\,\beta\Bigr)
条件协方差
V_i = \operatorname{Var}\bigl(\underline{b}_i \mid \underline{y}_i\bigr)
=\Sigma_{XX}-\Sigma_{XY}\Sigma_{YY}^{-1}\Sigma_{YX}
=\Sigma_b- \Sigma_b\,A^T\Bigl(A\,\Sigma_b\,A^T+\sigma^2 I\Bigr)^{-1}A\,\Sigma_b
证明完毕。
(3) 构建 Q 函数（E 步）: 利用完整数据的对数似然，对缺失的随机效应 \underline{b}_i 求条件期望，得到
Q\bigl(\theta\mid\theta^{(m)}\bigr)
=\sum_{i=1}^I E_{\underline{b}_i\mid \underline{y}_i,\theta^{(m)}}\Bigl[\log p(\underline{y}_i,\underline{b}_i\mid\theta)\Bigr]
对完整数据对数似然的分解，利用观测数据与随机效应独立（个体之间独立）的性质，Q 函数可以写成各被试贡献的和：
Q\bigl(\theta\mid\theta^{(m)}\bigr)=\sum_{i=1}^I Q_i\bigl(\theta\mid\theta^{(m)}\bigr)
(4) 最大化 Q 函数（M 步）: 在固定 Q\bigl(\theta\mid\theta^{(m)}\bigr) 的情况下，对 \theta 进行优化：
\theta^{(m+1)}=\arg\max_{\theta}\Bigl\{Q\bigl(\theta\mid\theta^{(m)}\bigr)+\log p(\theta)\Bigr\}
这里由于没有引入先验信息，所以不需要 \log p(\theta)。经过对 Q 函数求导，可以得到更新公式：
-
固定效应参数更新：
\beta^{(m+1)}=\Bigl(\sum_{i=1}^I A^T A\Bigr)^{-1}\sum_{i=1}^I A^T\Bigl(\underline{y}_i-A\,\mu_i^{(m)}\Bigr)
-
随机效应方差更新（分别对截距和斜率）：
\sigma_{0r}^{2(m+1)}=\frac{1}{I}\sum_{i=1}^I\Bigl[(\mu_i^{(m)})_1^2+\bigl(V_i^{(m)}\bigr)_{11}\Bigr]
\sigma_{1r}^{2(m+1)}=\frac{1}{I}\sum_{i=1}^I\Bigl[(\mu_i^{(m)})_2^2+\bigl(V_i^{(m)}\bigr)_{22}\Bigr]
-
残差方差更新：
\sigma^{2(m+1)}=\frac{1}{IJ}\sum_{i=1}^I\Biggl\{
\bigl\|\underline{y}_i-A\,\beta^{(m+1)}-A\,\mu_i^{(m)}\bigr\|^2
+\operatorname{tr}\Bigl(A\,V_i^{(m)}\,A^T\Bigr)
\Biggr\}
-
其中
- \mu_i^{(m)}=E\bigl[\underline{b}_i\mid \underline{y}_i,\theta^{(m)}\bigr]
- V_i^{(m)}=\operatorname{Var}\bigl[\underline{b}_i\mid \underline{y}_i,\theta^{(m)}\bigr]
证明: 固定效应参数 和 随机效应方差 和 残差 公式
(1) 固定效应参数 \beta 的更新
完整数据对数似然函数如下:
\begin{aligned}
\log p(\underline{y}_i, \underline{b}_i \mid \theta)
=&\; -\frac{J}{2}\log(2\pi\sigma^2)
-\frac{1}{2\sigma^2}\Bigl\|\underline{y}_i-A\,\beta-A\,\underline{b}_i\Bigr\|^2 -\frac{1}{2}\log\Bigl[(2\pi)^2|\Sigma_b|\Bigr]
-\frac{1}{2}\underline{b}_i^T\,\Sigma_b^{-1}\,\underline{b}_i
\end{aligned}
与 \beta 相关的部分出现在
-\frac{1}{2\sigma^2}\sum_{i=1}^I \Bigl\|\underline{y}_i-A\,\beta-A\,\underline{b}_i\Bigr\|^2
对 \underline{b}_i 求条件期望后
E\Bigl[\|\underline{y}_i-A\,\beta-A\,\underline{b}_i\|^2\Bigr]
=\Bigl\|\underline{y}_i-A\,\beta-A\,\mu_i^{(m)}\Bigr\|^2
+\operatorname{tr}\Bigl(A\,V_i^{(m)}\,A^T\Bigr)
其中
- \mu_i^{(m)}=E\bigl[\underline{b}_i\mid \underline{y}_i,\theta^{(m)}\bigr]
- V_i^{(m)}=\operatorname{Var}\bigl[\underline{b}_i\mid \underline{y}_i,\theta^{(m)}\bigr]
注意到 \operatorname{tr}\Bigl(A\,V_i^{(m)}\,A^T\Bigr) 与 \beta 无关，因此优化 Q 关于 \beta 就等价于最小化
\sum_{i=1}^I \Bigl\|\underline{y}_i-A\,\beta-A\,\mu_i^{(m)}\Bigr\|^2
这是一个标准的最小二乘问题，其最优解满足
\sum_{i=1}^I A^T\Bigl(\underline{y}_i-A\,\beta-A\,\mu_i^{(m)}\Bigr)=0
解得
\beta^{(m+1)}=\Bigl(\sum_{i=1}^I A^T A\Bigr)^{-1}\sum_{i=1}^I A^T\Bigl(\underline{y}_i-A\,\mu_i^{(m)}\Bigr)
(2) 随机效应方差 \sigma_{0r}^2 和 \sigma_{1r}^2 的更新
完整数据对数似然中关于 \underline{b}_i 的部分为
-\frac{1}{2}\sum_{i=1}^I\Bigl[\log|\Sigma_b|+\underline{b}_i^T\,\Sigma_b^{-1}\,\underline{b}_i\Bigr]
对每个个体取条件期望后，该部分变为
-\frac{1}{2}\sum_{i=1}^I\left[\log|\Sigma_b|+\operatorname{tr}\Bigl(\Sigma_b^{-1}\bigl(\mu_i^{(m)}\mu_i^{(m)T}+V_i^{(m)}\bigr)\Bigr)\right]
由于 \Sigma_b 是对角矩阵，即
\Sigma_b=\begin{pmatrix}\sigma_{0r}^2 & 0\\[1mm]0 & \sigma_{1r}^2\end{pmatrix}
所以
\log|\Sigma_b|=\log(\sigma_{0r}^2)+\log(\sigma_{1r}^2)
且
\operatorname{tr}\Bigl(\Sigma_b^{-1}\bigl(\mu_i^{(m)}\mu_i^{(m)T}+V_i^{(m)}\bigr)\Bigr)
=\frac{(\mu_i^{(m)})_1^2+(V_i^{(m)})_{11}}{\sigma_{0r}^2}
+\frac{(\mu_i^{(m)})_2^2+(V_i^{(m)})_{22}}{\sigma_{1r}^2}
将对所有个体求和后，针对 \sigma_{0r}^2 的部分
Q(\sigma_{0r}^2)=-\frac{I}{2}\log(\sigma_{0r}^2)
-\frac{1}{2\sigma_{0r}^2}\sum_{i=1}^I\Bigl[(\mu_i^{(m)})_1^2+(V_i^{(m)})_{11}\Bigr]
令对 \sigma_{0r}^2 求导并置零：
-\frac{I}{2\sigma_{0r}^2}+\frac{1}{2(\sigma_{0r}^2)^2}\sum_{i=1}^I\Bigl[(\mu_i^{(m)})_1^2+(V_i^{(m)})_{11}\Bigr]=0
乘以 2(\sigma_{0r}^2)^2 得
-I\,\sigma_{0r}^2+\sum_{i=1}^I\Bigl[(\mu_i^{(m)})_1^2+(V_i^{(m)})_{11}\Bigr]=0
因此，
\sigma_{0r}^{2(m+1)}=\frac{1}{I}\sum_{i=1}^I\Bigl[(\mu_i^{(m)})_1^2+(V_i^{(m)})_{11}\Bigr]
类似地，对 \sigma_{1r}^2 有
\sigma_{1r}^{2(m+1)}=\frac{1}{I}\sum_{i=1}^I\Bigl[(\mu_i^{(m)})_2^2+(V_i^{(m)})_{22}\Bigr]
(3) 残差方差 \sigma^2 的更新
完整数据对数似然中关于 \sigma^2 的部分主要为
-\frac{J}{2}\sum_{i=1}^I\log(2\pi\sigma^2)
-\frac{1}{2\sigma^2}\sum_{i=1}^I E\Bigl[\|\underline{y}_i-A\,\beta-A\,\underline{b}_i\|^2\Bigr]
在对 \underline{b}_i 求条件期望后，我们有
E\Bigl[\|\underline{y}_i-A\,\beta-A\,\underline{b}_i\|^2\Bigr]
=\|\underline{y}_i-A\,\beta-A\,\mu_i^{(m)}\|^2+\operatorname{tr}\Bigl(A\,V_i^{(m)}\,A^T\Bigr)
将对所有个体求和，并记总观测数为 IJ，则优化关于 \sigma^2 的部分为
-\frac{IJ}{2}\log(\sigma^2)
-\frac{1}{2\sigma^2}\sum_{i=1}^I\left\{\|\underline{y}_i-A\,\beta-A\,\mu_i^{(m)}\|^2+\operatorname{tr}\Bigl(A\,V_i^{(m)}\,A^T\Bigr)\right\}
对 \sigma^2 求导并令导数为零，得到
\sigma^{2(m+1)}=\frac{1}{IJ}\sum_{i=1}^I\left\{\|\underline{y}_i-A\,\beta^{(m+1)}-A\,\mu_i^{(m)}\|^2+\operatorname{tr}\Bigl(A\,V_i^{(m)}\,A^T\Bigr)\right\}
(5) 检查收敛：判断参数 \theta^{(m+1)} 与 \theta^{(m)} 是否足够接近（对数似然变化小于某一阈值
\epsilon
），若满足则停止迭代，否则返回 Step 2。
综上所述，我们在 E 步中先猜测隐藏数据（求出完整数据 \underline{b}_i 的条件分布或其期望），在 M 步中对参数 \theta 进行极大化，交替进行，尽管 EM 算法可能收敛慢、局部最优，但是它可以保证每一步的似然 p(\underline{y}\mid \theta) 不会减小.
TQ5.4：基于 EM 算法估计值对 个体化拟合的说明
经过 EM 算法迭代，我们得到固定效应参数和方差成分的估计值
\hat{\beta}=\begin{bmatrix}\hat{\beta}_0^{\text{LMM}}\\ \hat{\beta}_1^{\text{LMM}}\end{bmatrix},\quad \hat{\sigma}_{0r}^2,\quad \hat{\sigma}_{1r}^2,\quad \hat{\sigma}^2
第 i 个个体的线性拟合为固定效应与随机效应之和：
\hat{y}_{ij} = \Bigl(\hat{\beta}_0^{\text{LMM}} + \hat{\beta}_{0r,i}\Bigr) + \Bigl(\hat{\beta}_1^{\text{LMM}} + \hat{\beta}_{1r,i}\Bigr)\,t_j + \epsilon_{ij}
其中：
- \hat{\beta}_0^{\text{LMM}}, \hat{\beta}_1^{\text{LMM}} 为总体固定效应，反映全体个体的平均趋势
- \hat{\beta}_{0r,i}, \hat{\beta}_{1r,i} 为个体随机效应，表示该个体相对于总体趋势的偏差
个体随机效应向量服从分布
\begin{pmatrix}\hat{\beta}_{0r,i}\\ \hat{\beta}_{1r,i}\end{pmatrix} \sim \mathcal{N}\Bigl(0,\hat{\Sigma}_b\Bigr)
其协方差矩阵体现了截距与斜率之间的相关性
\hat{\Sigma}_b = \begin{pmatrix}\hat{\sigma}_{0r}^2 & \hat{\rho}_{0r,1r}\sqrt{\hat{\sigma}_{0r}^2\,\hat{\sigma}_{1r}^2}\\[1mm]\hat{\rho}_{0r,1r}\sqrt{\hat{\sigma}_{0r}^2\,\hat{\sigma}_{1r}^2} & \hat{\sigma}_{1r}^2\end{pmatrix}
这种模型在全局趋势的基础上，通过引入个体随机效应捕捉个体间的差异，允许不同个体的截距和斜率围绕总体均值波动，从而更灵活地描述非规则数据。
第 2 部分：处理鲁棒性
在第 1 部分中，我们将噪声 \epsilon 建模为高斯分布。为了应对可能出现的异常值，使模型更加鲁棒，可以将噪声和/或随机效应建模为 Student's t 分布。
Student's t 分布是一种具有较重尾部（heavy-tailed）的连续概率分布，其形式与正态分布类似，但在尾部衰减得较慢。它由自由度参数 \eta 控制，当 \eta 越小时，尾部越厚，而当 \eta \to \infty 时，Student's t 分布收敛到正态分布。
- 正态分布的尾部衰减很快，对于极端异常值的容忍度低，Student's t 分布尾部衰减较慢，如果存在异常值，其概率密度不会突然骤降，因此模型对异常数据的惩罚较小，使得异常值对整体模型的影响减弱，从而提升模型鲁棒性。
TQ6：观测值和随机效应的联合分布
在原模型中，我们有
\mathbf{y}_i = A \beta + A b_i + \epsilon_i
其中噪声 \epsilon_i 和随机效应 b_i 均假设服从自由度 \eta_i > 2 的 Student's t 分布
\epsilon_i \sim t_{\eta_i}(0,\sigma^2 I), \quad b_i \sim t_{\eta_i}(0,\Sigma_b)
假设一个 d 维随机向量 x 服从位置参数为 \mu 、尺度矩阵为 \Sigma 、自由度为 \eta_i 的 Student's t 分布，其密度函数可以写作
p(x) = \frac{\Gamma\Bigl(\frac{\eta_i+d}{2}\Bigr)}{\Gamma\Bigl(\frac{\eta_i}{2}\Bigr)(\eta_i\pi)^{d/2}|\Sigma|^{1/2}} \Biggl[1+\frac{1}{\eta_i}(x-\mu)^T\Sigma^{-1}(x-\mu)\Biggr]^{-\frac{\eta_i+d}{2}}
(1) 对 b_i 的密度推导
- 位置参数 \mu=0
- 尺度矩阵 \Sigma=\Sigma_b
- 向量维数 d=2
将这些量代入上式，得到随机效应的边缘分布
p(b_i) = \frac{\Gamma\Bigl(\frac{\eta_i+2}{2}\Bigr)}{\Gamma\Bigl(\frac{\eta_i}{2}\Bigr)(\eta_i\pi)^{1}|\Sigma_b|^{1/2}} \Biggl[1+\frac{1}{\eta_i}b_i^T\Sigma_b^{-1}b_i\Biggr]^{-\frac{\eta_i+2}{2}}
(2) 对条件下的 \mathbf{y}_i 密度推导
在原模型中，我们有
\mathbf{y}_i = A \beta + A b_i + \epsilon_i
有
- 位置参数 \mu = \mu_y = A \beta + A b_i
- 尺度矩阵 \Sigma = \sigma^2 I
- 向量维数 d = J
将这些量代入上式，得到给定随机效应 b_i 后观测值的条件分布
p(\mathbf{y}_i \mid b_i) = \frac{\Gamma\Bigl(\frac{\eta_i+J}{2}\Bigr)}{\Gamma\Bigl(\frac{\eta_i}{2}\Bigr)(\eta_i\pi)^{J/2}|\sigma^2I|^{1/2}} \Biggl[1+\frac{1}{\eta_i}(\mathbf{y}_i-\mu_y)^T(\sigma^2I)^{-1}(\mathbf{y}_i-\mu_y)\Biggr]^{-\frac{\eta_i+J}{2}}
化简一下，得到
p(\mathbf{y}_i \mid b_i) = \frac{\Gamma\Bigl(\frac{\eta_i+J}{2}\Bigr)}{\Gamma\Bigl(\frac{\eta_i}{2}\Bigr)(\eta_i\pi)^{J/2}(\sigma^2)^{J/2}} \Biggl[1+\frac{1}{\eta_i}\frac{\|\mathbf{y}_i-A\beta-A b_i\|^2}{\sigma^2}\Biggr]^{-\frac{\eta_i+J}{2}}
因此，观测值 \mathbf{y}_i 与随机效应 b_i 的联合分布为（链式法则）
p(\mathbf{y}_i, b_i) = p(b_i) \, p(\mathbf{y}_i \mid b_i)
TQ7：\mathbf{y}_i 的独立性及其边际分布
由于不同个体 i 之间的随机效应和误差均相互独立，则每个个体的完整数据 (\mathbf{y}_i, b_i) 也是独立的，从而联合分布之间也相互独立。
要得到观测值 \mathbf{y}_i 的边际分布，我们只需要对随机效应积分，即
p(\mathbf{y}_i) = \int p(\mathbf{y}_i, b_i) \, db_i = \int p(b_i) \, p(\mathbf{y}_i \mid b_i) \, db_i
TQ8：尺度混合表示及条件分布
我们知道，对于自由度 \eta_i 的 Student's t 分布，可以表示为尺度混合模型。引入隐藏变量（纹理参数） \tau_i ，满足
\tau_i \sim \operatorname{Gamma}\Bigl(\frac{\eta_i}{2}, \frac{\eta_i}{2}\Bigr)
则有以下条件分布说明：
(1) 观测值与随机效应联合条件于 \tau_i
条件下的随机效应
b_i \mid \tau_i \sim \mathcal{N}\Bigl(0, \frac{\Sigma_b}{\tau_i}\Bigr)
利用线性变换的性质
A b_i \mid \tau_i \sim \mathcal{N}\Bigl(0, \frac{A\Sigma_bA^T}{\tau_i}\Bigr)
条件下的随机噪音
\epsilon_i \mid \tau_i \sim \mathcal{N}\Bigl(0, \frac{\sigma^2 I}{\tau_i}\Bigr)
因此，由于
\mathbf{y}_i = A \beta + A b_i + \epsilon_i
观测值的条件分布为
\mathbf{y}_i \mid \tau_i \sim \mathcal{N}\Bigl(A \beta, \; \frac{1}{\tau_i}\Bigl(A\Sigma_bA^T+\sigma^2 I\Bigr)\Bigr)
因此，条件于 \tau_i 下的联合分布为
\begin{pmatrix}
b_i \\[1mm]
\mathbf{y}_i
\end{pmatrix} \Bigg| \tau_i \sim \mathcal{N}\Biggl(\begin{pmatrix}
0 \\[1mm]
A \beta
\end{pmatrix},\; \frac{1}{\tau_i}\begin{pmatrix}
\Sigma_b & \Sigma_b A^T \\[1mm]
A\Sigma_b & A\Sigma_bA^T+\sigma^2 I
\end{pmatrix}\Biggr)
(2) 观测值条件于 \tau_i
前面已得:
\mathbf{y}_i \mid \tau_i \sim \mathcal{N}\Bigl(A \beta, \; \frac{1}{\tau_i}\Bigl(A\Sigma_bA^T+\sigma^2 I\Bigr)\Bigr)
(3) 观测值条件于随机效应与 \tau_i
条件于 b_i 和 \tau_i ，有
\mathbf{y}_i \mid b_i, \tau_i \sim \mathcal{N}\Bigl(A \beta + A b_i, \; \frac{\sigma^2 I}{\tau_i}\Bigr)
TQ9：基于 EM 算法的最大似然估计
在该模型中，每个个体 i 的观测数据和隐藏数据构成完整数据
X_i = \{\mathbf{y}_i,\, b_i,\, \tau_i\}
其中
- \mathbf{y}_i 为第 i 个个体在所有时间点上的观测值
- b_i 为个体随机效应
- \tau_i 为尺度混合中引入的纹理参数，其先验分布为
\tau_i \sim \operatorname{Gamma}\Bigl(\frac{\eta_i}{2}, \frac{\eta_i}{2}\Bigr)
也就是说，我们不仅要对随机效应 b_i 进行处理，还引入了尺度混合表示中的隐藏变量（纹理参数） \tau_i，这体现在 E 步中，我们不仅需要计算
\mu_i = E\bigl(b_i \mid \mathbf{y}_i, \theta^{(m)}\bigr), \quad
V_i = \operatorname{Var}\bigl(b_i \mid \mathbf{y}_i, \theta^{(m)}\bigr)
还要计算纹理参数的条件期望
w_i = E\bigl(\tau_i \mid \mathbf{y}_i, \theta^{(m)}\bigr)
我们将模型假设为
\mathbf{y}_i = A \beta + A b_i + \epsilon_i
其中
\epsilon_i \mid \tau_i \sim \mathcal{N}\Bigl(0, \frac{\sigma^2 I}{\tau_i}\Bigr)
b_i \mid \tau_i \sim \mathcal{N}\Bigl(0, \frac{\Sigma_b}{\tau_i}\Bigr)
因此，完整数据的对数似然可写为
\log p(\mathbf{y}_i, b_i, \tau_i \mid \theta)
=\log p(\mathbf{y}_i \mid b_i, \tau_i, \beta, \sigma^2)
+\log p(b_i \mid \tau_i, \Sigma_b)
+\log p(\tau_i)
其中 \theta = (\beta, \Sigma_b, \sigma^2) 代表未知参数，具体来说:
- \beta 表示固定效应参数
- \Sigma_b 表示随机效应的协方差矩阵，描述 b_i 的方差结构
- \sigma^2 表示 \epsilon_i 噪声（残差）的方差
- 隐藏变量 \tau_i 的先验分布是固定形式 \tau_i \sim \operatorname{Gamma}\Bigl(\frac{\eta_i}{2}, \frac{\eta_i}{2}\Bigr) 与 \theta 无关
根据 EM 算法原理，我们将目标写成最大化完整数据对数似然的条件期望，即构造 Q 函数
Q(\theta \mid \theta^{(m)}) = \sum_{i=1}^I E_{(b_i, \tau_i) \mid \mathbf{y}_i, \theta^{(m)}}\Bigl[\log p(\mathbf{y}_i, b_i, \tau_i \mid \theta)\Bigr]
利用独立同分布的假设，该 Q 函数可以分解为各个个体的贡献之和
Q(\theta \mid \theta^{(m)}) = \sum_{i=1}^I Q_i(\theta \mid \theta^{(m)})
EM 算法步骤如下：
(1) 初始化 : 设初始参数 \theta^{(0)} = (\beta^{(0)}, \Sigma_b^{(0)}, \sigma^{2(0)})
(2) E 步: 在第 m 次迭代中，基于当前参数 \theta^{(m)}，需要计算每个个体 i 下隐藏数据 (b_i, \tau_i) 的条件后验分布，即
p(b_i, \tau_i \mid \mathbf{y}_i, \theta^{(m)})
更具体来说，计算以下分量
\mu_i^{(m)} = E\bigl(b_i \mid \mathbf{y}_i, \theta^{(m)}\bigr), \quad
V_i^{(m)} = \operatorname{Var}\bigl(b_i \mid \mathbf{y}_i, \theta^{(m)}\bigr)
w_i^{(m)} = E\bigl(\tau_i \mid \mathbf{y}_i, \theta^{(m)}\bigr)
(3) 构造 Q 函数:
将 E 步中求得的条件分布代入完整数据对数似然的期望，得到
Q(\theta \mid \theta^{(m)}) = \sum_{i=1}^I E_{(b_i, \tau_i) \mid \mathbf{y}_i, \theta^{(m)}}\Bigl[\log p(\mathbf{y}_i, b_i, \tau_i \mid \theta)\Bigr]
(4) M 步
固定 Q 函数后，对参数 \theta = (\beta, \Sigma_b, \sigma^2) 进行极大化，即求
\theta^{(m+1)} = \arg\max_{\theta} \Bigl\{ Q(\theta \mid \theta^{(m)}) \Bigr\}
经过对 Q 函数求导并令导数为零，可以得到更新公式，例如：
\beta^{(m+1)} = \Bigl(\sum_{i=1}^I A^T \, w_i^{(m)} A \Bigr)^{-1} \Bigl(\sum_{i=1}^I A^T \, w_i^{(m)} \Bigl(\mathbf{y}_i - A \mu_i^{(m)}\Bigr) \Bigr)
\Sigma_b^{(m+1)} = \frac{1}{I} \sum_{i=1}^I w_i^{(m)} \Bigl[\mu_i^{(m)} \bigl(\mu_i^{(m)}\bigr)^T + V_i^{(m)}\Bigr]
\sigma^{2(m+1)} = \frac{1}{IJ} \sum_{i=1}^I w_i^{(m)} \Biggl\{ \Bigl\|\mathbf{y}_i - A \beta^{(m+1)} - A \mu_i^{(m)}\Bigr\|^2 + \operatorname{tr}\Bigl(A \, V_i^{(m)} A^T\Bigr) \Biggr\}
(5) 检查收敛
基于完整数据 \{\mathbf{y}_i,\, b_i,\, \tau_i\} , 我们利用 EM 算法将原本难以直接求解的最大似然估计问题转化为交替求解“隐藏数据”（即随机效应和纹理参数）的后验期望和参数更新的问题。利用完整数据对数似然的分解及其在样本间的独立性，我们可以将 Q 函数写成各个个体贡献的和，从而大大简化计算。经过多次迭代后，算法保证每一步使观测数据的似然值不减，最终收敛到一个局部最优解，即为最大似然估计值.
第 3 部分：处理缺失值
参加实验的测试者不再来了，那么数据就缺失了，我们假设数据完全是 随机缺失 (MCAR) 的。
TQ10: 使用条件期望塔性质设计缺失数据的EM算法——理论推导
考虑每个个体 i 的观测向量分为观察部分与缺失部分：
\mathbf{y}_i = \begin{pmatrix}\mathbf{y}_{i,obs} \\ \mathbf{y}_{i,mis}\end{pmatrix}
其中模型为
\mathbf{y}_i = A\,\beta + A\,b_i + \epsilon_i,\quad \epsilon_i\sim\mathcal{N}(0,\sigma^2I)
将缺失数据看作隐变量，则完整数据为 \{\mathbf{y}_{i,obs},\mathbf{y}_{i,mis}, b_i\} ，则完整数据的对数似然为
\log p(\mathbf{y}_i, b_i\mid \theta) = -\frac{J}{2}\log(2\pi\sigma^2) -\frac{1}{2\sigma^2}\| \mathbf{y}_i - A\beta - A\,b_i \|^2 -\frac{1}{2}\log\Bigl[(2\pi)^2|\Sigma_b|\Bigr] -\frac{1}{2}b_i^T\Sigma_b^{-1}b_i
由于数据缺失，我们利用条件期望的塔性质（定理2.5.1, Anderson,1975），构造Q函数：
Q(\theta\mid\theta^{(m)})=\sum_{i=1}^I E_{\mathbf{y}_{i,mis},\,b_i\mid \mathbf{y}_{i,obs},\theta^{(m)}}\Bigl[\log p(\mathbf{y}_i, b_i\mid \theta)\Bigr]
在E步中，我们需要计算两部分隐变量的条件后验分布：
- 对于随机效应 b_i ，其条件分布已在TQ5中给出；
- 对于缺失观测值 \mathbf{y}_{i,mis} ，利用多元正态分布的性质，将其条件于观察部分和当前参数估计的条件期望与协方差计算出来：
E(\mathbf{y}_{i,mis}\mid \mathbf{y}_{i,obs},\theta^{(m)}) = \mu_{mis}^{(m)} + \Sigma_{mis,obs}^{(m)}\Bigl(\Sigma_{obs,obs}^{(m)}\Bigr)^{-1}(\mathbf{y}_{i,obs}-\mu_{obs}^{(m)})
接着，在M步中，将E步中计算得到的条件期望替代缺失值，构造 完整数据的对数似然，然后最大化以更新参数\theta=(\beta,\Sigma_b,\sigma^2)
该算法在E步中针对每个个体：
- 若存在缺失值，则利用当前参数计算缺失部分的条件期望进行 填补
- 同时对随机效应 b_i 计算条件均值与协方差
- M步中利用 完整 数据更新参数
EM 算法的完整步骤
(1) 初始化
\theta^{(0)} = (\beta^{(0)}, \Sigma_b^{(0)}, \sigma^{2(0)})
如果存在缺失数据，则利用插值方法（或其他）来进行初始填补，得到完整数据 \{\mathbf{y}_{i,obs},\mathbf{y}_{i,mis}^{(0)}\}
(2) E 步
在第 m 次迭代中，对于每个个体 i，将缺失数据 \mathbf{y}_{i,mis} 视为隐变量，其完整数据为
X_i = \{\mathbf{y}_{i,obs}, \mathbf{y}_{i,mis}, b_i\}
我们需要计算条件期望
E_{\mathbf{y}_{i,mis},\,b_i \mid \mathbf{y}_{i,obs}, \theta^{(m)}}\Bigl[\log p(\mathbf{y}_i, b_i \mid \theta)\Bigr]
这要求计算下列两个部分的后验条件期望：
(a) 随机效应 b_i 的后验均值和协方差
已知在完整数据模型中 b_i 的条件分布可得，记
\mu_i^{(m)} = E\bigl(b_i \mid \mathbf{y}_{i,obs}, \theta^{(m)}\bigr), \quad V_i^{(m)} = \operatorname{Var}\bigl(b_i \mid \mathbf{y}_{i,obs}, \theta^{(m)}\bigr)
(b) 缺失观测值 \mathbf{y}_{i,mis} 的条件期望与协方差
利用多元正态分布的性质，将 \mathbf{y}_i 分为观察部分与缺失部分
\mathbf{y}_i = \begin{pmatrix}\mathbf{y}_{i,obs} \\[1mm] \mathbf{y}_{i,mis}\end{pmatrix}
并记该正态分布的均值向量分为 \mu_{obs} 与 \mu_{mis}，协方差矩阵分块为
\Sigma = \begin{pmatrix} \Sigma_{obs,obs} & \Sigma_{obs,mis} \\[1mm] \Sigma_{mis,obs} & \Sigma_{mis,mis} \end{pmatrix}
则利用正态分布的条件分布公式，缺失部分的条件期望为
E\bigl(\mathbf{y}_{i,mis} \mid \mathbf{y}_{i,obs}, \theta^{(m)}\bigr) = \mu_{mis}^{(m)} + \Sigma_{mis,obs}^{(m)}\Bigl(\Sigma_{obs,obs}^{(m)}\Bigr)^{-1}\Bigl(\mathbf{y}_{i,obs}-\mu_{obs}^{(m)}\Bigr)
(3) 构造 Q 函数
将目标写成最大化完整数据对数似然的条件期望，即构造
Q(\theta \mid \theta^{(m)}) = \sum_{i=1}^I E_{\mathbf{y}_{i,mis},\,b_i \mid \mathbf{y}_{i,obs}, \theta^{(m)}}\Bigl[\log p(\mathbf{y}_i, b_i \mid \theta)\Bigr]
其中完整数据的对数似然为
\log p(\mathbf{y}_i, b_i\mid \theta) = -\frac{J}{2}\log(2\pi\sigma^2) -\frac{1}{2\sigma^2}\| \mathbf{y}_i - A\beta - A\,b_i \|^2 -\frac{1}{2}\log\Bigl[(2\pi)^2|\Sigma_b|\Bigr] -\frac{1}{2}b_i^T\Sigma_b^{-1}b_i
(4) M 步
固定 Q 函数后，对参数 \theta = (\beta, \Sigma_b, \sigma^2) 进行极大化，即求
\theta^{(m+1)} = \arg\max_\theta \Bigl\{ Q(\theta \mid \theta^{(m)}) \Bigr\}
其三大参数公式与 TQ9 结果完全相同:
-
固定效应更新：
\beta^{(m+1)} = \Bigl(\sum_{i=1}^I A^T\, w_i^{(m)} A\Bigr)^{-1}\Bigl(\sum_{i=1}^I A^T\, w_i^{(m)}\Bigl(\mathbf{y}_i-A\,\mu_i^{(m)}\Bigr)\Bigr),
其中 w_i^{(m)} 是缺失数据部分中（以及 b_i 部分）的权重因子，通常来自于 \tau_i 的后验条件期望。
-
随机效应协方差更新：
\Sigma_b^{(m+1)} = \frac{1}{I}\sum_{i=1}^I w_i^{(m)}\Bigl[\mu_i^{(m)}\bigl(\mu_i^{(m)}\bigr)^T + V_i^{(m)}\Bigr].
-
残差方差更新：
\sigma^{2(m+1)} = \frac{1}{IJ}\sum_{i=1}^I w_i^{(m)}\Biggl\{\Bigl\|\mathbf{y}_i-A\beta^{(m+1)}-A\,\mu_i^{(m)}\Bigr\|^2 + \operatorname{tr}\Bigl(A\,V_i^{(m)}A^T\Bigr)\Biggr\}
(5) 检查收敛
判断 \theta^{(m+1)} 与 \theta^{(m)} 是否满足收敛标准
TQ9（处理异常值） 和 TQ10（处理缺失数据） 本质都是用 EM 算法填补缺失数据，因此构造出的完整数据对数似然具有相同的形式，因此在构造 Q 函数时，依赖于参数 \theta=(\beta,\Sigma_b,\sigma^2) 的那部分对数似然（关于 \mathbf{y}_i-A\beta-A\,b_i 和关于 b_i 的部分）在形式上是一致的，因此，当我们对 Q 函数求极大化时，无论是处理异常值还是缺失数据，得到的关于 \beta、\Sigma_b 和 \sigma^2 的更新公式都是相同的。
TQ9（处理异常值） 和 TQ10（处理缺失数据） 问题的区别:
TQ9 的 EM 算法：
- 隐藏数据包括模型中的随机效应 b_i 以及尺度混合表示中引入的纹理参数 \tau_i，完整数据为 \{ \mathbf{y}_i,\, b_i,\, \tau_i \}。这主要用于处理异常值，使模型具有鲁棒性。
- E 步需要计算 E(b_i\mid \mathbf{y}_i,\theta^{(m)})、\operatorname{Var}(b_i\mid \mathbf{y}_i,\theta^{(m)}) 以及 E(\tau_i\mid \mathbf{y}_i,\theta^{(m)})
- M 步直接利用 完整数据 构造 Q 函数，对参数进行更新
- 本来就有完整的观测值，只不过有一些内部的隐藏变量（b_i 和 \tau_i），所以我们先直接写下完整数据的似然，然后用 E 步计算这些隐藏变量的期望，再在 M 步中求参数更新
TQ10 的 EM 算法：
- 隐藏数据包括模型中的随机效应 b_i 以及观测向量中的缺失部分 \mathbf{y}_{i,mis} ，完整数据为 \{ \mathbf{y}_{i,obs},\, \mathbf{y}_{i,mis},\, b_i \}，主要用于处理数据缺失问题。
- 除了计算 b_i 的后验均值与协方差之外，还需要对缺失观测值 \mathbf{y}_{i,mis} 求 条件期望（以及相应的条件协方差），然后把缺失值填补
- 观测数据是不完整的，必须先填补缺失的 \mathbf{y}_{i,mis}，才能写成完整数据的对数似然，也就是说比 TQ9 多了填补缺失值这一步
