Um “modelo de barreira” é realmente um modelo? Ou apenas dois modelos sequenciais separados?

Considere um modelo de obstáculos que prevê dados de contagem yde um preditor 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

Nesse caso, tenho dados de contagem com 69 zeros e 31 contagens positivas. Não importa, por enquanto, que este é, por definição do procedimento de geração de dados, um processo de Poisson, porque minha pergunta é sobre modelos de obstáculos.

Digamos que eu queira lidar com esses zeros em excesso por um modelo de barreira. Pela minha leitura sobre eles, parecia que os modelos de obstáculos não são modelos reais em si - eles estão apenas fazendo duas análises diferentes sequencialmente. Primeiro, uma regressão logística prevendo se o valor é positivo versus zero. Segundo, uma regressão de Poisson truncada com zero, incluindo apenas os casos diferentes de zero. Este segundo passo me pareceu errado porque é (a) jogar fora dados perfeitamente bons, o que (b) pode levar a problemas de energia, uma vez que muitos dados são zeros e (c) basicamente não é um "modelo" por si só , mas apenas executando sequencialmente dois modelos diferentes.

Então, tentei um "modelo de barreira" em vez de apenas executar a regressão de Poisson logística e truncada a zero separadamente. Eles me deram respostas idênticas (estou abreviando a saída, por uma questão de brevidade):

> # 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 ***
---

Isso me parece uma vez que muitas representações matemáticas diferentes do modelo incluem a probabilidade de que uma observação seja diferente de zero na estimativa de casos de contagem positiva, mas os modelos que eu executei acima se ignoram completamente. Por exemplo, isso é do Capítulo 5, página 128 dos Modelos Lineares Generalizados da Smithson & Merkle para Variáveis Dependentes Limitadas Categóricas e Contínuas :

... Segundo, a probabilidade de assumir qualquer valor (zero e números inteiros positivos) deve ser igual a um. Isso não é garantido na Equação (5.33). Para lidar com esta questão, multiplicamos a probabilidade de Poisson pela probabilidade de sucesso de Bernoulli . Esses problemas exigem que expressemos o modelo de barreira acima como que , , $y$ $\pi$

$\begin{matrix} (5.34) & P (Y = y | x, z, β, γ) = {\begin{cases} 1 - \hat{π} & for y = 0 \\ \hat{π} \times \frac{\exp (- \hat{λ}) {\hat{λ}}^{y} / y!}{1 - \exp (- \hat{λ})} & for y = 1, 2, \dots \end{cases} \end{matrix}$ $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}$ $\hat\lambda=\exp(\boldsymbol{x\beta})$ $\hat\pi = {\rm logit}^{-1}(\boldsymbol{z\gamma})$ $\boldsymbol x$ são as covariáveis do modelo de Poisson, são as covariáveis do modelo de regressão logística e e são os respectivos coeficientes de regressão ... . $\boldsymbol z$ $\hat{\boldsymbol{\beta}}$ $\hat{\boldsymbol{\gamma}}$

Fazendo os dois modelos completamente separados um do outro - o que parece ser o que os modelos de obstáculos fazem - não vejo como é incorporado na previsão de casos de contagem positiva. Mas com base em como consegui replicar a função executando apenas dois modelos diferentes, não vejo como desempenha um papel no Poisson truncado regressão. $\hat{\pi}$ hurdle $\text{logit}^{-1}(z\hat{\gamma})$

Estou entendendo os modelos de obstáculos corretamente? Eles parecem estar apenas executando dois modelos seqüenciais: primeiro, uma logística; Segundo, um Poisson, ignorando completamente os casos em que . Eu apreciaria se alguém pudesse esclarecer minha confusão com o negócio . $y = 0$ $\hat{\pi}$

Se estou certo de que é isso que são os modelos de obstáculos, qual é a definição de um modelo de "obstáculos", de maneira mais geral? Imagine dois cenários diferentes:

Imagine modelar a competitividade das raças eleitorais observando as pontuações de competitividade (1 - (proporção de votos do vencedor - proporção de votos do segundo colocado)). Isso é [0, 1), porque não há laços (por exemplo, 1). Um modelo de obstáculos faz sentido aqui, porque há um processo (a) a eleição não foi contestada? e (b) se não fosse, o que previa competitividade? Então, primeiro fazemos uma regressão logística para analisar 0 vs. (0, 1). Em seguida, fazemos regressão beta para analisar os casos (0, 1).
Imagine um estudo psicológico típico. As respostas são [1, 7], como uma escala Likert tradicional, com um enorme efeito teto em 7. Pode-se fazer um modelo de barreira com regressão logística de [1, 7) vs. 7 e, em seguida, uma regressão Tobit para todos os casos em que as respostas observadas são <7.

