J'ai des questions concernant la spécification et l'interprétation des GLMM. 3 questions sont définitivement statistiques et 2 sont plus spécifiquement sur R. Je poste ici parce que finalement je pense que le problème est l'interprétation des résultats GLMM.
J'essaie actuellement d'installer un GLMM. J'utilise les données du recensement américain de la Longitudinal Tract Database . Mes observations sont des secteurs de recensement. Ma variable dépendante est le nombre de logements vacants et je m'intéresse à la relation entre la vacance et les variables socio-économiques. L'exemple ici est simple, utilisant simplement deux effets fixes: le pourcentage de la population non blanche (race) et le revenu médian des ménages (classe), plus leur interaction. Je voudrais inclure deux effets aléatoires imbriqués: les tracts dans les décennies et les décennies, c'est-à-dire (décennie / tract). Je les considère aléatoires dans un effort pour contrôler l'autocorrélation spatiale (c'est-à-dire entre les voies) et temporelle (c'est-à-dire entre les décennies). Cependant, je suis intéressé par la décennie comme effet fixe aussi, donc je l'inclus aussi comme facteur fixe.
Étant donné que ma variable indépendante est une variable de nombre entier non négatif, j'ai essayé d'adapter le poisson et les GLMM binomiaux négatifs. J'utilise le journal du nombre total de logements comme compensation. Cela signifie que les coefficients sont interprétés comme l'effet sur le taux d'inoccupation et non sur le nombre total de logements vacants.
J'ai actuellement des résultats pour un Poisson et un GLMM binomial négatif estimés en utilisant glmer et glmer.nb de lme4 . L'interprétation des coefficients a du sens pour moi en fonction de ma connaissance des données et de la zone d'étude.
Si vous voulez que les données et le script soient sur mon Github . Le script comprend plus des investigations descriptives que j'ai faites avant de construire les modèles.
Voici mes résultats:
Modèle de Poisson
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)
Modèle binomial négatif
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
Tests Poisson 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
Tests DHARMa binomiaux négatifs
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 parcelles
Poisson
Binôme négatif
Questions statistiques
Comme je suis encore en train de trouver des GLMM, je ne me sens pas sûr de la spécification et de l'interprétation. J'ai quelques questions:
Il semble que mes données ne prennent pas en charge l'utilisation d'un modèle de Poisson et donc je suis mieux avec un binôme négatif. Cependant, je reçois constamment des avertissements que mes modèles binomiaux négatifs atteignent leur limite d'itération, même lorsque j'augmente la limite maximale. "Dans theta.ml (Y, mu, poids = objet @ resp $ poids, limite = limite,: limite d'itération atteinte." Cela se produit en utilisant plusieurs spécifications différentes (c'est-à-dire des modèles minimaux et maximaux pour les effets fixes et aléatoires). J'ai également essayé de supprimer les valeurs aberrantes dans ma personne à charge (brut, je sais!), Car les 1% supérieurs de valeurs sont très éloignés (plage 99% inférieure de 0-1012, 1% supérieure de 1013-5213). t avoir aucun effet sur les itérations et très peu d'effet sur les coefficients. Je n'inclus pas ces détails ici. Les coefficients entre Poisson et binôme négatif sont également assez similaires. Ce manque de convergence est-il un problème? Le modèle binomial négatif est-il un bon ajustement? J'ai également exécuté le modèle binomial négatif en utilisantAllFit et tous les optimiseurs ne lancent pas cet avertissement (bobyqa, Nelder Mead et nlminbw non).
La variance pour mon effet fixe de la décennie est toujours très faible ou nulle. Je comprends que cela pourrait signifier que le modèle est surajusté. La suppression de la décennie des effets fixes augmente la variance de l'effet aléatoire de la décennie à 0,2620 et n'a pas beaucoup d'effet sur les coefficients d'effet fixe. Y a-t-il quelque chose de mal à le laisser? Je suis bien en train de l'interpréter comme n'étant tout simplement pas nécessaire pour expliquer la variance d'observation.
Ces résultats indiquent-ils que je devrais essayer des modèles à gonflage nul? DHARMa semble suggérer que l'inflation zéro pourrait ne pas être le problème. Si vous pensez que je devrais essayer de toute façon, voir ci-dessous.
Questions R
Je serais prêt à essayer des modèles à gonflement nul, mais je ne sais pas quel package implémente des effets aléatoires imbriqués pour les GLMM de Poisson et les GLMM binomiaux négatifs. J'utiliserais glmmADMB pour comparer l'AIC avec des modèles gonflés à zéro, mais il est limité à un seul effet aléatoire, donc ne fonctionne pas pour ce modèle. Je pourrais essayer MCMCglmm, mais je ne connais pas les statistiques bayésiennes donc ce n'est pas non plus attrayant. D'autres options?
Puis-je afficher des coefficients exponentiels dans le résumé (modèle), ou dois-je le faire en dehors du résumé comme je l'ai fait ici?
la source
decade
à la fois fixe et aléatoire n'a pas de sens. Soit le fixer comme fixe et l'inclure uniquement(1 | decade:TRTID10)
comme aléatoire (ce qui équivaut à(1 | TRTID10)
supposer que votre niveauTRTID10
n'a pas les mêmes pendant des décennies différentes), soit le supprimer des effets fixes. Avec seulement 4 niveaux, vous feriez mieux de le faire réparer: la recommandation habituelle est d'ajuster des effets aléatoires si l'on a 5 niveaux ou plus.bobyqa
optimiseur et qu'il n'a produit aucun avertissement. Quel est le problème alors? Utilisez simplementbobyqa
.bobyqa
converge mieux que l'optimiseur par défaut (et je pense avoir lu quelque part qu'il va devenir par défaut dans les futures versions delme4
). Je ne pense pas que vous ayez à vous soucier de la non-convergence avec l'optimiseur par défaut s'il converge avecbobyqa
.Réponses:
Je pense qu'il y a des problèmes importants à résoudre avec votre estimation.
D'après ce que j'ai recueilli en examinant vos données, vos unités ne sont pas regroupées géographiquement, c'est-à-dire les secteurs de recensement au sein des comtés. Ainsi, l'utilisation de tracts comme facteur de regroupement n'est pas appropriée pour capturer l'hétérogénéité spatiale car cela signifie que vous avez le même nombre d'individus que de groupes (ou, pour le dire autrement, tous vos groupes n'ont qu'une seule observation chacun). L'utilisation d'une stratégie de modélisation à plusieurs niveaux nous permet d'estimer la variance au niveau individuel, tout en contrôlant la variance entre les groupes. Étant donné que vos groupes n'ont qu'un seul individu chacun, votre variance entre les groupes est la même que votre variance au niveau individuel, ce qui va à l'encontre de l'objectif de l'approche à plusieurs niveaux.
En revanche, le facteur de regroupement peut représenter des mesures répétées dans le temps. Par exemple, dans le cas d'une étude longitudinale, les scores "mathématiques" d'un individu peuvent être enregistrés annuellement, donc nous aurions une valeur annuelle pour chaque élève pendant n années (dans ce cas, le facteur de regroupement est l'élève comme nous avons n nombre d’observations «imbriquées» chez les élèves). Dans votre cas, vous avez répété des mesures de chaque secteur de recensement par
decade
. Ainsi, vous pouvez utiliser votreTRTID10
variable comme facteur de regroupement pour saisir "la variance entre les décennies". Cela conduit à 3142 observations nichées dans 635 secteurs, avec environ 4 et 5 observations par secteur de recensement.Comme mentionné dans un commentaire, l'utilisation en
decade
tant que facteur de regroupement n'est pas très appropriée, car vous ne disposez que d'environ 5 décennies pour chaque secteur de recensement, et leur effet peut être mieux saisi en introduisantdecade
comme covariable.Deuxièmement, pour déterminer si vos données doivent être modélisées en utilisant un modèle binomial poisson ou négatif (ou une approche gonflée zéro). Considérez la quantité de surdispersion dans vos données. La caractéristique fondamentale d'une distribution de Poisson est l'équidispersion, ce qui signifie que la moyenne est égale à la variance de la distribution. En regardant vos données, il est assez clair qu'il y a beaucoup de surdispersion. Les écarts sont beaucoup plus importants que les moyennes.
Néanmoins, pour déterminer si le binôme négatif est plus approprié statistiquement, une méthode standard consiste à effectuer un test de rapport de vraisemblance entre un modèle de Poisson et un modèle binomial négatif, ce qui suggère que le négbin est mieux adapté.
Après avoir établi cela, un prochain test pourrait déterminer si l'approche à plusieurs niveaux (modèle mixte) est justifiée en utilisant une approche similaire, ce qui suggère que la version à plusieurs niveaux offre un meilleur ajustement. (Un test similaire pourrait être utilisé pour comparer un ajustement glmer en supposant une distribution de poisson à l'objet glmer.nb, tant que les modèles sont par ailleurs les mêmes.)
En ce qui concerne les estimations des modèles poisson et nb, elles devraient en fait être très similaires, la principale distinction étant les erreurs standard, c'est-à-dire qu'en cas de surdispersion, le modèle poisson tend à fournir des erreurs standard biaisées. Prenons vos données comme exemple:
Remarquez comment les estimations des coefficients sont toutes très similaires, la principale différence n'étant que la signification de l'une de vos covariables, ainsi que la différence dans la variance des effets aléatoires, ce qui suggère que la variance au niveau de l'unité saisie par le paramètre de surdispersion dans le nb Le modèle (la
theta
valeur dans l'objet glmer.nb) capture une partie de la variance entre les voies capturée par les effets aléatoires.En ce qui concerne les coefficients exponentiels (et les intervalles de confiance associés), vous pouvez utiliser les éléments suivants:
Réflexions finales concernant une inflation zéro. Il n'y a pas d'implémentation à plusieurs niveaux (du moins que je sache) d'un modèle poisson ou negbin gonflé zéro qui vous permet de spécifier une équation pour le composant gonflé zéro du mélange. le
glmmADMB
modèle vous permet d'estimer un paramètre d'inflation zéro constant. Une autre approche serait d'utiliser lazeroinfl
fonction dans lepscl
package, bien que cela ne prenne pas en charge les modèles à plusieurs niveaux. Ainsi, vous pouvez comparer l'ajustement d'un binôme négatif à un seul niveau au binôme négatif gonflé à zéro à un seul niveau. Il est probable que si l'inflation nulle n'est pas significative pour les modèles à un seul niveau, il est probable qu'elle ne serait pas significative pour la spécification à plusieurs niveaux.Addenda
Si vous êtes préoccupé par l'autocorrélation spatiale, vous pouvez contrôler cela en utilisant une certaine forme de régression géographique pondérée (bien que je pense que cela utilise des données ponctuelles, pas des zones). Vous pouvez également regrouper vos secteurs de recensement en fonction d'un facteur de regroupement supplémentaire (États, comtés) et l'inclure comme effet aléatoire. Enfin, et je ne sais pas si cela est tout à fait possible, il peut être possible d'incorporer la dépendance spatiale en utilisant, par exemple, le nombre moyen de
R_VAC
voisins du premier ordre comme covariable. Dans tous les cas, avant de telles approches, il serait judicieux de déterminer si une autocorrélation spatiale est effectivement présente (en utilisant les tests Global Moran I, LISA et des approches similaires).la source
brms
peut s'adapter à des modèles binomiaux négatifs gonflés à zéro avec des effets aléatoires.brms
et de le comparer au modèle glmer.nb comme indiqué ci-dessus. J'essaierai également d'inclure un lieu défini par le recensement (essentiellement une municipalité, ~ 170 groupes) comme facteur de regroupement pour les effets aléatoires (seulement 5 comtés dans les données, donc je n'utiliserai pas cela). Je vais également tester l'autocorrélation spatiale des résidus en utilisant le Global Moran I. Je ferai rapport quand je l'aurai fait.brms
et en les comparant.