yes

Regresion Poisson.

Objetivos de aprendizaje

Después de terminar este documento, debería poder:

  • Describa por qué la regresión lineal simple no es ideal para los datos de Poisson.

  • Escriba un modelo de regresión de Poisson e identifique los supuestos para la inferencia.

  • Escriba la probabilidad de una regresión de Poisson y describa cómo podría usarse para estimar los coeficientes de un modelo.

  • Interprete los coeficientes estimados de una regresión de Poisson y construya intervalos de confianza para ellos.

  • Utilice las desviaciones de los modelos de regresión de Poisson para comparar y evaluar modelos.

  • Utilice una compensación para tener en cuenta el esfuerzo variable en la recopilación de datos.

  • Ajuste y utilice un modelo de Poisson (ZIP) inflado a cero.

# Packages required 
library(gridExtra)
library(knitr)
library(kableExtra)
library(mosaic)
## Registered S3 method overwritten by 'mosaic':
##   method                           from   
##   fortify.SpatialPolygonsDataFrame ggplot2
## 
## The 'mosaic' package masks several functions from core packages in order to add 
## additional features.  The original behavior of these functions should not be affected by this.
## 
## Attaching package: 'mosaic'
## The following objects are masked from 'package:dplyr':
## 
##     count, do, tally
## The following object is masked from 'package:Matrix':
## 
##     mean
## The following object is masked from 'package:ggplot2':
## 
##     stat
## The following objects are masked from 'package:stats':
## 
##     binom.test, cor, cor.test, cov, fivenum, IQR, median, prop.test,
##     quantile, sd, t.test, var
## The following objects are masked from 'package:base':
## 
##     max, mean, min, prod, range, sample, sum
library(xtable)
library(pscl) 
## Classes and Methods for R developed in the
## Political Science Computational Laboratory
## Department of Political Science
## Stanford University
## Simon Jackman
## hurdle and zeroinfl functions by Achim Zeileis
library(multcomp)
## Loading required package: mvtnorm
## Loading required package: survival
## Loading required package: TH.data
## Loading required package: MASS
## 
## Attaching package: 'MASS'
## The following object is masked from 'package:dplyr':
## 
##     select
## 
## Attaching package: 'TH.data'
## The following object is masked from 'package:MASS':
## 
##     geyser
library(pander)
library(MASS)
library(tidyverse)
## -- Attaching packages --------------------------------------- tidyverse 1.3.1 --
## v tibble  3.1.1     v purrr   0.3.4
## v tidyr   1.1.3     v stringr 1.4.0
## v readr   1.4.0     v forcats 0.5.1
## -- Conflicts ------------------------------------------ tidyverse_conflicts() --
## x dplyr::combine()           masks gridExtra::combine()
## x mosaic::count()            masks dplyr::count()
## x purrr::cross()             masks mosaic::cross()
## x mosaic::do()               masks dplyr::do()
## x tidyr::expand()            masks Matrix::expand()
## x dplyr::filter()            masks stats::filter()
## x ggstance::geom_errorbarh() masks ggplot2::geom_errorbarh()
## x dplyr::group_rows()        masks kableExtra::group_rows()
## x dplyr::lag()               masks stats::lag()
## x tidyr::pack()              masks Matrix::pack()
## x MASS::select()             masks dplyr::select()
## x mosaic::stat()             masks ggplot2::stat()
## x mosaic::tally()            masks dplyr::tally()
## x tidyr::unpack()            masks Matrix::unpack()

Introducción a la regresión de Poisson

Considere las siguientes preguntas:

  1. ¿El número de muertes de motocicletas en un año determinado está relacionado con las leyes estatales sobre cascos?

  2. ¿Difiere el número de empleadores que realizan entrevistas en el campus durante un año para las universidades públicas y privadas?

  3. ¿Difiere el número diario de visitas relacionadas con el asma a una sala de emergencias según los índices de contaminación del aire?

  4. ¿Se ha visto afectado el número de peces deformados en los lagos de Minnesota seleccionados al azar por cambios en oligoelementos en el agua durante la última década?

