Un «modèle d'obstacle» est-il vraiment un modèle? Ou simplement deux modèles séparés et séquentiels?

Considérons un modèle d'obstacles prédisant les données yde comptage à partir d'un prédicteur normal x:

set.seed(1839)
# simulate poisson with many zeros
x <- rnorm(100)
e <- rnorm(100)
y <- rpois(100, exp(-1.5 + x + e))

# how many zeroes?
table(y == 0)

FALSE  TRUE 
   31    69 

Dans ce cas, j'ai des données de comptage avec 69 zéros et 31 comptages positifs. Peu importe pour le moment qu'il s'agit, par définition de la procédure de génération de données, d'un processus de Poisson, car ma question concerne les modèles d'obstacles.

Disons que je veux gérer ces zéros en excès par un modèle d'obstacle. D'après mes lectures à leur sujet, il semblait que les modèles d'obstacles ne sont pas de véritables modèles en soi - ils ne font que deux analyses différentes séquentiellement. Premièrement, une régression logistique prédisant si la valeur est positive ou non par rapport à zéro. Deuxièmement, une régression de Poisson tronquée à zéro n'incluant que les cas non nuls. Cette deuxième étape me semblait mal parce qu'elle (a) jette des données parfaitement bonnes, ce qui (b) pourrait entraîner des problèmes d'alimentation car la plupart des données sont des zéros, et (c) fondamentalement pas un "modèle" en soi , mais exécutant simplement en séquence deux modèles différents.

J'ai donc essayé un «modèle d'obstacle» plutôt que d'exécuter séparément la régression de Poisson logistique et zéro tronquée. Ils m'ont donné des réponses identiques (j'abrège la sortie, par souci de concision):

> # hurdle output
> summary(pscl::hurdle(y ~ x))

Count model coefficients (truncated poisson with log link):
            Estimate Std. Error z value Pr(>|z|)  
(Intercept)  -0.5182     0.3597  -1.441   0.1497  
x             0.7180     0.2834   2.533   0.0113 *

Zero hurdle model coefficients (binomial with logit link):
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  -0.7772     0.2400  -3.238 0.001204 ** 
x             1.1173     0.2945   3.794 0.000148 ***

> # separate models output
> summary(VGAM::vglm(y[y > 0] ~ x[y > 0], family = pospoisson()))

Coefficients: 
            Estimate Std. Error z value Pr(>|z|)  
(Intercept)  -0.5182     0.3597  -1.441   0.1497  
x[y > 0]      0.7180     0.2834   2.533   0.0113 *

> summary(glm(I(y == 0) ~ x, family = binomial))

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)   0.7772     0.2400   3.238 0.001204 ** 
x            -1.1173     0.2945  -3.794 0.000148 ***
---

Cela me semble faux car de nombreuses représentations mathématiques différentes du modèle incluent la probabilité qu'une observation soit non nulle dans l'estimation des cas de comptage positifs, mais les modèles que j'ai exécutés ci-dessus s'ignorent complètement. Par exemple, cela provient du chapitre 5, page 128 des modèles linéaires généralisés de Smithson & Merkle pour les variables dépendantes catégoriques et continues limitées :

... Deuxièmement, la probabilité que assume une valeur (zéro et les entiers positifs) doit être égale à un. Cela n'est pas garanti dans l'équation (5.33). Pour résoudre ce problème, nous multiplions la probabilité de Poisson par la probabilité de succès de Bernoulli .      Ces problèmes nous obligent à exprimer le modèle d'obstacle ci-dessus comme où , ,π P ( Y = y | x , z , β , $y$$\pi$

