广义线性混合模型是什么？为什么用？在哪里用？如何工作？如何实现？需要多少数据？

广义线性混合模型（GLMM）是什么？

简单来说，广义线性混合模型（Generalized Linear Mixed Model, GLMM）是一种强大的统计建模工具，它结合了两种重要模型的特点：广义线性模型（Generalized Linear Model, GLM）和线性混合模型（Linear Mixed Model, LMM）。

它被设计用来分析那些既有非正态分布的响应变量（比如计数、比例、二元结果等），同时数据结构中又存在相关性或分组结构（比如重复测量数据、嵌套数据、面板数据等）的情况。

在GLMM中：

• 它继承了GLM处理非正态响应变量的能力，通过使用一个连接函数（Link Function）将响应变量的均值与线性预测器关联起来，并且允许响应变量服从指数族（Exponential Family）中的各种分布（如泊松分布、二项分布、伽马分布、负二项分布等）。
• 它继承了LMM处理数据相关性的能力，通过引入随机效应（Random Effects）来建模或解释这种相关性或分组带来的变异。

因此，GLMM的核心在于它能够同时处理：
非正态响应 + 数据相关性/分组结构

为什么要使用广义线性混合模型？

当你的数据满足以下两个条件：1. 响应变量不是连续且近似正态分布的（例如，是计数、是0/1结果、是比例等）；2. 数据点之间不是完全独立的（例如，同一受试者在不同时间的测量，来自同一家庭的成员，在同一学校的学生，在同一区域采集的样本等），那么你很可能就需要GLMM。

不使用GLMM，而使用更简单的模型（如标准GLM或忽略相关性的模型）会导致什么问题？

• 标准误差（Standard Errors）和P值不正确： 忽略数据内部的相关性会夸大你的样本量，导致标准误差被低估，P值过小，从而增加犯第一类错误（拒绝真实的原假设）的风险。
• 参数估计偏差： 在某些情况下，忽略相关性还会导致固定效应的参数估计产生偏差。
• 无法解释变异来源： 混合模型的一个重要优点是可以分解不同层级的变异（例如，个体内部的变异 vs 个体之间的变异，学校内部的变异 vs 学校之间的变异）。标准模型无法做到这一点。
• 处理复杂数据结构： GLMM可以更灵活地处理不平衡数据（每组测量次数不同）、缺失数据（在混合模型框架下处理更合理）以及更复杂的嵌套或交叉分组结构。

简单来说，使用GLMM是为了在响应变量非正态的情况下，更准确地建模和分析具有内在相关结构的数据，避免因忽略相关性而导致的错误推断。

广义线性混合模型在哪里使用？

GLMM的应用非常广泛，几乎涵盖了所有处理分组数据或重复测量数据，且响应变量是非连续或非正态的领域。以下是一些常见的应用场景：

• 医学和公共卫生：
  • 跟踪患者治疗效果（二元：改善/未改善；计数：发病次数）在不同时间点上的变化（重复测量），考虑患者个体之间的差异（随机效应）。
  • 分析多中心临床试验数据，考虑来自不同中心（分组）的患者之间的差异。
  • 研究疾病在不同家庭成员（分组）中的传播或关联性（二元：患病/未患病）。
• 心理学和行为科学：
  • 分析认知实验中被试在不同条件或不同时间点上的反应（二元：正确/错误；反应时间分布通常非正态）。
  • 调查学生在不同学校或班级（分组）的学习成绩（等级或计数）。
  • 研究家庭治疗对家庭成员互动频率（计数）的影响。
• 生态学和环境科学：
  • 分析不同样地（分组）内物种出现的频率或数量（计数）。
  • 跟踪同一批动物在不同时间点的位置或行为（可能需要特定分布）。
  • 研究污染对不同区域（分组）内植物生长（可能是非正态分布）的影响。
• 农业科学：
  • 评估不同处理方法对同一块地（分组）上作物产量（可能非正态）或病虫害数量（计数）的影响。
  • 跟踪同一批牲畜在不同生长阶段的体重增加（可能需要Gamma分布）。
• 教育研究：
  • 分析学生在不同学校（分组）的考试通过率（比例或二元结果）。
  • 研究教师培训对教师行为（计数或等级）的影响，考虑教师在不同学校的差异。
