在基本回归分析模型中， 假定回归误差项独立同分布， 方差相等， 还经常假定误差项服从正态分布。 在实际应用回归分析建模时， 还经常遇到模型误差项方差不相等或者误差项之间不独立的情形。 比如， 如果每个观测来自于一个群体的平均值， 设群体中的个体方差相同， 则每个观测的误差方差正比于群体中个体的个数的倒数， 不等于常数。 又比如， 医学研究中每个受试者有多次随访的数值型观测值， 则每个受试者的各次测量值是相关的， 不同受试者之间可以认为是独立的， 这样模型误差不是相互独立的。

考虑如下的线性模型： \[\begin{align} \boldsymbol Y =& X \boldsymbol \beta + \boldsymbol\varepsilon, \\ \boldsymbol\varepsilon \sim& \text{N}(\boldsymbol 0, \sigma^2 W^{-1}) . \tag{35.1} \end{align}\] 其中\(X\)为\(n\times p\)非随机的自变量矩阵， \(\boldsymbol\beta\)是\(p\)维回归系数向量， \(\boldsymbol\varepsilon\)为\(n\)维随机误差向量， \(W\)已知，\(\sigma^2\)未知。

因为标准的线性模型要求\(W = I\)， 所以这个模型是线性模型的推广。

考虑这个模型的估计问题。 设存在矩阵\(A\)使得\(W^{-1} = A A^T\)， 则 \[ A^{-1} \boldsymbol Y = A^{-1} X \boldsymbol \beta + A^{-1} \boldsymbol\varepsilon, \] 令\(\boldsymbol Y^* = A^{-1} \boldsymbol Y\), \(X^* = A^{-1} X\), \(\boldsymbol\varepsilon^* = A^{-1} \boldsymbol\varepsilon\)， 则 \[ \boldsymbol Y^* = X^* \boldsymbol\beta + \boldsymbol\varepsilon^*, \] 其中\(\text{Var}(\boldsymbol\varepsilon^*) = \sigma^2 I\)。 在矩阵可逆条件下\(\boldsymbol\beta\)的最小二乘估计为 \[\begin{aligned} \hat{\boldsymbol\beta} = (X^{*T} X^*)^{-1} X^{*T} \boldsymbol Y^* = (X^T W X)^{-1} X^T W \boldsymbol Y . \end{aligned}\] 这个公式称为加权最小二乘公式。

模型可以进一步推广。 考虑 \[\begin{align} \boldsymbol Y =& X \boldsymbol \beta + Z \boldsymbol \alpha + \boldsymbol\varepsilon, \\ \boldsymbol\varepsilon \sim& \text{N}(\boldsymbol 0, \sigma^2 W^{-1}), \\ \boldsymbol\alpha \sim& \text{N}(\boldsymbol 0, \Sigma) . \tag{35.2} \end{align}\] 其中\(\boldsymbol Y\)为可观测的\(n\)维因变量向量， \(X_{n\times p}\), \(Z_{n \times q}\)已知， \(\boldsymbol\beta_{p\times 1}\)为未知的非随机回归系数向量， 称为固定效应； \(\boldsymbol\alpha_{q\times 1}\)为未知的随机向量， 称为随机效应， \(\Sigma\)未知； \(\boldsymbol\varepsilon\)为未知的随机误差向量， \(\sigma^2\)未知，\(W\)已知。 称此模型为线性混合模型。

在估计模型时， \(\Sigma\)必须是对称非负定阵， 一般有一定的参数结构，所以一般会分解为 \[ \Sigma = \sigma^2 \Lambda_{\theta} \Lambda_{\theta}^T, \] 其中\(\theta\)为待估参数。

因为随机效应\(\boldsymbol\alpha\)未知但一般不需要估计， 所以可以将其作用混合到因变量\(\boldsymbol Y\)的方差结构中， 得到 \[\begin{aligned} E \boldsymbol Y =& X \boldsymbol\beta, \\ \text{Var}(\boldsymbol Y) =& \sigma^2 W^{-1} + Z \Sigma Z^T = \sigma^2(W^{-1} + (Z \Lambda_{\theta})(Z \Lambda_{\theta})^T) . \end{aligned}\]

线性混合效应模型适用的数据包括：

