如何从Gamma分布GLMER模型提取参数计算未观测数据的对数似然？

问题核心 #

使用lme4拟合Gamma分布的广义线性混合效应模型（GLMER）后，需要计算未观测数据的对数似然，但无法正确提取Gamma分布的形状（shape）和尺度（scale）参数，原代码得到的参数与真实值差距较大。

原代码的关键问题 #

1. 尺度参数计算错误：用训练集的预测均值计算scale，但测试集的scale应基于对应测试样本的预测均值，而非训练集数据。
2. 参数逻辑混淆：未明确Gamma分布在GLMER中的参数化关系，导致参数提取和计算的对应性出错。

修正方案与代码实现 #

关键参数关系明确 #

在lme4的Gamma族模型中：

• 响应变量的方差满足 Var(y) = μ² * φ，其中μ是预测均值，φ是分散参数。
• 对应Gamma分布的方差公式 Var(y) = μ² / shape，因此可得 shape = 1/φ。
• Gamma分布的均值与参数关系为 μ = shape * scale，因此每个样本的 scale = μ / shape。

修正后的完整代码 #

library(lme4)

set.seed(321)

# 构建训练集
df_train <- data.frame(
  y = rgamma(100, shape = 3, scale = 2),
  x = rnorm(100),
  group = gl(10, 10)
)

# 构建测试集
df_test <- data.frame(
  y = rgamma(50, shape = 3, scale = 2),
  x = rnorm(50),
  group = sample(gl(10, 5), 50, replace = TRUE)
)

# 拟合Gamma GLMER模型（恒等链接）
my_model <- glmer(y ~ x + (1|group), family = Gamma(link = "identity"), data = df_train)

# 计算分散参数φ（两种方法二选一即可）
# 方法1：手动计算（无需额外包）
res_pearson <- residuals(my_model, type = "pearson")
phi <- sum(res_pearson^2) / df.residual(my_model)

# 方法2：使用blmeco包的dispersion_glmer
# library(blmeco)
# phi <- dispersion_glmer(my_model)

# 推导全局形状参数
shape <- 1 / phi

# 计算测试集的预测均值μ
pred_mean_test <- predict(my_model, newdata = df_test, type = "response")

# 计算每个测试样本的尺度参数
scale_test <- pred_mean_test / shape

# 计算测试集的对数似然
log_likelihoods_test <- dgamma(df_test$y, shape = shape, scale = scale_test, log = TRUE)

# 查看结果示例
head(log_likelihoods_test)
# 查看平均对数似然
mean(log_likelihoods_test)

重要说明 #

• 形状参数是全局固定值：Gamma GLMER中，shape由模型的分散参数推导，不随单个样本变化；而scale是随样本变化的，由每个样本的预测均值和全局shape计算得出。
• 样本量影响参数精度：原代码中训练集仅100个样本，导致shape估计值与真实值（3）有偏差，增大样本量会显著提升估计精度。

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