Cada ejemplo implica predecir una respuesta utilizando una o más variables explicativas, aunque en estas los ejemplos tienen variables de respuesta que son recuentos por alguna unidad de tiempo o espacio. Un Poisson al azar la variable se utiliza a menudo para modelar recuentos; consulte para conocer las propiedades de la distribución de Poisson. Dado que una variable aleatoria de Poisson es un recuento, su valor mínimo es cero y, en teoría, el máximo es ilimitado. Nos gustaría modelar nuestro parámetro principal \(\lambda\), el número promedio de ocurrencias por unidad de tiempo o espacio, como un función de una o más covariables. Por ejemplo, en la primera pregunta anterior, \(\lambda_{i}\) representa el promedio número de muertes de motocicletas en un año para el estado \(i_{1}\) y esperamos mostrar que la variabilidad de un estado a otro en \(\lambda_{i}\) puede ser explicado por las leyes estatales sobre cascos Para un modelo de regresión lineal de mínimos cuadrados, el parámetro de interés es la respuesta promedio, \(\mu_{i}\), para asunto \(i_{1}\) y \(\mu_{i}\) se modela como una línea en el caso de una variable explicativa. Por analogía, podría parecer Es razonable intentar modelar el parámetro de Poisson \(\lambda_{i}\) como una función lineal de una variable explicativa, pero hay Hay algunos problemas con este enfoque. De hecho, un modelo como \(\lambda_{i} = \beta_{0} + \beta_{1} x_{i}\) no funciona bien para Poisson datos. Es seguro que una línea producirá valores negativos para ciertos \(x_{i,}\) pero \(\lambda_{i}\) solo puede tomar valores de 0 a \(\infty\). Además, el supuesto de igual varianza en la inferencia de regresión lineal se viola porque a medida que la tasa media para una variable de Poisson aumenta, la varianza también aumenta (recuerde del Capítulo 3 que si \(Y\) es el recuento observado, entonces \(E (Y) = var(Y) = \lambda)\).

Una forma de evitar estos problemas es modelar \(\log \left(\lambda_{i} \right)\) en lugar de \(\lambda_{i}\) como una función de las covariables. El \(\log\left(\lambda_{j} \right)\) toma valores de \(-\infty\) a \(\infty\). También podemos tener en cuenta el aumento de la varianza con un aumento de la media utilizando este enfoque. (Tenga en cuenta que en Beyond Multiple Linear Regression usamos log para representar el logaritmo natural.) Por lo tanto, consideraremos el modelo de regresión de Poisson: \[ \log \left(\lambda_{i} \right) = \beta_{0} + \beta_{1} x_{i} \] donde los valores observados \(Y_{i}\sim\)Poisson con \(\lambda = \lambda_{i}\) para un \(x{i}\) dado. Por ejemplo, cada estado \(i\) puede potencialmente tener un \(\lambda\) diferente dependiendo de su valor de \(x_{i}\), donde \(x_{i}\) podría representar la presencia o ausencia de una ley particular del casco. Tenga en cuenta que el modelo de regresión de Poisson no contiene un término de error separado como \(\epsilon\) we ver en regresión lineal, porque \(\lambda\) determina tanto la media como la varianza de una variable aleatoria Poisson.

Supuestos de regresión de Poisson

Al igual que la regresión lineal de mínimos cuadrados (LLSR), el uso de la regresión de Poisson para hacer inferencias requiere supuestos del modelo.

  1. Respuesta de Poisson La variable de respuesta es un recuento por unidad de tiempo o espacio, descrito por una distribución.

  2. Independencia Las observaciones deben ser independientes entre sí.

  3. Media = Varianza Por definición, la media de una variable aleatoria de Poisson debe ser igual a su varianza.

  4. Linealidad El logaritmo de la tasa media, \(\log(\lambda)\), debe ser una función lineal de \(x\).

