如何使用 lme4 或 glmmTMB 拟合准泊松模型？

Question

我正在尝试在 R 中拟合混合效应拟泊松模型。特别是，我正在尝试复制通过 ppml 命令在 stata 中获得的结果。 lme4 不支持准家庭。 link: https://bbolker.github.io/mixedmodels-misc/glmmFAQ.html#fitting-models-with-overdispersion 有拟合拟二项式模型的说明，但没有拟泊松模型。

谢谢！

Answer 1

link 中给出的方法对（准）泊松模型和（准）二项式模型同样有效。关键在于，准似然模型确实代表了对参数的标准误差和相关统计量的post-拟合调整；他们不会（或不应该...）改变模型的拟合方式。

glmer 对“离散响应”（二项式、泊松等）实际上是离散的有点挑剔，但 glmmTMB 是 looser/more 宽容。

这种方法会尽可能多地放置随机效应可以解释的方差，然后对任何剩余的部分进行 post hoc 调整（或下）分散。

我们将使用 grouseticks data set from Elston et al 2001（最初的分析使用观察级随机效应而不是准似然来处理个体观察级的过度分散（=一只小鸡的蜱数，嵌套在 brood 中，嵌套在 location 中。

library(lme4)
g <- transform(grouseticks, sHEIGHT = drop(scale(HEIGHT)))
form <- TICKS~YEAR+sHEIGHT+(1|BROOD)+(1|LOCATION)
full_mod1  <- glmer(form, family="poisson", data=g)

中度过度离散：deviance(full_mod1)/df.residual(full_mod1) 为 1.86。（计算（皮尔逊平方和 residuals/residual df 的比率），正如我们将在下面做的那样，稍微更稳健。

未调整系数table:

printCoefmat(coef(summary(full_mod1)), digits=2)
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)     0.47       0.19     2.4     0.02 *  
YEAR96          1.17       0.23     5.1    4e-07 ***
YEAR97         -0.98       0.25    -3.8    1e-04 ***
sHEIGHT        -0.85       0.12    -6.8    1e-11 ***

现在定义准似然调整函数（如link）：

quasi_table <- function(model,ctab=coef(summary(model))) {
    phi <- sum(residuals(model, type="pearson")^2)/df.residual(model)
    qctab <- within(as.data.frame(ctab),
    {   `Std. Error` <- `Std. Error`*sqrt(phi)
        `z value` <- Estimate/`Std. Error`
        `Pr(>|z|)` <- 2*pnorm(abs(`z value`), lower.tail=FALSE)
    })
    return(qctab)
}

应用它：

printCoefmat(quasi_table(full_mod1),digits=2)
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)     0.47       0.25     1.8    0.065 .  
YEAR96          1.17       0.30     3.8    1e-04 ***
YEAR97         -0.98       0.34    -2.9    0.004 ** 
sHEIGHT        -0.85       0.16    -5.2    3e-07 ***

如前所述，估计值是相同的；标准误差和 p 值已适当放大，z 值已适当缩小。

如果您更喜欢“整洁”的统计信息：

library(tidyverse)
library(broom.mixed)
phi <- sum(residuals(full_mod1, type="pearson")^2)/df.residual(full_mod1)
(full_mod1 
  %>% tidy(effect="fixed")
  %>% transmute(term = term, estimate = estimate,
                std.error = std.error * sqrt(phi),
                statistic = estimate/std.error,
                p.value = 2*pnorm(abs(statistic), lower.tail=FALSE))
)
 term        estimate std.error statistic     p.value
  <chr>          <dbl>     <dbl>     <dbl>       <dbl>
1 (Intercept)    0.467     0.253      1.84 0.0654     
2 YEAR96         1.17      0.303      3.84 0.000121   
3 YEAR97        -0.978     0.336     -2.91 0.00363    
4 sHEIGHT       -0.847     0.164     -5.15 0.000000260

我不知道人们为处理混合模型（例如 merTools、sjstats 等）构建的任何 'downstream' 包是否具有这些功能。可以说 broom.mixed should/could 内置了它。
上面的所有代码应该用 glmmTMB() 代替 glmer 工作得同样好，但我还没有测试它。

如何使用 lme4 或 glmmTMB 拟合准泊松模型？

How do I fit a quasi-poisson model with lme4 or glmmTMB?

r

glm