Seria seguro chamar essas duas situações de modelos de "obstáculo" , mesmo se eu os estimar com dois modelos seqüenciais (logística e beta no primeiro caso, logística e Tobit no segundo)?

r count-data zero-inflation Mark White
fonte

Acredito que os modelos de obstáculos são equivalentes à execução dos dois modelos separados (binários + zero-truncados). A razão técnica para isso funcionar é que o primeiro modelo usa apenas zero / não zero para estimar ; o segundo modelo condiciona em uma resposta diferente de zero para estimar .

π

$\pi$

λ

$\lambda$

precisa

Então seria então para cada cujo ?

\hat{π}

$\hat{\pi}$

1

$1$

i

$i$

y > 0

$y > 0$

Mark White

Não. O modelo condicional exclui o termo , ou seja, ...

\hat{π}

$\hat \pi$

P (Y = y | Y > 0) = \exp (- \hat{λ}) etc.

$P(Y=y|Y>0) = \exp(-\hat \lambda) \textrm{etc.}$

Ben Bolker

Ah obrigado. Portanto, suponho que a equação de Smithson e Merkle descreva um modelo diferente do que é implementado pscl::hurdle, mas parece o mesmo na Equação 5 aqui: cran.r-project.org/web/packages/pscl/vignettes/countreg.pdf Ou talvez eu ainda falta algo básico que faça clique para mim?

Mark1:

É o mesmo modelo. Mike e Ed focam no caso mais simples (logit + Poisson), que é o padrão hurdle(). Em nossa pareada / vinheta, tentamos enfatizar os blocos de construção mais gerais.

Achim Zeileis

Respostas:

Separando a probabilidade de log

É correto que a maioria dos modelos de obstáculos possa ser estimada separadamente (eu diria, em vez de sequencialmente ). O motivo é que a probabilidade do log pode ser decomposta em duas partes que podem ser maximizadas separadamente. Isso ocorre porque é apenas um fator de escala em (5.34) que se torna um termo aditivo na probabilidade de log. $\hat \pi$

Na notação de Smithson & Merkle: que é a densidade da distribuição Poisson (não truncada) e é o fator do truncamento zero.