# - UKaccident.csv is modified from builtin data Seatbelts
acc = read.csv("E:/UKaccident.csv")
#- driverskilled: number of death
#- law: before seatbelt law = 0, after law = 1
str(acc)
## 'data.frame':    122 obs. of  2 variables:
##  $ driverskilled: int  107 97 102 87 119 106 110 106 107 125 ...
##  $ law          : int  0 0 0 0 0 0 0 0 0 0 ...
head(acc); tail(acc)
# - some descriptives
tapply(acc$driverskilled, acc$law, sum)  # total death before vs after
##     0     1 
## 11826  1294
table(acc$law)  # num of observations before vs after
## 
##   0   1 
## 107  15
# - mean count, manually
11826/107  # 110.5234, count before law
## [1] 110.5234
1294/15  # 86.26667, count after law
## [1] 86.26667
model.acc = glm(driverskilled ~ law, data = acc, family = poisson)
summary(model.acc)  # significant p based on Wald test
## 
## Call:
## glm(formula = driverskilled ~ law, family = poisson, data = acc)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -3.16127  -0.72398   0.04531   0.77308   1.89182  
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  4.705227   0.009196 511.681   <2e-16 ***
## law         -0.247784   0.029281  -8.462   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 219.17  on 121  degrees of freedom
## Residual deviance: 142.64  on 120  degrees of freedom
## AIC: 940.7
## 
## Number of Fisher Scoring iterations: 4
# - to get CI
cbind(coef(model.acc), confint(model.acc))
## Waiting for profiling to be done...
##                             2.5 %     97.5 %
## (Intercept)  4.7052269  4.6871495  4.7231960
## law         -0.2477837 -0.3056189 -0.1908312
# - ln(count) = 4.71 - 0.25*LAW
4.71 - 0.25  # = 4.46
## [1] 4.46
exp(4.71)  # 111.0522, count before law
## [1] 111.0522
exp(4.46)  # 86.48751, count after law
## [1] 86.48751
# - Model fit
library(epiDisplay)
## Loading required package: foreign
## Loading required package: nnet
## 
## Attaching package: 'epiDisplay'
## The following object is masked from 'package:ggplot2':
## 
##     alpha
## The following object is masked from 'package:lattice':
## 
##     dotplot
poisgof(model.acc)  # fit well, based on chi-square test on the residual deviance
## $results
## [1] "Goodness-of-fit test for Poisson assumption"
## 
## $chisq
## [1] 142.6436
## 
## $df
## [1] 120
## 
## $p.value
## [1] 0.07764771
# - Diagnostics
#   - standardized residuals
sr = rstandard(model.acc)
sr[abs(sr) > 1.96]
##         4        54        55        91       113 
## -2.335861 -3.176147 -2.857937 -2.647896 -3.098644
#   - predicted count vs fitted values
fitted.acc = model.acc$fitted
data.frame(acc, fitted.acc)[names(sr[abs(sr) > 1.96]),]  # look at the discrepancies
# Summary with RR
idr.display(model.acc)  # easier, also view LR test
## 
## Poisson regression predicting driverskilled 
##  
##              IDR(95%CI)        P(Wald's test) P(LR-test)
## law: 1 vs 0  0.78 (0.74,0.83)  < 0.001        < 0.001   
##                                                         
## Log-likelihood = -468.3481
## No. of observations = 122
## AIC value = 940.6963

Regresion Poisson Ejemplo 2

dataset<-read.csv("https://stats.idre.ucla.edu/stat/data/poisson_sim.csv")

dataset$prog <- factor(dataset$prog,
  levels = 1:3,
  labels = c("General", "Academic", "Vocational"))

Estimando el modelo.

m1 <- glm(num_awards ~ prog + math, family = poisson, data = dataset)
m1
## 
## Call:  glm(formula = num_awards ~ prog + math, family = poisson, data = dataset)
## 
## Coefficients:
##    (Intercept)    progAcademic  progVocational            math  
##       -5.24712         1.08386         0.36981         0.07015  
## 
## Degrees of Freedom: 199 Total (i.e. Null);  196 Residual
## Null Deviance:       287.7 
## Residual Deviance: 189.4     AIC: 373.5
m2 <- glm(num_awards ~ math, family = poisson, data = dataset)
summary(m1)
## 
## Call:
## glm(formula = num_awards ~ prog + math, family = poisson, data = dataset)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.2043  -0.8436  -0.5106   0.2558   2.6796  
## 
## Coefficients:
##                Estimate Std. Error z value Pr(>|z|)    
## (Intercept)    -5.24712    0.65845  -7.969 1.60e-15 ***
## progAcademic    1.08386    0.35825   3.025  0.00248 ** 
## progVocational  0.36981    0.44107   0.838  0.40179    
## math            0.07015    0.01060   6.619 3.63e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 287.67  on 199  degrees of freedom
## Residual deviance: 189.45  on 196  degrees of freedom
## AIC: 373.5
## 
## Number of Fisher Scoring iterations: 6

La interpretación aquí es bastante sencilla. La tabla Coeficientes tiene la columna Estimación que nos muestra los coeficientes del Intercepto, del prog cuando es Académico, o Vocacional (no se muestra el General, ya que es implícitamente cero), y de Matemáticas.