...λ=exp(

$$P(Y=y|\boldsymbol{x,z,\beta,\gamma}) = \begin{cases} 1-\hat\pi &\text{for } y=0 \\ \hat\pi\times\frac{\exp(-\hat\lambda)\hat\lambda^y/y!}{1-\exp(-\hat\lambda)} &\text{for } y=1,2,\ldots \end{cases} \tag{5.34}$$ π = l o g i t - 1 ( z γ ) x z β $\hat\lambda=\exp(\boldsymbol{x\beta})$$\hat\pi = {\rm logit}^{-1}(\boldsymbol{z\gamma})$$\boldsymbol x$sont les covariables pour le modèle de Poisson, sont les covariables pour le modèle de régression logistique et et sont les coefficients de régression respectifs ... . $\boldsymbol z$$\hat{\boldsymbol{\beta}}$$\hat{\boldsymbol{\gamma}}$

En faisant les deux modèles complètement séparés l'un de l'autre - ce qui semble être ce que font les modèles d'obstacles - je ne vois pas comment est incorporé dans la prédiction des cas de comptage positifs. Mais d'après la façon dont j'ai pu répliquer la fonction en exécutant simplement deux modèles différents, je ne vois pas comment joue un rôle dans le Poisson tronqué régression du tout. logit-1(z γ$\hat{\pi}$hurdle$\text{logit}^{-1}(z\hat{\gamma})$

Suis-je en train de comprendre correctement les modèles d'obstacles? Ils semblent deux exécuter simplement deux modèles séquentiels: premièrement, une logistique; Deuxièmement, un Poisson, ignorant complètement les cas où . J'apprécierais si quelqu'un pouvait éclaircir ma confusion avec l'entreprise .$y = 0$$\hat{\pi}$

Si j'ai raison de dire que c'est ce que sont les modèles d'obstacles, quelle est la définition d'un modèle «d'obstacles», plus généralement? Imaginez deux scénarios différents:

• Imaginez modéliser la compétitivité des courses électorales en examinant les scores de compétitivité (1 - (proportion de votes des gagnants - proportion de votes des finalistes)). C'est [0, 1), car il n'y a pas de liens (par exemple, 1). Un modèle d'obstacle est logique ici, car il y a un processus (a) l'élection était-elle incontestée? et b) si ce n'était pas le cas, quelle était la compétitivité prévue? Nous faisons donc d'abord une régression logistique pour analyser 0 vs (0, 1). Ensuite, nous faisons une régression bêta pour analyser les (0, 1) cas.

• Imaginez une étude psychologique typique. Les réponses sont [1, 7], comme une échelle de Likert traditionnelle, avec un énorme effet de plafond à 7. On pourrait faire un modèle d'obstacle qui est une régression logistique de [1, 7) vs 7, puis une régression Tobit pour tous les cas où les réponses observées sont <7.

Serait-il prudent d'appeler ces deux situations des modèles «d'obstacles» , même si je les estime avec deux modèles séquentiels (logistique puis bêta dans le premier cas, logistique puis Tobit dans le second)?

Séparer la log-vraisemblance

Il est vrai que la plupart des modèles d'obstacles peuvent être estimés séparément (je dirais plutôt que de manière séquentielle ). La raison en est que la log-vraisemblance peut être décomposée en deux parties qui peuvent être maximisées séparément. C'est parce que est juste un facteur d'échelle dans (5.34) qui devient un terme additif dans la log-vraisemblance.$\hat \pi$

Dans la notation de Smithson & Merkle: où est la densité de la distribution de Poisson (non tronquée) et est le facteur de la troncature zéro.

Il devient alors évident que (modèle logit binaire) et (modèle de Poisson tronqué) peuvent être maximisés séparément, conduisant aux mêmes estimations de paramètres, covariances, etc. comme dans le cas où ils sont maximisés conjointement.$\ell_1(\gamma)$$\ell_2(\beta)$

La même logique fonctionne également si la probabilité d'obstacle zéro n'est pas paramétrée via un modèle logit mais tout autre modèle de régression binaire, par exemple, une distribution de comptage censurée à droite à 1. Et, bien sûr, pourrait également être une autre distribution de comptage, par exemple, binôme négatif. La séparation entière ne se décompose que s'il existe des paramètres partagés entre l'obstacle zéro et la partie de comptage tronquée.$\pi$$f(\cdot)$

Un exemple frappant serait l'utilisation de distributions binomiales négatives avec des paramètres séparés mais communs dans les deux composantes du modèle. (Ceci est disponible dans le package de R-Forge, le successeur de l' implémentation.)$\mu$$\theta$hurdle(..., separate = FALSE, dist = "negbin", zero.dist = "negbin")countregpscl

Des questions concrètes

(a) Jeter des données parfaitement bonnes: Dans votre cas, oui, en général non. Vous disposez de données provenant d'un seul modèle de Poisson sans excès de zéros (bien que plusieurs zéros ). Par conséquent, il n'est pas nécessaire d'estimer des modèles distincts pour les zéros et les non-zéros. Cependant, si les deux parties sont réellement motivées par des paramètres différents, il faut en tenir compte.