\begin{array}{rcl} ℓ (β, γ; y, x, z) & = & ℓ_{1} (γ; y, z) + ℓ_{2} (β; y, x) \\ = & \sum_{i : y_{i} = 0} \log {1 - {l o g i t}^{- 1} (z_{i}^{⊤} γ)} + \sum_{i : y_{i} > 0} \log {{l o g i t}^{- 1} (z_{i}^{⊤} γ)} + \\ \sum_{i : y_{i} > 0} [\log {f (y_{i}; \exp (x_{i}^{⊤} β)} - \log {1 - f (0; \exp (x_{i}^{⊤} β)}] \end{array}

$\begin{eqnarray*} \ell(\beta, \gamma; y, x, z) & = & \ell_1(\gamma; y, z) + \ell_2(\beta; y, x) \\ & = & \sum_{i: y_i = 0} \log\left\{1 - \mathrm{logit}^{-1}(z_i^\top \gamma)\right\} + \sum_{i: y_i > 0} \log\left\{\mathrm{logit}^{-1}(z_i^\top \gamma)\right\} + \\ & & \sum_{i: y_i > 0} \left[ \log \left\{f(y_i; \exp(x_i^\top \beta)\right\} - \log\left\{ 1 - f(0; \exp(x_i^\top \beta)\right\}\right] \end{eqnarray*}$

f (y; λ) = \exp (- λ) λ^{y} / y!

$f(y; \lambda) = \exp(-\lambda) \lambda^y/y!$

1 - f (0; λ) = 1 - \exp (- λ)

$1 - f(0; \lambda) = 1 - \exp(-\lambda)$

Torna-se óbvio que (modelo de logit binário) e (modelo de Poisson truncado com zero) podem ser maximizados separadamente, levando às mesmas estimativas de parâmetros, covariâncias etc. como no caso onde eles são maximizados em conjunto. $\ell_1(\gamma)$ $\ell_2(\beta)$

A mesma lógica também funciona se a probabilidade de obstáculo zero não for parametrizada por meio de um modelo de logit, mas qualquer outro modelo de regressão binária, por exemplo, uma distribuição de contagem censurada à direita em 1. E, é claro, também pode ser outra distribuição de contagem, por exemplo, binomial negativo. Toda a separação é interrompida apenas se houver parâmetros compartilhados entre o obstáculo zero e a parte da contagem truncada. $\pi$ $f(\cdot)$

Um exemplo proeminente seria se distribuições binomiais negativas com parâmetros mas comuns fossem empregados nos dois componentes do modelo. (Isto está disponível em no pacote de R-Forge, o sucessor da implementação.) $\mu$ $\theta$ hurdle(..., separate = FALSE, dist = "negbin", zero.dist = "negbin")countregpscl

Questões concretas

(a) Jogando fora dados perfeitamente bons: no seu caso, sim, em geral não. Você tem dados de um único modelo de Poisson sem excesso de zeros (embora muitos zeros ). Portanto, não é necessário estimar modelos separados para zeros e não zeros. No entanto, se as duas partes são realmente orientadas por parâmetros diferentes, é necessário levar isso em consideração.

(b) Pode levar a problemas de energia, já que muitos dados são zeros: não necessariamente. Aqui, você tem um terço das observações que são "sucessos" (cruzamentos de obstáculos). Isso não seria considerado muito extremo em um modelo de regressão binária. (Obviamente, se não for necessário estimar modelos separados, você ganhará energia.)

(c) Basicamente, não é um "modelo" por si só, mas apenas executa sequencialmente dois modelos diferentes: isso é mais filosófico e não tentarei dar "uma" resposta. Em vez disso, vou apontar pontos de vista pragmáticos. Para estimativa de modelo , pode ser conveniente enfatizar que os modelos são separados porque - como você mostra - talvez você não precise de uma função dedicada para a estimativa. Para a aplicação do modelo , por exemplo, para previsões ou resíduos, etc., pode ser mais conveniente ver isso como um modelo único.

(d) Seria seguro chamar essas duas situações de modelos de 'obstáculo': em princípio, sim. No entanto, o jargão pode variar entre as comunidades. Por exemplo, a regressão beta com zero obstáculo é mais comum (e muito confusa) chamada regressão beta com zero inflado. Pessoalmente, acho a última muito enganadora, porque a distribuição beta não possui zeros que poderiam ser inflados - mas, de qualquer maneira, é o termo padrão na literatura. Além disso, o modelo tobit é um modelo censurado e, portanto, não é um modelo de obstáculos. No entanto, poderia ser estendido por um modelo probit (ou logit) mais um modelo normal truncado . Na literatura econométrica, isso é conhecido como modelo de duas partes de Cragg.

Comentários de software

O countregpacote no R-Forge em https://R-Forge.R-project.org/R/?group_id=522 é a implementação sucessora de hurdle()/ zeroinfl()para pscl. A principal razão pela qual ainda não está no CRAN é que queremos revisar a predict()interface, possivelmente de uma maneira que não seja totalmente compatível com versões anteriores. Caso contrário, a implementação é bastante estável. Comparado a psclele, vem com alguns recursos interessantes, por exemplo:

Uma zerotrunc()função que usa exatamente o mesmo código hurdle()da parte truncada em zero do modelo. Assim, oferece uma alternativa para VGAM.
Além disso, ele funciona como d / p / q / r para as distribuições de contagem com zero truncado, obstáculo e zero. Isso facilita olhar para eles como um modelo "único" em vez de modelos separados.
Para avaliar a qualidade do ajuste, estão disponíveis exibições gráficas como rootogramas e parcelas residuais quantílicas aleatórias. (Ver Kleiber & Zeileis, 2016, The American Statistician , 70 (3), 296–303. Doi: 10.1080 / 00031305.2016.1173590 .)

Dados simulados

Seus dados simulados são provenientes de um único processo Poisson. Se efor tratado como um regressor conhecido, seria um Poisson GLM padrão. Se eé um componente de ruído desconhecido, então existe alguma heterogeneidade não observada causando um pouco de superdispersão que poderia ser capturado por um modelo binomial negativo ou algum outro tipo de mistura contínua ou efeito aleatório etc. No entanto, como o efeito de eé bastante pequena aqui , nada disso faz uma grande diferença. Abaixo, estou tratando ecomo um regressor (ou seja, com coeficiente verdadeiro de 1), mas você também pode omitir isso e usar binomial negativo ou modelos de Poisson. Qualitativamente, tudo isso leva a insights semelhantes.

## 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

Isso reflete que todos os três modelos podem estimar consistentemente os parâmetros verdadeiros. Observar os erros padrão correspondentes mostra que, neste cenário (sem a necessidade de uma parte do obstáculo), o Poisson GLM é mais eficiente:

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

Os critérios de informação padrão selecionariam o verdadeiro Poisson GLM como o melhor modelo:

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

E um teste de Wald detectaria corretamente que os dois componentes do modelo de barreira não são significativamente diferentes:

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

Finalmente ambos rootogram(p)e qqrplot(p)mostram que o de Poisson GLM se ajusta aos dados muito bem e que não há excesso de zeros ou sugestões em mais misspecifications.

Achim Zeileis
fonte

Qual é a diferença entre excesso de zeros e muitos zeros?

tatami

Um exemplo: Uma distribuição de Poisson com expectativa tem uma probabilidade de cerca de . Certamente são muitos zeros . No entanto, se você tiver uma distribuição com a forma de um Poisson (0,5), mas com mais zeros, haverá excesso de zeros .

λ = 0.5

$\lambda = 0.5$

f (0; λ = 0.5) \approx 60 %

$f(0; \lambda = 0.5) \approx 60\%$

Achim Zeileis

Concordo que é difícil entender a diferença entre os modelos inflado a zero e os obstáculos. Ambos são um tipo de modelo de mistura. Pelo que sei, a diferença importante é que, em um modelo inflado a zero, você mistura uma massa a zero com uma distribuição \ textit {que também pode assumir o valor zero}. Para um modelo de obstáculo, você mistura uma massa em zero com uma distribuição que aceita apenas valores maiores que 0. Portanto, no modelo inflado a zero, é possível distinguir entre 'zeros estruturais' (correspondente à massa em zero) e 'amostragem de zeros 'correspondente à ocorrência casual de 0 do modelo em que você está misturando. Naturalmente, essa identificação depende fortemente de fazer a escolha certa da distribuição! Mas, se você tem um Poisson inflado com zero, por exemplo, você pode distinguir entre zeros que vêm do componente Poisson (zeros de amostragem) e zeros que vêm da massa em zero (zeros estruturais). Se você tiver um modelo inflado a zero e a distribuição em que estiver misturando não tiver massa em zero, isso poderá ser interpretado como um modelo de barreira.

Glen Satten
fonte

Embora a distinção entre os dois tipos de zeros seja uma necessidade que sai diretamente da especificação do modelo, é possível calcular o mesmo tipo de quantidade para um modelo de barreira. Os chamados zeros estruturais também podem ser calculados a partir da distribuição de contagem não truncada (por exemplo, Poisson), embora seus parâmetros tenham sido baseados em uma amostra truncada . A probabilidade de zeros estruturais é então a diferença entre a probabilidade de zero (em geral, da parte do obstáculo zero) e de amostragem de zeros.

Achim Zeileis

Em relação ao aspecto filosófico, "quando devemos considerar algo como um modelo único e quando dois modelos separados" , pode ser interessante notar que as estimativas amostrais dos parâmetros-modelo estão correlacionadas.

No gráfico abaixo, com uma simulação, você vê principalmente a correlação entre a inclinação e a interceptação da parte de contagens. Mas há também uma ligeira relação entre a parte das contagens e a parte do obstáculo. Se você alterar os parâmetros, por exemplo, diminuir o lambda na distribuição de Poisson ou diminuir o tamanho da amostra, a correlação se tornará mais forte.

Então, eu diria que você não deve considerá-lo como dois modelos separados . Ou pelo menos há alguma relação, mesmo quando na prática você pode calcular as duas estimativas independentemente uma da outra.

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]))

Não faz muito sentido que haja uma correlação entre as duas partes. Mas isso provavelmente pode ser devido a níveis discretos das estimativas para os parâmetros no modelo de Poisson e como eles se relacionam com o número de zero.

Sextus Empiricus
fonte

Não consigo replicar isso. Para mim: truncdist::rtrunc(n = 100, spec = 'pois', a = 0, b = Inf, lambda = exp(-1.5 + rnorm(100)))produz um erro (usando a versão 1.0.2): Error in if (G.a == G.b) { : the condition has length > 1. De qualquer forma, o uso rhpois()do pacote countregno R-Forge é mais fácil para simular a partir de um modelo de Poisson de obstáculos com probabilidade de cruzamento de obstáculos pie expectativa de Poisson subjacente (não truncada) lambda. Se eu usá-las, obtenho apenas pequenas correlações empíricas entre as partes do obstáculo zero e da contagem truncada.

Achim Zeileis

Processo de geração de dados:

dgp <- function(n = 100, b = c(-0.5, 2), g = c(0.5, -2)) { x <- runif(n, -1, 1) ; y <- rhpois(n, lambda = exp(b[1] + b[2] * x), pi = plogis(g[1] + g[2] * x)); data.frame(x = x, y = y) }

Simulação: set.seed(1); cf <- t(replicate(3000, coef(hurdle(y ~ x, data = dgp())))). Avaliação: pairs(cf)e cor(cf). A verificação colMeans(cf)também mostra que a estimativa funcionou razoavelmente bem.

Achim Zeileis

@AchimZeileis no momento, não tenho possibilidade de analisar seu erro e comentar sobre ele. De qualquer forma, a correlação não é mais do que muito pequena na imagem que mostrei. A questão era mais filosófica / teórica. Na prática, você provavelmente terá poucos problemas ao tratar o modelo como duas etapas separadas e não integradas.

Sextus Empiricus