• 社会科学：
  • 分析调查数据中，来自同一家庭、社区或公司（分组）的个体的态度或行为（等级、二元、计数）。
  • 研究政策变化对不同地区（分组）失业率（比例）的影响。

只要你的数据有分组或重复测量，并且响应变量不是通常的正态连续变量，GLMM就可能是合适的选择。

广义线性混合模型如何工作？

GLMM的工作原理比标准线性模型或标准GLM更复杂，因为它需要同时处理连接函数、非正态分布和随机效应。其核心思想是将线性预测器分解为固定效应和随机效应两部分，然后通过连接函数将这个完整的线性预测器映射到响应变量的均值。

模型的构成

一个典型的GLMM结构可以概括为：

1. 线性预测器（Linear Predictor）： 这部分包含了固定效应和随机效应。
   • 固定效应 (Fixed Effects)：
     这些是你感兴趣的主要自变量，它们对响应变量的影响是固定的，不随个体或组变化。例如，治疗方式、性别、年龄等。用矩阵形式表示为 $X\beta$，其中 $X$ 是固定效应的设计矩阵，$\beta$ 是对应的参数向量。
   • 随机效应 (Random Effects)：
     这些用来建模由于数据结构（分组、重复测量）引入的变异性或相关性。随机效应是随机变量，通常假设它们服从均值为零的多元正态分布，具有一个需要估计的协方差矩阵 $G$。例如，个体之间的随机差异、不同学校之间的随机差异。用矩阵形式表示为 $Z\gamma$，其中 $Z$ 是随机效应的设计矩阵，$\gamma$ 是对应的随机效应向量 ($\gamma \sim N(0, G)$)。

   完整的线性预测器通常表示为 $\eta = X\beta + Z\gamma$。注意，随机效应 $\gamma$ 依赖于你定义的分组结构。

2. 连接函数（Link Function）：
   一个可逆函数 $g$，它将响应变量的期望值（均值）$\mu$ 与线性预测器 $\eta$ 联系起来：$g(\mu) = \eta$。选择哪个连接函数取决于响应变量的分布。例如，对于泊松回归（计数），通常用对数连接函数 $g(\mu) = \log(\mu)$；对于二元回归，通常用logit连接函数 $g(\mu) = \text{logit}(\mu) = \log(\frac{\mu}{1-\mu})$。
3. 条件分布（Conditional Distribution）：
   在给定随机效应 $\gamma$ 的条件下，响应变量 $Y$ 服从一个特定的指数族分布，其均值由 $\mu = g^{-1}(\eta) = g^{-1}(X\beta + Z\gamma)$ 决定。这个分布是针对单个观测的，假设在条件于随机效应的情况下，不同观测是独立的。例如，对于二元数据可能是 Bernoulli 分布，对于计数数据可能是 Poisson 或 Negative Binomial 分布。

模型估计的挑战与方法

与线性混合模型（LMM）不同，GLMM的条件分布通常不是正态分布，这使得计算其边际似然（通过对随机效应积分得到）变得非常复杂，往往没有解析解。积分的形式如下：

$L(\beta, G | Y) = \int L(\beta | Y, \gamma) f(\gamma | G) d\gamma$

其中 $L(\beta | Y, \gamma)$ 是在给定随机效应下的条件似然，$f(\gamma | G)$ 是随机效应的概率密度函数。

因此，GLMM的参数估计（寻找使边际似然最大或后验概率最大的参数值）通常需要依赖于数值近似方法或模拟方法：

• 伪似然方法（Pseudo-Likelihood, PL 或 PQL）： 这是早期使用较多的方法，通过对响应变量或其函数进行泰勒展开来近似处理非线性。虽然计算速度快，但在离散分布（如二元、泊松）且随机效应方差较大时，可能产生显著的偏差。
• 拉普拉斯近似（Laplace Approximation）： 对积分进行拉普拉斯近似，比PQL更精确，但仍然是近似方法。
• 自适应高斯-赫米特求积（Adaptive Gaussian-Hermite Quadrature, AGQ）： 一种更精确的数值积分方法，通过在每个随机效应的条件众数周围进行积分。对于随机效应维度较低（比如只有1-2个随机效应变量）时效果很好，但当随机效应维度增加时计算量呈指数增长。
• 最大似然（Maximum Likelihood, ML）和限制最大似然（Restricted Maximum Likelihood, REML）： 对于某些特定的GLMM（如泊松对数线性混合模型），可以推导出精确或近似的ML/REML方法。但通常计算复杂。
• 马尔科夫链蒙特卡洛（Markov Chain Monte Carlo, MCMC）： 基于模拟的方法，通过构建一个马尔科夫链来从后验分布中采样，从而估计参数。计算量大，但可以提供更完整的参数分布信息，尤其适用于复杂模型。

