У меня есть несколько вопросов, касающихся спецификации и интерпретации GLMM. 3 вопроса, безусловно, являются статистическими, а 2 - более конкретно о R. Я пишу здесь, потому что, в конечном счете, я думаю, что проблема заключается в интерпретации результатов GLMM.

В настоящее время я пытаюсь соответствовать GLMM. Я использую данные переписи США из базы данных продольных трактов . Мои наблюдения - это переписные участки. Моя зависимая переменная - это количество вакантных жилищных единиц, и меня интересует связь между вакансией и социально-экономическими переменными. Пример здесь прост, просто используются два фиксированных эффекта: процент небелого населения (раса) и средний доход домохозяйства (класс), а также их взаимодействие. Я хотел бы включить два вложенных случайных эффекта: тракты в течение десятилетий и десятилетий, т. Е. (Декада / тракт). Я рассматриваю их как случайные, чтобы контролировать пространственную (то есть между участками) и временную (то есть между десятилетиями) автокорреляцию. Тем не менее, меня интересует десятилетие как фиксированный эффект, поэтому я включаю его и как фиксированный фактор.

Поскольку моя независимая переменная является неотрицательной переменной целого числа, я пытался подогнать пуассоновские и отрицательные биномиальные GLMM. Я использую журнал общих единиц жилья в качестве зачета. Это означает, что коэффициенты интерпретируются как влияние на уровень вакантных площадей, а не на общее количество вакантных домов.

В настоящее время у меня есть результаты для Пуассона и отрицательного биномиального GLMM, оцененные с использованием glmer и glmer.nb из lme4 . Интерпретация коэффициентов имеет смысл для меня, основываясь на моих знаниях данных и области исследования.

Если вам нужны данные и скрипт, они есть на моем 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

Участки дхарма

Пуассон

Отрицательный бином

Статистические вопросы

Так как я все еще выясняю GLMM, я не уверен в спецификации и интерпретации. У меня есть несколько вопросов:

1. Похоже, мои данные не поддерживают использование модели Пуассона, и поэтому мне лучше с отрицательным биномом. Однако я постоянно получаю предупреждения о том, что мои отрицательные биномиальные модели достигают своего предела итерации, даже когда я увеличиваю максимальный предел. «В theta.ml (Y, mu, weights = object @ resp $ weights, limit = limit,: достигнут предел итерации.» Это происходит с использованием довольно многих различных спецификаций (то есть минимальных и максимальных моделей для фиксированных и случайных эффектов). Я также попытался удалить выбросы в моем иждивенце (брутто, я знаю!), Поскольку верхние 1% значений очень сильно отличаются (нижние 99% находятся в диапазоне от 0 до 1012, верхние 1% - от 1013 до 5213). не влияет на итерации и очень мало влияет на коэффициенты. Я не включаю эти детали здесь. Коэффициенты между Пуассоном и отрицательным биномом также очень похожи. Является ли это отсутствие конвергенции проблемой? Подходит ли отрицательная биноминальная модель? Я также запустил отрицательную биномиальную модель, используяAllFit и не все оптимизаторы выдают это предупреждение (bobyqa, Nelder Mead и nlminbw не сделали).

2. Дисперсия для моего фиксированного эффекта десятилетия постоянно очень низкая или равна 0. Я понимаю, что это может означать, что модель подходит больше. Удаление десятилетия из фиксированных эффектов действительно увеличивает дисперсию случайного эффекта десятилетия до 0,2620 и не оказывает большого влияния на коэффициенты фиксированного эффекта. Что-то не так с этим? Я прекрасно понимаю, что это просто не нужно для объяснения различий между наблюдениями.

3. Эти результаты указывают на то, что я должен попробовать модели с нулевым уровнем инфляции? DHARMa, похоже, предполагает, что нулевая инфляция не может быть проблемой. Если вы думаете, что я все равно должен попробовать, см. Ниже.

R вопросы

4. Я хотел бы попробовать модели с нулевым раздувом, но я не уверен, какой пакет подразумевает вложенные случайные эффекты для нулевых раздуваемых пуассоновских и отрицательных биномиальных GLMM. Я бы использовал glmmADMB для сравнения AIC с моделями с нулевым раздувом, но он ограничен одним случайным эффектом, поэтому не работает для этой модели. Я мог бы попробовать MCMCglmm, но я не знаю байесовской статистики, так что это тоже не привлекательно. Есть еще варианты?

5. Могу ли я отображать экспоненциальные коэффициенты в сводке (модели), или мне нужно делать это вне сводки, как я это сделал здесь?

Я считаю, что есть некоторые важные проблемы, которые необходимо решить с вашей оценкой.

Исходя из того, что я собрал, изучив ваши данные, ваши единицы не сгруппированы географически, т.е. участки переписи внутри округов. Таким образом, использование трактов в качестве группирующего фактора не подходит для захвата пространственной неоднородности, поскольку это означает, что у вас такое же количество людей, что и у групп (или, другими словами, у всех ваших групп есть только одно наблюдение в каждой). Использование многоуровневой стратегии моделирования позволяет нам оценивать дисперсию на индивидуальном уровне, контролируя при этом дисперсию между группами. Поскольку в каждой группе есть только один человек, ваша межгрупповая дисперсия такая же, как и ваша индивидуальная дисперсия, что противоречит цели многоуровневого подхода.

С другой стороны, фактор группировки может представлять повторные измерения во времени. Например, в случае продольного обучения баллы по математике для отдельного человека могут записываться ежегодно, поэтому мы будем иметь годовое значение для каждого учащегося в течение 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.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 Модель ( thetaзначение в объекте glmer.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

Заключительные мысли о нулевой инфляции. Не существует многоуровневой реализации (по крайней мере, мне известно) модели нуассона или Пуассона с нулевым завышением, которая позволяет задавать уравнение для компонента с нулевым завышением. glmmADMBмодель позволяет оценить параметр постоянной нулевой инфляции. Альтернативным подходом было бы использование zeroinflфункции в psclпакете, хотя это не поддерживает многоуровневые модели. Таким образом, вы можете сравнить подгонку отрицательного бинома с одним уровнем к отрицательному биному с завышенным нулевым уровнем. Скорее всего, если нулевая инфляция незначительна для одноуровневых моделей, вполне вероятно, что она не будет существенной для многоуровневой спецификации.

добавление

Если вас беспокоит пространственная автокорреляция, вы можете контролировать это, используя некоторую форму географической взвешенной регрессии (хотя я считаю, что для этого используются точечные данные, а не области). Кроме того, вы можете сгруппировать свои переписные участки по дополнительному фактору группировки (штаты, округа) и включить его в качестве случайного эффекта. Наконец, и я не уверен, что это вполне осуществимо, можно включить пространственную зависимость, используя, например, среднее число R_VACсоседей первого порядка в качестве ковариации. В любом случае, до таких подходов было бы разумно определить, действительно ли присутствует пространственная автокорреляция (с использованием глобальных тестов Морана I, LISA и аналогичных подходов).