Nota: todos los coeficientes están en la escala logarítmica (así que cuando dije “implícitamente cero”, en realidad es log (1) = 0).

Antes de continuar, y en aras de la exhaustividad, volvamos a poner estos coeficientes en la fórmula para que podamos comprenderla mejor.

Lo que dice, en palabras simples, es que si elevamos el número de Euler al resultado de sumar todos los coeficientes, obtendremos como salida el valor ajustado del modelo, verifique esto:

num_awards<-exp(-5.24712 + 1.08386 + 0.07015*dataset$math)
head(num_awards)
## [1] 0.2760675 0.2760675 0.3407317 0.2961290 0.2573650 0.2961290
# or: num_awards = e^(intercept_coef + progAcademic_coef + math_coef * math)
# or: log(num_awards) = intercept_coef + progAcademic_coef + math_coef * math
# or even: num_awards = exp(-5.24712) * exp(1.08386) * exp(0.07015)^math

Incrementemos ahora la robustez del resumen utilizando algunas funciones disponibles en el paquete sándwich.

# The default summary function assumes constant variances, thus
# computing (X'X)^(-1) from the QR decomposition of X.
# Here we recompute the SE assuming the variances are not constant, increasing
# robustness.
# The 'HC0' means we are using the White's estimator.
cov.m1 <- sandwich::vcovHC(m1, type = "HC0")
std.err <- sqrt(diag(cov.m1))
r.est <- cbind(
  Estimate = coef(m1), "Robust SE" = std.err,
  "Pr(>|z|)" = 2 * pnorm(abs(coef(m1) / std.err), lower.tail = FALSE),
  LL = coef(m1) - 1.96 * std.err,
  UL = coef(m1) + 1.96 * std.err
)
r.est
##                  Estimate  Robust SE     Pr(>|z|)          LL          UL
## (Intercept)    -5.2471244 0.64599839 4.566630e-16 -6.51328124 -3.98096756
## progAcademic    1.0838591 0.32104816 7.354745e-04  0.45460476  1.71311353
## progVocational  0.3698092 0.40041731 3.557157e-01 -0.41500870  1.15462716
## math            0.0701524 0.01043516 1.783975e-11  0.04969947  0.09060532

Aquí vemos algunas pequeñas variaciones en los resultados, pero lo más importante es que también tenemos el intervalo de confianza LL y UL.

Ratios de tasa de incidentes

A veces, queremos mostrar los coeficientes como razones de tasas de incidencia (las transformamos de la escala logarítmica) y sus errores estándar. Para lograr eso, usamos el método delta () del paquete msm

# drop the p-value and exp the coefficients
rexp.est <- exp(r.est[, -3])

# compute the SE for the exp coefficients
s <- msm::deltamethod(list(
  ~ exp(x1), ~ exp(x2), ~ exp(x3), ~ exp(x4)
), coef(m1), cov.m1)
# replace with the new SE's
rexp.est[, "Robust SE"] <- s
rexp.est[, "LL"] <- rexp.est[, 1] - 1.96 * s
rexp.est[, "UL"] <- rexp.est[, 1] + 1.96 * s

rexp.est
##                  Estimate  Robust SE           LL         UL
## (Intercept)    0.00526263 0.00339965 -0.001400685 0.01192594
## progAcademic   2.95606545 0.94903937  1.095948293 4.81618261
## progVocational 1.44745846 0.57958742  0.311467108 2.58344980
## math           1.07267164 0.01119351  1.050732371 1.09461091

Aquí vemos que para prog, usando el valor General como referencia, cuando el valor es Académico, la tasa de incidencia es 2.956 veces la referencia, y cuando el valor es Vocacional, la tasa de incidencia es 1.447 veces la referencia. Para la variable matemática continua, vemos que por cada aumento de unidad en matemáticas, la tasa de incidencia de num_awards aumenta en \(~ 7.3\%\).

Evaluación del modelo

Bondad de ajuste

Otro factor importante en todo modelo GLM a analizar es la relación entre la Desviación Residual (RD) y los Grados de Libertad (DF). El RD nos informa la bondad de ajuste para el modelo en su conjunto. El RD es la diferencia entre la desviación del modelo y la desviación de los datos (el mejor modelo teórico). Idealmente, el RD debería tender a cero. Podemos probar esto usando la prueba de chi-cuadrado. Si p no es significativo, aceptamos que el modelo se ajusta a los datos. Si no es así, podemos revisar nuestro modelo, buscando predictores omitidos, demasiados predictores, problemas de linealidad o dispersión excesiva (un problema común en las distribuciones de tipo “Poisson”).