现代统计软件通常提供了多种估计算法选项，用户需要根据数据和模型特点进行选择，并在可能的情况下对结果进行验证。

如何使用/实现广义线性混合模型？

实现GLMM通常需要专门的统计软件和相应的函数库。整个流程大致包括数据准备、模型设定、模型拟合、诊断与评估以及结果解释。

1. 数据准备与探索

• 确保你的数据包含了响应变量、固定效应变量以及至少一个标识分组或重复测量个体的变量（用于定义随机效应）。
• 对数据进行描述性统计和可视化，了解响应变量的分布特征，以及各组/个体内部和之间的变异性。这有助于选择合适的响应变量分布和随机效应结构。

2. 模型设定

这是GLMM建模的核心，你需要明确指定以下几个部分：

• 响应变量： 你要建模的因变量。
• 响应变量分布（Distribution）： 根据响应变量的类型和分布特点选择。例如：
  • 二元（0/1）：Bernoulli 或 Binomial（如果是N次试验中的成功次数）。
  • 计数（非负整数）：Poisson 或 Negative Binomial（如果存在过度分散）。
  • 比例（0到1之间）：Binomial（需要指定试验总次数）或 Beta Binomial，有时也用 Beta 分布（如果不是基于固定次数的计数）。
  • 连续非负数据（如时间、金额）：Gamma 或 Inverse Gaussian。
• 连接函数（Link Function）： 与所选分布相匹配的标准连接函数，或根据需要选择其他合适的函数。
  • Bernoulli/Binomial: Logit (最常用), Probit, Cloglog
  • Poisson/Negative Binomial: Log
  • Gamma/Inverse Gaussian: Log, Inverse
• 固定效应（Fixed Effects）： 包含在模型线性预测器中的自变量，它们对响应变量的影响是固定不变的。类似于标准回归模型中的预测变量。
• 随机效应（Random Effects）： 建模分组或重复测量带来的变异性。常见的随机效应结构包括：
  • 随机截距 (Random Intercept)： 假设每个组/个体有一个不同的基线水平，但所有组/个体对固定效应变量的反应斜率是相同的。例如，`(1 | GroupID)` 表示 GroupID 标识的每个组有一个随机截距。
  • 随机斜率 (Random Slope)： 假设不同组/个体对某个或多个固定效应变量的反应斜率是不同的。例如，`(Variable | GroupID)` 表示 Variable 对响应变量的影响斜率在 GroupID 各组之间是随机变化的。`(1 + Variable | GroupID)` 表示同时有随机截距和随机斜率，且它们之间可能存在协方差。
  • 更复杂的结构： 如嵌套随机效应（学生嵌套在班级中，班级嵌套在学校中）或交叉随机效应（同一个评估者评估不同的被试，同一个被试被不同的评估者评估）。

3. 模型拟合

使用统计软件中的GLMM函数来拟合模型。你需要指定数据、公式（包含响应变量、固定效应和随机效应）以及响应变量的分布和连接函数。

例如，在 R 语言中使用 lme4 或 glmmTMB 包：

library(lme4)
model <- glmer(Response ~ FixedVar1 + FixedVar2 + (1 | GroupID), data = mydata, family = poisson(link = "log"))

或者使用 glmmTMB 包，它支持更多的分布和结构：

library(glmmTMB)
model <- glmmTMB(Response ~ FixedVar1 + FixedVar2 + (1 | GroupID), data = mydata, family = poisson(link = "log"))

在 SAS 中使用 PROC GLIMMIX，在 Stata 中使用 gllamm 或 mixed 命令，Python 的 statsmodels 库也在不断完善这方面的功能。

