使用 lme4 glmer 和 glmer.nb 帮助解释计数数据 GLMM - 负二项式与泊松

机器算法验证 r lme4-nlme 泊松分布负二项分布咕噜咕噜

2022-03-18 13:40:22

我对 GLMM 的规范和解释有一些疑问。3 个问题绝对是统计问题，2 个问题更具体是关于 R。我在这里发帖是因为最终我认为问题是对 GLMM 结果的解释。

我目前正在尝试安装 GLMM。我正在使用来自Longitudinal Tract Database的美国人口普查数据。我的观察是人口普查区。我的因变量是空置住房单元的数量，我对空置率和社会经济变量之间的关系感兴趣。这里的例子很简单，只使用了两个固定效应：非白人人口百分比（种族）和家庭收入中位数（阶级），加上它们的交互作用。我想包括两个嵌套的随机效应：几十年和几十年内的大片，即（十年/大片）。我正在考虑这些随机性以控制空间（即区域之间）和时间（即几十年之间）自相关。但是，我也对作为固定效应的十年感兴趣，所以我也将它作为固定因素包括在内。

由于我的自变量是一个非负整数计数变量，我一直在尝试拟合泊松和负二项式 GLMM。我使用总住房单位的对数作为偏移量。这意味着系数被解释为对空置率的影响，而不是空置房屋的总数。

我目前使用 lme4 中的 glmer 和 glmer.nb 估计了泊松和负二项式 GLMM 的结果。基于我对数据和研究领域的了解，系数的解释对我来说是有意义的。

如果你想要他们在我的Github上的数据和脚本。该脚本包括我在构建模型之前所做的更多描述性调查。

这是我的结果：

泊松模型

Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
 Family: poisson  ( log )
Formula: R_VAC ~ decade + P_NONWHT + a_hinc + P_NONWHT * a_hinc + offset(HU_ln) +      (1 | decade/TRTID10)
   Data: scaled.mydata

     AIC      BIC   logLik deviance df.resid 
 34520.1  34580.6 -17250.1  34500.1     3132 

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-2.24211 -0.10799 -0.00722  0.06898  0.68129 

Random effects:
 Groups         Name        Variance Std.Dev.
 TRTID10:decade (Intercept) 0.4635   0.6808  
 decade         (Intercept) 0.0000   0.0000  
Number of obs: 3142, groups:  TRTID10:decade, 3142; decade, 5

Fixed effects:
                 Estimate Std. Error z value Pr(>|z|)    
(Intercept)     -3.612242   0.028904 -124.98  < 2e-16 ***
decade1980       0.302868   0.040351    7.51  6.1e-14 ***
decade1990       1.088176   0.039931   27.25  < 2e-16 ***
decade2000       1.036382   0.039846   26.01  < 2e-16 ***
decade2010       1.345184   0.039485   34.07  < 2e-16 ***
P_NONWHT         0.175207   0.012982   13.50  < 2e-16 ***
a_hinc          -0.235266   0.013291  -17.70  < 2e-16 ***
P_NONWHT:a_hinc  0.093417   0.009876    9.46  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) dc1980 dc1990 dc2000 dc2010 P_NONWHT a_hinc
decade1980  -0.693                                            
decade1990  -0.727  0.501                                     
decade2000  -0.728  0.502  0.530                              
decade2010  -0.714  0.511  0.517  0.518                       
P_NONWHT     0.016  0.007 -0.016 -0.015  0.006                
a_hinc      -0.023 -0.011  0.023  0.022 -0.009  0.221         
P_NONWHT:_h  0.155  0.035 -0.134 -0.129  0.003  0.155   -0.233
convergence code: 0
Model failed to converge with max|grad| = 0.00181132 (tol = 0.001, component 1)

负二项式模型

Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
 Family: Negative Binomial(25181.5)  ( log )
Formula: R_VAC ~ decade + P_NONWHT + a_hinc + P_NONWHT * a_hinc + offset(HU_ln) +      (1 | decade/TRTID10)
   Data: scaled.mydata

     AIC      BIC   logLik deviance df.resid 
 34522.1  34588.7 -17250.1  34500.1     3131 

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-2.24213 -0.10816 -0.00724  0.06928  0.68145 

Random effects:
 Groups         Name        Variance  Std.Dev. 
 TRTID10:decade (Intercept) 4.635e-01 6.808e-01
 decade         (Intercept) 1.532e-11 3.914e-06
