泊松广义线性混合模型(GLMMs)：lme4与glmmADMB结果对比及模型选择

随机系数泊松模型拟合难度较高，lme4和glmmADMB两个工具输出的参数估计结果往往存在一定差异，本次案例如下：

# 加载依赖包
library(lme4)
library(glmmADMB)

# 读取数据集
myds <- read.csv("https://raw.githubusercontent.com/Leprechault/trash/main/my_glmm_dataset.csv")
str(myds)

# 'data.frame': 526 obs. of  10 variables:
#  $ Bioma      : chr  "Pampa" "Pampa" "Pampa" "Pampa" ...
#  $ estacao    : chr  "verao" "verao" "verao" "verao" ...
#  $ ciclo      : chr  "1°" "1°" "1°" "1°" ...
#  $ Hour       : int  22 23 0 1 2 3 4 5 6 7 ...
#  $ anthill    : num  23.5 23.5 23.5 23.5 23.5 ...
#  $ formigueiro: int  2 2 2 2 2 2 2 2 2 2 ...
#  $ ladenant   : int  34 39 29 25 20 31 16 28 21 12 ...
#  $ unladen    : int  271 258 298 317 316 253 185 182 116 165 ...
#  $ UR         : num  65.7 69 71.3 75.8 78.1 ...
#  $ temp       : num  24.3 24.3 24 23.7 23.1 ...

本次研究的响应变量为携带负载的蚂蚁数量ladenant，自变量为生物群系Bioma、温度temp、湿度UR，蚁穴formigueiro为伪重复项，因此分别采用lme4和glmmADMB构建泊松广义线性混合模型：

lme4拟合结果 #

m.laden.1  <- glmer(ladenant ~ Bioma + poly(temp,2) + UR + (1 | formigueiro), data = myds, family = poisson(link = "log"))
summary(m.laden.1)
# Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
#  Family: poisson  ( log )
# Formula: ladenant ~ Bioma + poly(temp, 2) + UR + (1 | formigueiro)
#    Data: myds

#      AIC      BIC   logLik deviance df.resid 
#  21585.9  21615.8 -10786.0  21571.9      519 

# Scaled residuals: 
#     Min      1Q  Median      3Q     Max 
# -10.607  -4.245  -1.976   2.906  38.242 

# Random effects:
#  Groups      Name        Variance Std.Dev.
#  formigueiro (Intercept) 0.02049  0.1432  
# Number of obs: 526, groups:  formigueiro, 5

# Fixed effects:
#                   Estimate Std. Error z value Pr(>|z|)    
# (Intercept)      0.7379495  0.0976701   7.556 4.17e-14 ***
# BiomaTransition  1.3978383  0.0209623  66.684  < 2e-16 ***
# BiomaPampa      -0.1256759  0.0268164  -4.687 2.78e-06 ***
# poly(temp, 2)1   7.1035195  0.2079550  34.159  < 2e-16 ***
# poly(temp, 2)2  -7.2900687  0.2629908 -27.720  < 2e-16 ***
# UR               0.0302810  0.0008029  37.717  < 2e-16 ***
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

# Correlation of Fixed Effects:
#             (Intr) BmTrns BimPmp p(,2)1 p(,2)2
# BiomaTrnstn -0.586                            
# BiomaPampa  -0.199  0.352                     
# ply(tmp,2)1 -0.208  0.267  0.312              
# ply(tmp,2)2 -0.191  0.085 -0.175 -0.039       
# UR          -0.746  0.709  0.188  0.230  0.316
# optimizer (Nelder_Mead) convergence code: 0 (OK)
# Model is nearly unidentifiable: very large eigenvalue
#  - Rescale variables?

glmmADMB拟合结果 #

m.laden.2 <- glmmadmb(ladenant ~ Bioma + poly(temp,2) + UR + (1 | formigueiro), data = myds, family = "poisson", link = "log")
summary(m.laden.2)
# Call:
# glmmadmb(formula = ladenant ~ Bioma + poly(temp, 2) + UR + (1 | 
#     formigueiro), data = myds, family = "poisson", link = "log")

# AIC: 12033.9 

# Coefficients:
#                 Estimate Std. Error z value Pr(>|z|)    
# (Intercept)      1.52390    0.26923    5.66  1.5e-08 ***
# BiomaTransition  0.23967    0.08878    2.70   0.0069 ** 
# BiomaPampa       0.09680    0.05198    1.86   0.0626 .  
# poly(temp, 2)1  -0.38754    0.55678   -0.70   0.4864    
# poly(temp, 2)2  -1.16028    0.39608   -2.93   0.0034 ** 
# UR               0.01560    0.00261    5.97  2.4e-09 ***
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

# Number of observations: total=526, formigueiro=5 
# Random effect variance(s):
# Group=formigueiro
#             Variance StdDev
# (Intercept)  0.07497 0.2738

结果比对与模型选择方法 #

两个工具拟合得到的模型中，Bioma变量的显著性水平存在极大差异，除直接对比两个包的输出结果外，还可以通过以下方法比对结果、选择最优模型：

• 先排查基础输入与数值稳定性问题：首先确认两个模型调用的数据集、变量编码完全一致，避免输入数据差异导致的结果偏差。另外lme4输出已经提示模型接近不可识别，建议先对温度、湿度等连续自变量做标准化处理，消除量纲导致的数值计算问题后重新拟合两个模型，再做对比。
• 验证模型假设是否成立：泊松模型的核心假设是方差等于均值，你可以计算离散参数：dispersion = sum(resid(m.laden.1, type = "pearson")^2) / df.residual(m.laden.1)，如果离散参数远大于1，说明存在过离散，普通泊松模型本身就不适用，两个工具的估计都会存在偏差，建议换成负二项混合模型重新拟合后再对比。
• 对比拟合优度与实际匹配度：统一口径计算两个模型的AIC、BIC、对数似然值，优先选择指标更优的模型。另外可以做后验预测检查：用拟合模型生成模拟的响应变量数据，和真实数据的分布（均值、分位数、极端值分布）做对比，哪个模型生成的模拟数据和真实数据匹配度更高，适用性更强。
• 验证估计稳健性：用bootstrap重抽样方法对数据集多次有放回抽样，分别用两个工具拟合模型，看参数估计的分布是否稳定，优先选择参数估计置信区间覆盖合理、波动更小的工具结果。另外本次随机效应分组只有5个蚁穴，组数较少会放大不同积分方法的差异，可以用glmmTMB作为第三方参考工具，使用自适应高斯-埃尔米特积分并增加积分点数重新拟合模型，看结果和哪个工具的输出更接近。
• 结合研究逻辑验证：可以先绘制分Bioma的ladenant和温度、湿度的散点拟合线，看原始数据的趋势和哪个模型的参数估计方向、效应大小更匹配，优先选择符合实际数据规律和研究常识的结果。

内容的提问来源于stack exchange，提问作者Leprechault