4. 模型诊断与评估

拟合模型后，进行诊断以检查模型的假设是否得到满足，并评估模型的拟合优度：

• 残差分析： 检查残差图，看是否存在模式、异方差或异常值。对于GLMM，残差定义和解释可能比线性模型更复杂。
• 随机效应分布： 检查估计的随机效应是否近似服从正态分布（这是一个常见的假设）。
• 模型拟合优度： 使用信息准则（如 AIC, BIC）比较不同模型；对于嵌套模型，可以使用似然比检验（LRT），但要注意自由度的计算。对于非计数数据，也可以使用残差图等。
• 过分散检查： 如果响应变量是计数或二项数据，检查是否存在过度分散（Observed variance > Theoretical variance for Poisson/Binomial），如果存在，可能需要改用负二项分布或 Beta Binomial 分布。

5. 结果解释

解释模型拟合结果：

• 固定效应： 解释固定效应系数的含义。记住它们是在连接函数的尺度上解释的（例如，泊松回归中，系数表示对数平均计数的变化；二元回归中，系数表示 log-odds 的变化）。通常需要转换回原始尺度（例如，计算赔率比 Odds Ratio 或率比 Rate Ratio）以便于理解。
• 随机效应： 解释随机效应的方差分量。这量化了不同分组或个体带来的变异程度。例如，随机截距的方差表示个体基线水平的异质性；随机斜率的方差表示个体对某个预测变量反应斜率的异质性。如果随机截距和斜率的协方差也被估计，解释它们的关联。
• 组内相关性（Intraclass Correlation Coefficient, ICC）： 对于随机截距模型，可以计算 ICC 来表示总变异中由分组结构解释的比例。

使用广义线性混合模型需要多少数据？

这是一个常见且关键的问题，但没有一个简单的数字答案。所需数据量取决于模型的复杂性、数据的变异性以及你想要检测到的效应大小。然而，对于GLMM，一个特别重要的考虑因素是分组/簇（Groups/Clusters）的数量，而不仅仅是总的观测数量。

关键因素：分组数量

估计随机效应的方差成分需要足够多的分组。如果只有很少量的分组，模型很难准确估计组间的变异性，拟合可能不稳定，甚至无法收敛。

• 经验法则： 虽然没有严格的阈值，但通常建议至少有 10-20 个分组。拥有 30-50 个或更多分组通常可以提供更稳定的随机效应方差估计。如果分组数量非常少（例如，少于5-10个），GLMM可能不适用，或者随机效应的方差估计会非常不可靠。
• 组内观测数量： 每组内的观测数量也很重要。如果每组只有很少的观测，即使分组数量多，模型也可能难以捕捉组内模式或区分组内/组间变异。但相对于分组数量，组内观测数量的要求通常弹性更大。

模型复杂性

模型的复杂性也会影响所需的数据量：

• 随机效应结构： 如果模型包含多个随机效应（例如，随机截距加随机斜率，或多个层级的随机效应），或者随机效应之间允许协方差，那么需要更多的数据（尤其是更多分组）来稳定地估计这些额外的方差和协方差参数。
• 固定效应数量： 像所有回归模型一样，固定效应越多，通常需要越大的样本量来精确估计这些效应。

响应变量分布

对于某些非正态分布（尤其是二元或泊松分布），当事件发生率非常低或非常高时，模型拟合可能更具挑战性，需要更多数据（或更多分组）来确保足够的事件发生数。

总结关于数据量的考量：

1. 最重要：保证有足够多的分组/簇。 这是GLMM成功拟合和可靠推断的基础。
2. 拥有适度的组内观测数量。
3. 总观测数量越多越好，但其重要性次于分组数量（在一定范围内）。
4. 复杂的模型需要更多数据，尤其是更多分组。
5. 罕见事件（如二元数据中的低发生率）需要更大的样本量。

在实际应用中，如果分组数量较少，需要对随机效应的估计持谨慎态度，或者考虑将随机效应简化为固定效应（如果合理且可行），但这会改变模型的解释和推断范围。最好的做法是在设计研究时就考虑到数据结构，并确保采集足够多的分组数据。对于复杂的研究设计，可能需要进行模拟来评估所需样本量和分组数量。