线性混合效应模型

线性混合效应模型用于涉及依赖数据的回归分析。当处理纵向研究设计或其他涉及对每个受试者进行多次观察的研究设计时，会出现此类数据。一些特定的线性混合效应模型包括

  • 随机截距模型，其中组内的所有响应都通过一个特定于该组的值进行加性偏移。

  • 随机斜率模型，其中组内的响应遵循一个（条件）均值轨迹，该轨迹在观测到的协变量中线性变化，斜率（可能还有截距）随组而变化。

  • 方差成分模型，其中一个或多个分类协变量的水平与从分布中抽取的值相关联。这些随机项根据每个观测值的协变量值加性确定其条件均值。

statsmodels 实现的 LME 主要是基于组的，这意味着随机效应必须为不同组中的响应独立实现。在我们的混合模型实现中有两种随机效应：(i) 具有未知协方差矩阵的随机系数（可能是向量），以及 (ii) 从共同的一元分布中独立抽取的随机系数。对于 (i) 和 (ii)，随机效应通过与组特定设计矩阵的矩阵/向量乘积影响组的条件均值。

一个简单的随机系数示例，如上文（i）所述，是：

\[Y_{ij} = \beta_0 + \beta_1X_{ij} + \gamma_{0i} + \gamma_{1i}X_{ij} + \epsilon_{ij}\]

在这里，\(Y_{ij}\) 是受试者 \(i\) 的第 \(j^\rm{th}\) 次测量的响应，而 \(X_{ij}\) 是该响应的协变量。 “固定效应参数” \(\beta_0\) 和 \(\beta_1\) 由所有受试者共享，误差 \(\epsilon_{ij}\) 与其他所有因素独立，并且同分布（均值为零）。 “随机效应参数” \(\gamma_{0i}\) 和 \(\gamma_{1i}\) 遵循均值为零的二元分布，由三个参数描述：\({\rm var}(\gamma_{0i})\)、\({\rm var}(\gamma_{1i})\) 和 \({\rm cov}(\gamma_{0i}, \gamma_{1i})\)。 还有一个参数用于 \({\rm var}(\epsilon_{ij})\)。

一个简单的方差分量示例，如上文(ii)所述，是：

\[Y_{ijk} = \beta_0 + \eta_{1i} + \eta_{2j} + \epsilon_{ijk}\]

在这里，\(Y_{ijk}\) 是在条件 \(i, j\) 下的第 \(k^\rm{th}\) 次测量的响应。唯一的“均值结构参数”是 \(\beta_0\)。\(\eta_{1i}\) 是独立同分布的，均值为零，方差为 \(\tau_1^2\)，而 \(\eta_{2j}\) 也是独立同分布的，均值为零，方差为 \(\tau_2^2\)。

statsmodels MixedLM 处理大多数非交叉随机效应模型，以及一些交叉模型。要在模型中包含交叉随机效应，需要将整个数据集视为一个单一的组。然后可以使用模型中的方差分量参数来定义具有各种交叉和非交叉随机效应组合的模型。

statsmodels LME 框架目前支持通过 Wald 检验和系数上的置信区间进行事后推断，轮廓似然分析，似然比检验，以及 AIC。

示例

In [1]: import statsmodels.api as sm

In [2]: import statsmodels.formula.api as smf

In [3]: data = sm.datasets.get_rdataset("dietox", "geepack").data

In [4]: md = smf.mixedlm("Weight ~ Time", data, groups=data["Pig"])

In [5]: mdf = md.fit()

In [6]: print(mdf.summary())
         Mixed Linear Model Regression Results
========================================================
Model:            MixedLM Dependent Variable: Weight    
No. Observations: 861     Method:             REML      
No. Groups:       72      Scale:              11.3669   
Min. group size:  11      Log-Likelihood:     -2404.7753
Max. group size:  12      Converged:          Yes       
Mean group size:  12.0                                  
--------------------------------------------------------
             Coef.  Std.Err.    z    P>|z| [0.025 0.975]
--------------------------------------------------------
Intercept    15.724    0.788  19.952 0.000 14.179 17.268
Time          6.943    0.033 207.939 0.000  6.877  7.008
Group Var    40.394    2.149                            
========================================================

详细的示例可以在这里找到

Wiki上有一些笔记本示例： 混合线性模型的Wiki笔记本

技术文档

数据被划分为不相交的组。 第\(i\)组的概率模型为：

\[Y = X\beta + Z\gamma + Q_1\eta_1 + \cdots + Q_k\eta_k + \epsilon\]

哪里

  • \(n_i\) 是组 \(i\) 中的观测数量

  • \(Y\) 是一个 \(n_i\) 维的响应向量

  • \(X\) 是一个 \(n_i * k_{fe}\) 维的固定效应系数矩阵

  • \(\beta\) 是一个 \(k_{fe}\) 维的固定效应斜率向量

  • \(Z\) 是一个 \(n_i * k_{re}\) 维的随机效应系数矩阵

  • \(\gamma\) 是一个 \(k_{re}\) 维的随机向量，均值为0，协方差矩阵为 \(\Psi\)；注意，每个组都有其独立的gamma实现。

  • \(Q_j\) 是一个 \(n_i \times q_j\) 维的设计矩阵，用于第 \(j^\rm{th}\) 个方差分量。

  • \(\eta_j\) 是一个 \(q_j\) 维的随机向量，包含独立且同分布的值，方差为 \(\tau_j^2\)。

  • \(\epsilon\) 是一个 \(n_i\) 维的独立同分布正态误差向量，均值为0，方差为 \(\sigma^2\)；\(\epsilon\) 值在组内和组间都是独立的

\(Y, X, \{Q_j\}\) 和 \(Z\) 必须完全被观测到。\(\beta\)、 \(\Psi\) 和 \(\sigma^2\) 使用 ML 或 REML 估计进行估计， 而 \(\gamma\)、\(\{\eta_j\}\) 和 \(\epsilon\) 是 随机的，因此定义了概率模型。

边际均值结构为 \(E[Y|X,Z] = X*\beta\)。如果仅对边际均值结构感兴趣，GEE 是混合模型的一个良好替代方案。

符号：

  • \(cov_{re}\) 是随机效应协方差矩阵（如上所述为 \(\Psi\)），\(scale\) 是（标量）误差方差。每个方差分量也有一个单一的估计方差参数 \(\tau_j^2\)。对于单个组，给定外生变量的内生变量的边际协方差矩阵为 \(scale*I + Z * cov_{re} * Z\)，其中 \(Z\) 是单个组中随机效应的设计矩阵。

参考文献

实现细节的主要参考资料是：

  • MJ Lindstrom, DM Bates (1988). 用于重复测量数据的线性混合效应模型的Newton Raphson和EM算法。 美国统计协会杂志。第83卷，第404期，第1014-1022页。

请参阅这份更新的文档：

所有似然、梯度和Hessian计算都严格遵循Lindstrom和Bates的方法。

以下两份文档更多是从用户的角度撰写的：

模块参考

模型类是：

MixedLM(endog, exog, groups[, exog_re, ...])

线性混合效应模型

结果类是：

MixedLMResults(model, params, cov_params)

包含拟合线性混合效应模型结果的类。
Last update: Oct 16, 2024