(b) Peut entraîner des problèmes d'alimentation car la plupart des données sont des zéros: pas nécessairement. Ici, vous avez un tiers des observations qui sont des "succès" (franchissement d'obstacles). Cela ne serait pas considéré comme très extrême dans un modèle de régression binaire. (Bien sûr, s'il n'est pas nécessaire d'estimer des modèles séparés, vous gagneriez en puissance.)

(c) Fondamentalement, pas un «modèle» en soi, mais simplement l'exécution séquentielle de deux modèles différents: c'est plus philosophique et je n'essaierai pas de donner «une» réponse. Au lieu de cela, je soulignerai des points de vue pragmatiques. Pour l' estimation du modèle , il peut être pratique de souligner que les modèles sont séparés car - comme vous le montrez - vous n'avez peut-être pas besoin d'une fonction dédiée pour l'estimation. Pour l' application du modèle , par exemple, pour les prévisions ou les résidus, etc., il peut être plus pratique de voir cela comme un modèle unique.

(d) Serait-il prudent d'appeler ces deux situations des modèles «d'obstacles»: en principe oui. Cependant, le jargon peut varier d'une communauté à l'autre. Par exemple, la régression bêta zéro obstacle est plus communément (et de manière très confuse) appelée régression bêta gonflée zéro. Personnellement, je trouve ce dernier très trompeur car la distribution bêta n'a pas de zéros qui pourraient être gonflés - mais c'est quand même le terme standard dans la littérature. De plus, le modèle tobit est un modèle censuré et donc pas un modèle obstacle. Il pourrait cependant être étendu par un modèle probit (ou logit) plus un modèle normal tronqué . Dans la littérature économétrique, il s'agit du modèle en deux parties de Cragg.

Commentaires sur le logiciel

Le countregpackage sur R-Forge à https://R-Forge.R-project.org/R/?group_id=522 est l'implémentation successive de hurdle()/ zeroinfl()vers pscl. La principale raison pour laquelle il n'est (toujours) pas sur CRAN est que nous voulons réviser l' predict()interface, peut-être d'une manière qui n'est pas entièrement rétrocompatible. Sinon, l'implémentation est assez stable. Par rapport à psclcela, il est livré avec quelques fonctionnalités intéressantes, par exemple:

• Une zerotrunc()fonction qui utilise exactement le même code que hurdle()pour la partie tronquée à zéro du modèle. Ainsi, il offre une alternative à VGAM.

• De plus, il fonctionne en tant que d / p / q / r pour les distributions de comptage tronquées, hachées et gonflées à zéro. Cela facilite de les considérer comme «un» modèle plutôt que comme des modèles séparés.

• Pour évaluer la qualité de l'ajustement, des affichages graphiques comme des rootogrammes et des tracés résiduels quantiles randomisés sont disponibles. (Voir Kleiber & Zeileis, 2016, The American Statistician , 70 (3), 296–303. Doi: 10.1080 / 00031305.2016.1173590 .)

Données simulées

Vos données simulées proviennent d'un seul processus de Poisson. Si eest traité comme un régresseur connu, ce serait un GLM de Poisson standard. S'il es'agit d'une composante de bruit inconnue, il y a une hétérogénéité non observée provoquant un peu de surdispersion qui pourrait être capturée par un modèle binomial négatif ou un autre type de mélange continu ou d'effet aléatoire, etc. Cependant, comme l'effet de eest plutôt faible ici , rien de tout cela ne fait une grande différence. Ci-dessous, je traite ecomme un régresseur (c'est-à-dire avec un vrai coefficient de 1) mais vous pouvez également l'omettre et utiliser des modèles binomiaux ou de Poisson négatifs. Qualitativement, tous ces éléments conduisent à des idées similaires.

## Poisson GLM
p <- glm(y ~ x + e, family = poisson)
## Hurdle Poisson (zero-truncated Poisson + right-censored Poisson)
library("countreg")
hp <- hurdle(y ~ x + e, dist = "poisson", zero.dist = "poisson")
## all coefficients very similar and close to true -1.5, 1, 1
cbind(coef(p), coef(hp, model = "zero"), coef(hp, model = "count"))
##                   [,1]       [,2]      [,3]
## (Intercept) -1.3371364 -1.2691271 -1.741320
## x            0.9118365  0.9791725  1.020992
## e            0.9598940  1.0192031  1.100175

Cela montre que les trois modèles peuvent toujours estimer les vrais paramètres. L'examen des erreurs standard correspondantes montre que dans ce scénario (sans avoir besoin d'une pièce de haie), le GLM de Poisson est plus efficace:

serr <- function(object, ...) sqrt(diag(vcov(object, ...)))
cbind(serr(p), serr(hp, model = "zero"), serr(hp, model = "count"))
##                  [,1]      [,2]      [,3]
## (Intercept) 0.2226027 0.2487211 0.5702826
## x           0.1594961 0.2340700 0.2853921
## e           0.1640422 0.2698122 0.2852902

Les critères d'information standard sélectionneraient le véritable Poisson GLM comme meilleur modèle:

AIC(p, hp)
##    df      AIC
## p   3 141.0473
## hp  6 145.9287

Et un test de Wald permettrait de détecter correctement que les deux composantes du modèle d'obstacle ne sont pas significativement différentes:

hurdletest(hp)
## Wald test for hurdle models
## 
## Restrictions:
## count_((Intercept) - zero_(Intercept) = 0
## count_x - zero_x = 0
## count_e - zero_e = 0
## 
## Model 1: restricted model
## Model 2: y ~ x + e
## 
##   Res.Df Df  Chisq Pr(>Chisq)
## 1     97                     
## 2     94  3 1.0562     0.7877

Enfin, les deux rootogram(p)et qqrplot(p)montrent que le GLM de Poisson correspond très bien aux données et qu'il n'y a pas de zéros en excès ou d'indices sur d'autres erreurs de spécification.

Je suis d'accord que la différence entre les modèles zéro gonflé et les obstacles est difficile à comprendre. Les deux sont une sorte de modèle de mélange. D'après ce que je peux dire, la différence importante est que, dans un modèle gonflé à zéro, vous mélangez une masse à zéro avec une distribution \ textit {qui peut également prendre la valeur zéro}. Pour un modèle d'obstacle, vous mélangez une masse à zéro avec une distribution qui ne prend que des valeurs supérieures à 0. Ainsi, dans le modèle gonflé à zéro, vous pouvez distinguer les «zéros structurels» (correspondant à la masse à zéro) et les «zéros d'échantillonnage» 'correspondant à l'occurrence fortuite de 0 du modèle dans lequel vous mixez. Bien sûr, cette identification dépend fortement du bon choix de distribution! Mais, si vous avez un Poisson gonflé à zéro, par exemple, vous pouvez distinguer les zéros provenant de la composante de Poisson (échantillonnage des zéros) et les zéros provenant de la masse à zéro (zéros structurels). Si vous avez un modèle gonflé à zéro et que la distribution dans laquelle vous mixez n'a pas de masse à zéro, cela pourrait être interprété comme un modèle à obstacle.

En ce qui concerne l'aspect philosophique, "quand devrions-nous considérer quelque chose comme un modèle unique et quand deux modèles distincts" , il peut être intéressant de noter que les estimations d'échantillon des paramètres du modèle sont corrélées.

Dans le graphique ci-dessous avec une simulation, vous voyez principalement la corrélation entre la pente et l'ordonnée à l'origine de la partie des comptes. Mais il existe également une légère relation entre la partie comptes et la partie obstacle. Si vous modifiez les paramètres, par exemple, réduisez le lambda dans la distribution de Poisson ou la taille de l'échantillon, la corrélation devient plus forte.

Je dirais donc que vous ne devriez pas le considérer comme deux modèles distincts . Ou du moins, il existe une certaine relation, même si dans la pratique, vous pouvez calculer les deux estimations indépendamment l'une de l'autre.

set.seed(1839)

Nrep <- 3000
Ns <- 100
pars <- matrix(rep(0,3*Nrep),Nrep)
colnames(pars) <- c("count_intercept","count_slope","hurdle_intercept")

# simulation-loop
# Note that a truncated poisson is used to generate data
# this will make the parameters from the hurdle function easier to interpret and compare
for (i in 1:Nrep) {
  x <- rnorm(Ns,0,1)
  e <- rbinom(Ns,1,exp(-0.7))
  y <- e*truncdist::rtrunc(n=Ns,spec='pois',a=0,b=Inf,lambda=exp(-1.5 + x))
  mod <- pscl::hurdle(y ~ 1+x|1, link="log")
  pars[i,1]<-mod$coefficients$count[1]
  pars[i,2]<-mod$coefficients$count[2]
  pars[i,3]<-mod$coefficients$zero[1]
}  

# viewing data
plotpars <- pars[pars[,1]>-7,] #clipping
pairs(plotpars,cex=0.7,pch=21,
      col= rgb(0,0,0,0.03),
      bg = rgb(0,0,0,0.03))

# demonstrating linear relation / significant correlation
summary(lm(pars[,1] ~ pars[,3]))

Cela n'a pas beaucoup de sens qu'il y ait une corrélation entre les deux parties. Mais cela pourrait probablement être dû à des niveaux discrets d'estimations pour les paramètres du modèle de Poisson et à leur lien avec le nombre de zéro.