GLMM（广义线性混合模型）原理以及分析步骤

生态和进化研究中，研究者收集的数据往往 “不规矩”：比如调查某个物种在不同区域 “有或没有”（二元数据）、统计鸟窝中雏鸟存活比例（比例数据）、数一片草地里幼苗的数量（计数数据）。这些数据都不符合基础统计学要求的 “正态分布”，而且还常存在 “随机差异”—— 比如同一种植物在不同地块的生长情况不同（地块的随机效应）、同一棵树结的果实数量每年有波动（时间的随机效应）。

过去处理这些数据的方法都有漏洞：要么强行把数据 “掰成” 正态分布（比如给数据做对数转换，但含大量零值的计数数据根本转不动）；要么假装没有随机差异（导致统计结果不准，出现 “伪重复” 错误）；要么把随机差异当成固定不变的影响（比如把某几个地块的差异当成所有地块的规律，结论没法推广）。而 GLMMs 就是为解决这些问题而生的 —— 它既能处理 “不规矩” 的非正态数据，又能把随机差异纳入统计，堪称这类研究的 “定制化工具”。

GLMM到底是什么？

广义线性混合模型（GLMM） 是一种统计模型，它通过一个链接函数，建立响应变量的数学期望与线性预测变量之间的关系，其中线性预测变量既包含固定效应，也包含随机效应，并且响应变量可以服从指数族分布（如二项分布、泊松分布等）。

可以把它拆解成三个部分来理解，它其实就是三个经典统计模型的“超级合体”：

广义线性模型 (GLM)：处理非正态数据。
线性模型 (LM)：建立线性关系。
混合模型 (MM)：同时包含固定效应和随机效应。

是不是有点晕？别急，我们用一个比喻来彻底理解它。

假设你的任务是研究“学习时间”如何影响“数学考试成绩”。

响应变量 (Y)：数学考试成绩。这是一个计数数据（0到100分），通常不是完美的正态分布。
固定效应 (X)：学习时间。这是你核心关心的因素，你想知道它的影响。
随机效应：学生来自不同的学校。你知道不同学校的教学质量、资源完全不同，这会影响成绩。但你并不关心“A校比B校具体高几分”，你只关心“学校之间差异有多大”。

现在，我们看看不同类型的模型会怎么处理这个问题：

1. 如果用普通线性模型 (LM)

成绩 = β₀ + β₁ * 学习时间 + ε
问题：成绩可能不是正态分布，而且它完全忽略了“学校”的影响。结果会把所有学生混在一起，误差ε会非常大，模型不准。

2. 如果用广义线性模型 (GLM)

g(成绩) = β₀ + β₁ * 学习时间 （比如用泊松分布和log链接函数）
改进：解决了成绩非正态的问题。
问题：依然忽略了“学校”效应。模型假设所有学生都是在同一条件下学习的，这显然不合理。

3. 如果用线性混合模型 (LMM)

成绩 = β₀ + β₁ * 学习时间 + γₖ（学校随机效应） + ε
改进：加入了“学校”作为随机效应（γₖ），承认了学校间的差异。
问题：成绩可能仍然不服从正态分布，模型假设 residuals (ε) 是正态的，这可能不成立。

4. 最终方案：广义线性混合模型 (GLMM) —— 合体！

GLMM 把上面所有的优点结合起来：

g(E(成绩)) = β₀ + β₁ * 学习时间 + γₖ（学校随机效应）

广义 (Generalized)：通过一个链接函数 g()（例如log函数）来处理非正态的考试成绩，使其符合模型要求。
线性 (Linear)：预测部分 β₀ + β₁ * 学习时间 + γₖ 仍然是线性的。
混合 (Mixed)：模型同时包含了：
- 固定效应：β₁ * 学习时间（我们想知道的核心效应）。
- 随机效应：γₖ（代表第k个学校的效应，我们假设它来自一个正态分布，关心的是其方差）。

这样，GLMM就能既处理非正态数据，又考虑到了数据的分层/分组结构，从而得出更准确、更可靠的结论。

一个GLMM就像一台精密的机器，由以下几个关键零件构成：

组件	是什么？	在比喻中代表	目的
响应变量 (Y)	要分析的结果	数学考试成绩	我们想解释的东西
误差分布	Y的数据类型	泊松分布	尊重数据的本质属性
链接函数 (g)	连接Y与预测值的桥梁	log() 函数	让线性预测值落在Y的合理范围内
固定效应	研究者关心的主要因素	学习时间	估计其具体效应大小
随机效应	代表分组结构的变异	学校	估计组间差异（方差），而非具体效应

GLMM完整分析流程

一、明确分析目标

GLMM的核心目标有两个：