Number of obs: 3142, groups:  TRTID10:decade, 3142; decade, 5

Fixed effects:
                 Estimate Std. Error z value Pr(>|z|)    
(Intercept)     -3.612279   0.028946 -124.79  < 2e-16 ***
decade1980       0.302897   0.040392    7.50 6.43e-14 ***
decade1990       1.088211   0.039963   27.23  < 2e-16 ***
decade2000       1.036437   0.039884   25.99  < 2e-16 ***
decade2010       1.345227   0.039518   34.04  < 2e-16 ***
P_NONWHT         0.175216   0.012985   13.49  < 2e-16 ***
a_hinc          -0.235274   0.013298  -17.69  < 2e-16 ***
P_NONWHT:a_hinc  0.093417   0.009879    9.46  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) dc1980 dc1990 dc2000 dc2010 P_NONWHT a_hinc
decade1980  -0.693                                            
decade1990  -0.728  0.501                                     
decade2000  -0.728  0.502  0.530                              
decade2010  -0.715  0.512  0.517  0.518                       
P_NONWHT     0.016  0.007 -0.016 -0.015  0.006                
a_hinc      -0.023 -0.011  0.023  0.022 -0.009  0.221         
P_NONWHT:_h  0.154  0.035 -0.134 -0.129  0.003  0.155   -0.233

泊松 DHARMa 检验

    One-sample Kolmogorov-Smirnov test

data:  simulationOutput$scaledResiduals
D = 0.044451, p-value = 8.104e-06
alternative hypothesis: two-sided

    DHARMa zero-inflation test via comparison to expected zeros with simulation under H0 = fitted model

data:  simulationOutput
ratioObsExp = 1.3666, p-value = 0.159
alternative hypothesis: more

负二项式 DHARMa 检验

    One-sample Kolmogorov-Smirnov test

data:  simulationOutput$scaledResiduals
D = 0.04263, p-value = 2.195e-05
alternative hypothesis: two-sided

    DHARMa zero-inflation test via comparison to expected zeros with simulation under H0 = fitted model

data:  simulationOutput2
ratioObsExp = 1.376, p-value = 0.174
alternative hypothesis: more

DHARMa 地块

泊松

负二项式

统计问题

由于我仍在研究 GLMM，我对规范和解释感到不安全。我有一些疑问：