本章理论和例子参考了(Andrzej Galecki 2013)。

本章主要使用nmleU包的armd数据。 ARMD（年龄相关性黄斑变性）是一种老年人眼科疾病， 会使得患病的人逐渐失明。 为了研究某种新药的疗效， 收录了240位病人， 随机分为治疗组与对照组， 并在入组时以及4、12、24、52周时对每位病人测量视力指标。 这样的数据称为纵向数据， 每位病人都测量了5个不同时间的因变量值， 但这5个时间是共同的。 有些病人的随访有缺失， 并且有些是在中间缺失的。 每个病人的各次视力指标都与入组时的基础测量值有关。

加权回归也可以用nlme包的gls()函数计算， gls()可以进行误差项方差可变或者相关的广义最小二乘估计。 如：

loads和slippage的关系存在非线性， 我们用二次项来表示非线性。 模型为： \[\begin{aligned} \text{loads}_{ij} =& \beta_0 + b_i + \beta_1 \text{slippage} + \beta_2 \text{slippage}^2 + \epsilon_{ij}, \\ & i=1,2,\dots, 8, \ j=1,2,\dots, 15. \end{aligned}\] 其中\(\epsilon_{ij}\)独立同\(\text{N}(0, \sigma_e^2)\)分布， \(b_i\)独立同\(\text{N}(0, \sigma_b^2)\)分布， \(\{ b_i \}\)与\(\{ \epsilon_{ij} \}\)相互独立。

于是， \[\begin{aligned} \text{loads}_{ij} | b_i \sim& \text{N}(\beta_0 + b_i, \sigma_e^2), \\ \text{loads}_{ij} \sim& \text{N}(\beta_0, \sigma_b^2 + \sigma_e^2), \\ \text{Cov}(\text{loads}_{ij}, \text{loads}_{ik}) =& \sigma_b^2, \ j \neq k , \\ \text{Corr}(\text{loads}_{ij}, \text{loads}_{ik}) =& \frac{\sigma_b^2}{\sigma_b^2 + \sigma_e^2}, \\ \text{Cov}(\text{loads}_{ij}, \text{loads}_{mk}) =& 0, \ i \neq m . \end{aligned}\]

即同一板材样品针对不同滑动程度的力测量值之间的相关系数为常数， 不同样品之间独立， 称这样的方差结构为复合对称(compound symmetry)。 这样的方差结构比较简单， 在纵向数据情形下不够准确， 因为时间相近的同样品测量值的相关系数更高。

\(\beta_0\)是slippage等于0时loads的基准值， 应该等于0，但是因为每个样品在slippage等于0时的值是\(\beta_0 + b_i\)， 所以我们不限制\(\beta_0 = 0\)。 \(b_i\)是第\(i\)样品与其它样品所需力的差异的表现， 我们不关心\(b_i\)的具体值。

用nlme包建模如下：

程序中random =是随机效应的设置， | specimen表明按照样品分组， 各组之间独立， 每组内部有随机效应， ~ 1表示随机效应中仅有一个随机的截距项(\(b_i\))。

输出中(Intercept), slippage和I(slippage^2)部分是固定效应估计， 与线性回归模型中类似项的解释相同。

输出结果中随机效应部分， 给出了随机截距项的方差估计\(\sigma_b^2 = 0.9467^2\), 随机误差的方差估计\(\sigma_e^2 = 0.5334236^2\)。

计算模型预测值并作图：

进一步地， 考虑每个样品的slippage一次项系数也增加一个随机效应， 二次项系数仍仅有固定效应， 模型变成：

\[\begin{aligned} \text{loads}_{ij} =& \beta_0 + b_i + (\beta_1 + c_i) \text{slippage} + \beta_2 \text{slippage}^2 + \epsilon_{ij}, \\ & i=1,2,\dots, 8, \ j=1,2,\dots, 15. \end{aligned}\] 其中\(\epsilon_{ij}\)独立同\(\text{N}(0, \sigma_e^2)\)分布， \(b_i\)独立同\(\text{N}(0, \sigma_b^2)\)分布， \(c_i\)独立同\(\text{N}(0, \sigma_c^2)\)分布， \(\{ b_i, c_i \}\)与\(\{ \epsilon_{ij} \}\)相互独立。

建模程序：