算 “固定效应的效应量”：比如 “施肥量每增加 10kg，植物结果数能多多少”“温度升高 1℃，物种存活率会变多少”—— 这些是你研究要回答的核心问题（如 “施肥是否促进植物繁殖”）；
算 “随机效应的方差”：比如 “不同地块间植物生长差异有多大”“不同基因型的后代存活差异有多大”—— 这些是数据中的随机波动，决定你的结论能否推广到更多地块、基因型。

简单说：先想清楚 “我关注的因素（固定效应）有什么影响”“数据中的随机差异（随机效应）有多大”，再开始后续步骤。

二、选择变量，确定模型里包含什么

1. 选 “固定效应”：只留 “和目标相关的关键变量”

选什么：优先纳入 “研究目标直接相关的变量”（如实验中的 “施肥处理”“剪枝处理”，调查中的 “降水”“温度”），以及 “重要的交互作用”（如 “施肥 × 剪枝”—— 判断施肥效果是否受剪枝影响）；
不选什么：别加无关变量（如研究植物结果数，没必要加 “土壤颜色” 这种无生物学意义的变量），避免模型太复杂。

2. 选 “随机效应”：必须满足 “样本量要求”，不然算不出来

选什么：生态与进化研究里最常用的是 “实验区块”（如不同地块）、“个体”（如同一植株的多个后代）、“基因型 / 种群”、“时间周期”（如不同年份）；
关键要求：每个随机效应至少有5-6 个水平（比如 “地块” 作为随机效应，至少选 5 个不同地块），每个 “处理组 + 随机效应水平” 至少有10-20 个样本（比如 “施肥 + 地块 1” 这个组合，至少测 10 株植物）。

3. 选择误差分布与链接函数

这是GLMM处理“非正态数据”的核心。

3.1 误差分布

误差分布描述了你的响应变量（y）天生具有的随机变异模式。

它回答了这个问题： “我的数据是怎么分散开的？”
为什么它重要？ 不同的数据类型，其变异的“规矩”完全不同。你不能用对待身高的方式去对待计数数据。

常见的数据类型及其对应的分布：

数据类型	例子	合适的分布	直观理解
连续数据	身高、体重、浓度、温度	正态分布	数据对称地分布在均值两侧，极端值较少。
计数数据	动物数量、发病次数、叶片数	泊松分布	数据是非负整数（0, 1, 2...），且方差与均值相等。
二元数据	生存/死亡、成功/失败、有/无	二项分布	结果只有两种可能，关心的是其中一种发生的概率（p）。
比例数据	死亡率、出苗率、性别比	二项分布	本质是多次二元试验的成功比例（成功次数/总试验次数）。

选错了会怎样？
如果你对计数数据（比如“每只羊身上的虱子数”）使用正态分布，模型会错误地估计不确定性，可能预测出负的虱子数（这显然不可能），导致结论完全错误。

3.2 链接函数

一句话解释：链接函数是一个数学变换，它把模型的线性预测结果“映射”到响应变量分布所允许的范围内。

它回答了这个问题： “我如何确保模型预测的值符合现实？”
为什么它重要？ 模型的右边（β₀ + β₁x₁ + ...）理论上可以从 -∞ 到 +∞，但你的响应变量（y）可能被限制在一个范围内（比如概率必须在0到1之间）。链接函数就是架在这两者之间的桥。

常见的链接函数：

分布	常用链接函数	作用	直观理解
正态分布	恒等链接	不做任何变换。`预测值 = β₀ + β₁x₁`	线性关系，最直接。
泊松分布	对数链接	`log(预测值) = β₀ + β₁x₁`	确保预测值永远为正数。取指数后变回原尺度。
二项分布	Logit链接	`log( p / (1-p) ) = β₀ + β₁x₁`	将概率p从(0,1)映射到(-∞, +∞)。`p / (1-p)`是发生比（Odds）。

举个例子（Logit链接）：
我们想用温度（x）预测企鹅蛋孵化的概率（p）。

概率p只能在0到1之间，但线性组合 β₀ + β₁*温度 可能算出 -10 或 100，这显然不是概率。
Logit函数 将概率p转换成 “发生比的对数”（也叫Log-Odds）。发生比（Odds）是 (概率p) / (1-概率p)，它的范围是(0, +∞)，取对数后范围就是(-∞, +∞)——这下就能和 β₀ + β₁*温度 完美对接了！
当我们得到模型结果后，再通过逆Logit函数（即逻辑函数）转换回来，就能得到在特定温度下介于0和1之间的预测概率。

二、检查数据是否合格

1.为什么要检查数据？

核心原因：所有的统计模型都建立在一系列假设之上。如果数据严重违背这些假设，你的模型结果（参数估计、p值、置信区间）将是不可靠的，甚至是误导性的。