看来我的数据不支持使用泊松模型，因此我最好使用负二项式。但是，我一直收到警告说我的负二项式模型达到了它们的迭代限制，即使我增加了最大限制。“在 theta.ml(Y, mu, weights = object@resp$weights, limit = limit, : 达到迭代限制。” 使用相当多的不同规范（即固定和随机效应的最小和最大模型）会发生这种情况。我也尝试过删除我的依赖项中的异常值（我知道，总额！），因为前 1% 的值非常异常（底部 99% 的范围从 0-1012，顶部的 1% 从 1013-5213）。那没有对迭代没有任何影响，对系数也几乎没有影响。我在这里不包括这些细节。泊松和负二项式之间的系数也非常相似。这种缺乏收敛性是一个问题吗？负二项式模型是否合适？我还使用了运行负二项式模型AllFit并不是所有的优化器都会抛出这个警告（bobyqa、Nelder Mead 和 nlminbw 没有）。
我的十年固定效应的方差始终非常低或为 0。我理解这可能意味着模型过度拟合。从固定效应中取出十年确实会将十年随机效应方差增加到 0.2620，并且对固定效应系数没有太大影响。把它留在里面有什么问题吗？我很好地将其解释为根本不需要在观察方差之间进行解释。
这些结果是否表明我应该尝试零膨胀模型？DHARMa 似乎暗示零通胀可能不是问题。如果您认为我还是应该尝试，请参见下文。

问题

我愿意尝试零膨胀模型，但我不确定哪个包对零膨胀泊松和负二项式 GLMM 实施嵌套随机效应。我会使用 glmmADMB 将 AIC 与零膨胀模型进行比较，但它仅限于单个随机效应，因此不适用于该模型。我可以尝试 MCMCglmm，但我不知道贝叶斯统计，所以这也没有吸引力。还有其他选择吗？
我可以在摘要（模型）中显示指数系数，还是必须像我在这里所做的那样在摘要之外进行？

1个回答

我相信你的估计有一些重要的问题需要解决。

根据我通过检查您的数据收集的信息，您的单位没有按地理分组，即县内的人口普查区。因此，使用区域作为分组因子不适合捕捉空间异质性，因为这意味着您拥有与组相同数量的个体（或者换句话说，您的所有组每个只有一个观察值）。使用多级建模策略使我们能够估计个体级方差，同时控制组间方差。由于您的组每个只有一个人，因此您的组间差异与您的个人级别差异相同，从而违背了多级别方法的目的。

另一方面，分组因子可以表示随时间的重复测量。例如，在纵向研究的情况下，个人的“数学”分数可能每年记录一次，因此我们将为每个学生设置 n 年的年度值（在这种情况下，分组因素是学生，因为我们有 n “嵌套”在学生中的观察数量）。在您的情况下，您对每个人口普查区重复测量了decade. 因此，您可以使用您的TRTID10变量作为分组因子来捕捉“十年间的差异”。这导致 3142 个观测值嵌套在 635 个区域中，每个人口普查区域大约有 4 和 5 个观测值。

正如评论中提到的，decade用作分组因素不是很合适，因为每个人口普查区只有大约 5 个十年，并且可以更好地捕捉它们的影响decade作为协变量引入。

其次，确定您的数据是否应该使用泊松或负二项式模型（或零膨胀方法）进行建模。考虑数据中过度分散的数量。泊松分布的基本特征是等分散性，这意味着均值等于分布的方差。查看您的数据，很明显存在过度分散。方差远大于均值。

library(dplyr)    
 dispersionstats <- scaled.mydata %>%
 + group_by(decade) %>%
 + summarise(
 + means = mean(R_VAC),
 + variances = var(R_VAC),
 + ratio = variances/means)

##   dispersionstats
##   # A tibble: 5 x 5
##   decade     means variances     ratio 
##    <int>     <dbl>     <dbl>     <dbl> 
## 1   1970  45.43513   4110.89  90.47822 
## 2   1980 103.52365  17323.34 167.33707 
## 3   1990 177.68038  62129.65 349.67087 
## 4   2000 190.23150  91059.60 478.67784 
## 5   2010 247.68246 126265.60 509.78821

尽管如此，为了确定负二项式在统计上是否更合适，标准方法是在泊松模型和负二项式模型之间进行似然比检验，这表明负二项式模型更适合。

library(MASS)
library(lmtest)

modelformula <- formula(R_VAC ~ factor(decade) + P_NONWHT * a_hinc + offset(HU_ln))

poismodel <- glm(modelformula, data = scaled.mydata, family = "poisson")   
nbmodel <- glm.nb(modelformula, data = scaled.mydata)

lrtest(poismodel, nbmodel)

## Likelihood ratio test

##  Model 1: R_VAC ~ factor(decade) + P_NONWHT * a_hinc + offset(HU_ln)  
## Model 2: R_VAC ~ factor(decade) + P_NONWHT * a_hinc + offset(HU_ln)
##   #Df  LogLik Df  Chisq Pr(>Chisq)
## 1   8 -154269
## 2   9  -17452  1 273634  < 2.2e-16 ***
##  ---
## Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

在确定这一点之后，下一个测试可以考虑使用类似方法是否需要多级（混合模型）方法，这表明多级版本提供了更好的拟合。（可以使用类似的测试来比较 glmer 拟合假设泊松分布到 glmer.nb 对象，只要模型在其他方面相同。）

library(lme4)

glmmformula <- update(modelformula, . ~ . + (1|TRTID10))

nbglmm <- glmer.nb(glmmformula, data = scaled.mydata)

lrtest(nbmodel, nbglmm)

## Model 1: R_VAC ~ factor(decade) + P_NONWHT * a_hinc + offset(HU_ln)
## Model 2: R_VAC ~ factor(decade) + P_NONWHT + a_hinc + (1 | TRTID10) +
##     P_NONWHT:a_hinc + offset(HU_ln)
##   #Df LogLik Df Chisq Pr(>Chisq)
## 1   9 -17452
## 2  10 -17332  1 239.3  < 2.2e-16 ***
## ---
## Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

关于泊松模型和 nb 模型的估计，它们实际上预计彼此非常相似，主要区别是标准误差，即如果存在过度离散，泊松模型往往会提供有偏差的标准误差。以您的数据为例：

poissonglmm <- glmer(glmmformula, data = scaled.mydata)
summary(poissonglmm)

## Random effects:
##  Groups  Name        Variance Std.Dev.
## TRTID10 (Intercept) 0.2001   0.4473
## Number of obs: 3142, groups:  TRTID10, 635

## Fixed effects:
##                     Estimate Std. Error z value Pr(>|z|)
## (Intercept)        -2.876013   0.020602 -139.60   <2e-16 ***
## factor(decade)1980  0.092597   0.007602   12.18   <2e-16 ***
## factor(decade)1990  0.903543   0.007045  128.26   <2e-16 ***
## factor(decade)2000  0.854821   0.006913  123.65   <2e-16 ***
## factor(decade)2010  0.986126   0.006723  146.67   <2e-16 ***
## P_NONWHT           -0.125500   0.014007   -8.96   <2e-16 ***
## a_hinc             -0.107335   0.001480  -72.52   <2e-16 ***
## P_NONWHT:a_hinc     0.160937   0.003117   51.64   <2e-16 ***
## ---
## Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

summary(nbglmm)
## Random effects:
##  Groups  Name        Variance Std.Dev.
##  TRTID10 (Intercept) 0.09073  0.3012
## Number of obs: 3142, groups:  TRTID10, 635

## Fixed effects:
##                     Estimate Std. Error z value Pr(>|z|)
## (Intercept)        -2.797861   0.056214  -49.77  < 2e-16 ***
## factor(decade)1980  0.118588   0.039589    3.00  0.00274 **
## factor(decade)1990  0.903440   0.038255   23.62  < 2e-16 ***
## factor(decade)2000  0.843949   0.038172   22.11  < 2e-16 ***
## factor(decade)2010  1.068025   0.037376   28.58  < 2e-16 ***
## P_NONWHT            0.020012   0.089224    0.22  0.82253
## a_hinc             -0.129094   0.008109  -15.92  < 2e-16 ***
## P_NONWHT:a_hinc     0.149223   0.018967    7.87 3.61e-15 ***
## ---
## Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

请注意系数估计值非常相似，主要区别仅在于您的协变量之一的显着性，以及随机效应方差的差异，这表明 nb 中的过度分散参数捕获的单位级方差模型（thetaglmer.nb 对象中的值）捕获随机效应捕获的一些区域间方差。

关于指数系数（和相关的置信区间），您可以使用以下内容：

fixed <- fixef(nbglmm)
confnitfixed <- confint(nbglmm, parm = "beta_", method = "Wald") # Beware: The Wald method is less accurate but much, much faster.

# The exponentiated coefficients are also known as Incidence Rate Ratios (IRR)
IRR <- exp(cbind(fixed, confintfixed)
IRR
##                         fixed      2.5 %     97.5 %
## (Intercept)        0.06094028 0.05458271 0.06803835
## factor(decade)1980 1.12590641 1.04184825 1.21674652
## factor(decade)1990 2.46807856 2.28979339 2.66024515
## factor(decade)2000 2.32553168 2.15789585 2.50619029
## factor(decade)2010 2.90962703 2.70410073 3.13077444
## P_NONWHT           1.02021383 0.85653208 1.21517487
## a_hinc             0.87889172 0.86503341 0.89297205
## P_NONWHT:a_hinc    1.16093170 1.11856742 1.20490048

最后的想法，关于零通胀。没有多级实现（至少我知道）零膨胀泊松或 negbin 模型允许您为混合物的零膨胀分量指定方程。该glmmADMB模型可让您估计一个恒定的零通胀参数。另一种方法是使用包zeroinfl中的函数pscl，尽管这不支持多级模型。因此，您可以将单级负二项式的拟合与单级零膨胀负二项式进行比较。如果零通胀对单级模型不显着，则很可能对多级规范不显着。

附录

如果您担心空间自相关，您可以使用某种形式的地理加权回归来控制这一点（尽管我相信这使用点数据，而不是区域）。或者，您可以按附加分组因素（州、县）对人口普查区域进行分组，并将其作为随机效应包含在内。最后，我不确定这是否完全可行，可以使用例如R_VAC一阶邻居的平均计数作为协变量来合并空间依赖性。无论如何，在这种方法之前，确定空间自相关是否确实存在（使用 Global Moran's I、LISA 测试和类似方法）是明智的。

其它你可能感兴趣的问题

上一篇随机与固定效应：为什么斜率的标准误差相同，但截距却大不相同？下一篇矩生成函数的实际用途