# Check if the model fits the data, this is the null hypothesis.
with(
  m1,
  cbind(
    res.deviance = deviance,
    df = df.residual,
    p = pchisq(deviance, df.residual, lower.tail = FALSE)
  )
)
##      res.deviance  df         p
## [1,]     189.4496 196 0.6182274

Ahora, ¿recuerdas el segundo ajuste o modelo que hicimos? Usemos el ANOVA y verifiquemos qué modelo parece ser mejor para nuestros datos:

# Attention on the output: Model 1 is actually m2
anova(m2, m1, test = "Chisq")

Primero, atención, el Modelo 1 se refiere al modelo m2, el modelo que hicimos sin la variable prog.

Aquí vemos que el modelo ajustado que incluye la variable prog es un predictor significativamente mejor de num_awards.

Medias marginales

También podemos inspeccionar cuáles son los resultados esperados de nuestra variable dependiente num_awards para cada uno de los prog de predictores categóricos que mantienen matemáticas en su media global. Esto se llama los medias marginales. Para resolver esto, usamos nuestro modelo ajustado y hacemos predicciones con cada valor progresivo, manteniendo las matemáticas iguales a la media general:

# Create a dummy data frame with mean(math) and all prog levels.
(s1 <- data.frame(
  math = mean(dataset$math),
  prog = factor(1:3, levels = 1:3, labels = levels(dataset$prog))
))
# Then predict the values using type = "response".
(marginal <- predict(m1, s1, type = "response", se.fit = TRUE))
## $fit
##         1         2         3 
## 0.2114109 0.6249446 0.3060086 
## 
## $se.fit
##          1          2          3 
## 0.07050108 0.08628117 0.08833706 
## 
## $residual.scale
## [1] 1

Con esta salida, podemos desinfectar nuestros hallazgos en Razones de tasas de incidentes: manteniendo las matemáticas en su media, el número o eventos predichos para progGeneral (1) es 0.211, para progAcademic (2) es 0.625 y para progVocational (3) es 0.306. Por lo tanto, al hacer las proporciones obtenemos progAcademic / progGeneral igual a 2.956 y progAcademic / progVocational igual a 1.447. Las mismas proporciones estimadas.

Mejor que simplemente mirar un número es mirar un gráfico. Así que construyamos una combinación de los valores predichos de nuestro modelo junto con las medias marginales.

# backup the dataset, so we don't mess with the original
data <- dataset
# compute the fitted values
data$phat <- predict(m1, type = "response")

# order by program and then by math
data <- data %>% arrange(prog, math)

# create the plot using math as X axis, the marginal means as Y and prog as colors
ggplot(data, aes(x = math, y = phat, color = prog)) +
  # plots the line of the marginal means
  geom_line(size = 1) +
  # plots the fitted values, we use 'jitter' to make it prettier
  geom_point(aes(y = num_awards), alpha = .5, position = position_jitter(h = .2)) +
  labs(title = "Fitted + Marginal means", x = "Math Score", y = "Expected number of awards")

A tener en cuenta

  • Cuando los datos parecen estar dispersos en exceso, debemos mirar si nuestro modelo es realmente el modelo apropiado. Es posible que necesitemos agregar o eliminar variables, por ejemplo. Arriba vimos cómo la eliminación de prog creaba un modelo peor.

  • Suponiendo que el modelo se especifica correctamente, se debe verificar el supuesto principal de igual media y varianza. Hay varias pruebas para eso, fuera del alcance de este artículo.

  • El problema más común que podríamos enfrentar es tener un proceso de generación de datos adicional que contamina nuestros datos. Y para Poisson, es común tener un proceso de generación cero (distribución de Poisson más muchos ceros). En ese caso tenemos el modelo inflado a cero que es más apropiado. Si los datos no permiten ceros, use el modelo truncado por cero.

  • Recuerde verificar si sus datos tienen una variable de exposición (como eventos y tamaño_población). Esta variable debe estar presente en el modelo con la opción de compensación.

  • La variable de resultado en una regresión de Poisson no puede tener números negativos y la exposición no puede tener ceros.

  • La regresión de Poisson se estima mediante la estimación de máxima verosimilitud. Por lo general, requiere un tamaño de muestra grande.