具体来说，检查是为了：

避免模型失败（GIGO原则）：垃圾进，垃圾出。数据中的问题（如异常值、结构缺失）会导致模型无法收敛（算不出结果）或计算出错误的结果。
验证你的选择：你为模型选择的“零件”（如泊松分布）是否真的适合你的数据？检查可以帮助你确认或质疑最初的选择。
提高统计功效：干净、符合假设的数据能让模型更准确地捕捉到真实存在的效应，从而提高发现显著结果的能力。
确保结论的可解释性：只有当你理解了数据的结构和潜在问题后，你对模型结果的解释才是有意义的。

2.如何检查数据？

图表类型	什么时候用	检验目标（数据 “合格” 的标准）
散点图	固定效应中有 “连续变量”（如施肥量、温度）	看 “经链接函数转换后的因变量” 和连续变量是否呈直线关系—— 比如计数数据经对数转换后，和施肥量的散点要大致排成直线，不能是曲线
箱线图	固定效应中有 “分类变量”（如处理组 / 对照组、不同地块）	看不同分类组间，“转换后的因变量” 的方差是否差不多—— 比如对照组和处理组的箱线图宽度一致，不能一个宽一个窄
残差图	先拟合一个 “只含固定效应的 GLM”（暂时不算随机效应）	看残差（实际值 - 预测值）是否随机分布（没有明显趋势，比如残差不随预测值增大而变大），同时找 “异常值”（如某数据的残差远超其他点，可能是记录错误）

三、估算参数——按照模型复杂度选方法

在GLMM中，我们既要估计固定效应（比如实验处理的影响、连续变量的斜率），也要估计随机效应（比如不同个体、种群、时间点之间的变异）。

固定效应参数（比如处理效应、协变量系数）：描述“平均趋势”。
随机效应参数（其实主要是它们的方差 σ²）：描述“组间差异的大小”。

所以我们写出来的目标函数就是：

最大化它，就能得到固定效应和随机效应的估计。

“最大化似然”意思是：

我假设某个参数值（比如固定效应 β=2，随机效应标准差 σ=1）。
然后问：在这个假设下，我手头的观测数据出现的概率有多大？
如果换一组参数（比如 β=3，σ=0.5），算出来的概率更大，那就说明后者更符合数据。

最终找到那组参数，让“数据出现的可能性”最大。

在普通回归里，似然函数很好算。但在 GLMM 中有 随机效应，它们并不是固定的常数，而是从一个分布（通常是正态分布）里来的。
因此似然函数要把随机效应“积分掉”：

其中 u是随机效应。这个积分往往没有解析解，尤其是随机效应多的时候，直接算不动。

所以常常是使用各种方法近似计算这个积分，以下是方法的总结。

估算方法	优势	劣势	适用场景
PQL（惩罚拟似然）	计算快、软件易实现（SAS PROC GLIMMIX、R glmmPQL）	随机效应标准差大或数据均值小时（如泊松均值 < 5、二元数据）估算偏倚，无法用于似然推断	简单模型（单随机效应）、数据符合 “均值 / 成功数> 5”
拉普拉斯近似	精度高于 PQL，能近似真实似然，支持似然推断	比 PQL 慢、灵活性低	中等复杂度模型（1-2 个随机效应）、需精准结果
GHQ（高斯 - 埃尔米特求积）	精度最高	计算极慢，随机效应 > 2-3 个时不可用	随机效应少、追求极致精度
MCMC（马尔可夫链蒙特卡洛）	灵活，能处理多随机效应，可结合先验信息（贝叶斯框架）	计算慢、技术门槛高（需设先验、验收敛）	复杂模型（多随机效应）、数据量大

另外，算固定效应时用 “ML（最大似然）”，方便比较不同固定效应的模型，但会低估随机效应的标准差（即随机效应的变异大小），仅在 “样本量极大” 时，ML 对随机效应的估算偏倚才会减小，小 / 中等样本下偏倚明显，因此不适合精准估算随机效应参数；算随机效应时用 “REML（限制性最大似然）”，能减少方差估算的偏倚，但无法直接比较 “固定效应不同的模型”—— 因 REML 会将固定效应的信息纳入似然计算，不同固定效应模型的 “剩余似然” 不具有可比性，若强行比较会导致结果错误。

四、验证与优化模型

算完参数不代表结束，最初的 “全模型” 可能包含多余效应（如某随机效应实际无变异），需要验证和优化，避免结论片面：

1. 验证模型假设（过离散+残差异常值）：确保参数估算的可靠性

检验过离散：针对计数数据（泊松分布）或比例数据（二项分布），计算 “Pearson 残差平方和”，若该值远大于 “残差自由度”（样本量 - 模型参数个数），说明数据存在过离散（实际波动比模型预测大），需更换误差分布（如泊松换负二项）或用 “准分布”（准泊松、准二项）修正，修正后需重新估算参数；若残差平方和与残差自由度大致相等，说明无过离散
检验残差异常值与趋势：绘制 “标准化残差 - 拟合值” 散点图，观察残差是否呈 “随机分布”（无明显趋势，如残差不随拟合值增大而增大 / 减小），同时识别极端异常值（残差绝对值远超其他点的观测）—— 若存在异常值，需核查数据记录（如是否为录入错误），必要时删除极端值并重新估算参数；若残差有明显趋势，需调整模型（如增加协变量、更换链接函数）

2.模型选择

4.1 信息理论方法

核心思想：比较多个模型的预期预测能力，而不是检验某个效应是否为0。它允许比较非嵌套模型。
如何操作：
1. 建立一组候选模型（这些模型是全局模型的不同简化版本）。
2. 用最大似然（ML） 而不是REML来拟合所有模型（因为REML不能用于比较固定效应不同的模型）。
3. 为每个模型计算Akaike信息准则（AIC）。AIC = -2 × 对数似然值 + 2 × 模型参数个数。A值越小，模型越好。
4. 计算每个模型与最佳模型（AIC最小）的差值（ΔAIC）。
5. 根据ΔAIC评估模型支持度：
  - ΔAIC < 2：该模型与最佳模型同样好。
  - 4 < ΔAIC < 7：该模型的支持度显著降低。
  - ΔAIC > 10：该模型基本上没有支持。
变体：
- AICc：用于小样本情况。
- QAIC：用于存在过离散的数据（与准似然方法配套使用）。
模型平均（Model Averaging）：如果几个模型的ΔAIC都很接近（例如< 2），说明没有唯一的最佳模型。建议不要只选一个，而是对这些模型的参数估计进行加权平均，从而得到一个更稳健的结论。

4.2 频率假设检验方法

核心思想：使用统计检验来判断某个效应是否“显著”不为零。

1.用 Wald 检验筛选固定效应，判断核心变量的显著性

待模型假设验证通过（无过离散、残差正常）后，再用 Wald 检验筛选固定效应，

需根据 “是否存在过离散” 选择检验形式：无过离散时用Wald Z 检验或 χ² 检验，有过离散（已用准分布修正）时用Wald t 检验或 F 检验（t/F 检验能校正过离散带来的误差不确定性）；
计算固定效应的 “参数估计值 ÷ 标准误” 得到检验统计量，结合残差自由度判断 P 值 —— 若 P<0.05，认为固定效应显著，保留在模型中；若 P≥0.05，可结合生物学意义决定是否删除（避免仅依赖统计结果忽略重要生物学效应）。

2. 用似然比检验筛选随机效应，判断随机变异的必要性

作用：判断随机效应所代表的变异是否重要（例如，个体间的差异是否显著大于0？）。

做法：

- 拟合两个模型：一个包含该随机效应（完整模型），一个不包含（简化模型）。
- 使用anova(model_full, model_reduced)比较两个模型。
- 注意：由于随机效应的方差不能为负，其零假设（方差=0）位于参数空间的边界。因此，传统的*p*值需要调整（通常是除以2），或者使用更复杂的混合χ²分布。
判断标准：校正后若 P<0.05，认为随机效应有必要纳入模型；若 P≥0.05，可删除该随机效应（尤其当随机效应的方差估计值极小时）

缺点：这种方法通常涉及逐步回归（stepwise regression），而统计学家强烈反对盲目使用逐步回归，因为它会放大假阳性率（I类错误），并且模型选择结果很不稳定。

4.3 贝叶斯方法（作为替代）

核心思想：使用MCMC采样，并计算DIC（Deviance Information Criterion） 作为AIC的贝叶斯替代品来进行模型比较。DIC会自动惩罚模型复杂度。
优点：特别适合复杂模型，能很自然地处理不确定性

优先级排序：

若研究目标是 “明确单个效应是否必要”（如 “地块随机效应是否需保留”“施肥固定效应是否显著”），优先用频率学派假设检验（似然比检验筛随机效应 + Wald 检验筛固定效应）；
若研究目标是 “比较多候选模型、注重预测准确性”（如 “施肥 + 温度” 与 “施肥 + 降水” 哪个更能解释植物结果数），优先用信息论方法（AIC/AICc/QAIC + 模型平均）；
若模型复杂（含 3 个及以上随机效应）或有可靠先验信息，优先用贝叶斯方法（MCMC 后验概率 / DIC）。

文献参考

Bolker, Benjamin M., et al. "Generalized linear mixed models: a practical guide for ecology and evolution." Trends in ecology & evolution 24.3 (2009): 127-135.