Variables de Conteo

Una variable de conteo nos da información sobre el número de eventos que ocurren en un periodo determinado de tiempo. Una característica importante de estas variables es que nunca pueden ser negativas y deben tomar un valor discreto (e.g., 0, 1, 2, 3, 4, …). Las variables dependientes de conteo por lo tanto tienden a tener un sesgo (hacia la izquierda o derecha) y una distribución de probabilidad discreta del tipo Poisson o Binomial Negativa.

Como un ejemplo de una variable de conteo vamos a usar el número de historias de televisión sobre política energética en un determinado mes. Esta variable proviene de la base de datos de Peake y Eshbaugh-Soha (2008) denotada como PESenergy.csv en la sección Files del Repositorio en Github de nuestra clase.

Vamos primero a darnos una idea de la distribución de la variable descrita anteriormente.

library(foreign)
energy_data <- read.csv("/Users/sergiobejar/Downloads/PESenergy.csv") ##carga base de datos

Con un histograma podemos ver la distribución de la variable dependiente.

library(ggplot2)
ggplot(energy_data, aes(x = Energy)) + 
        geom_histogram(binwidth = 10, col='black', fill='green', alpha=0.5) + 
        theme_bw() ## estoy agrupando en bins de tamaño 10 para propósitos ilustrativos

El número de noticias en un mes es un evento de conteo. Aquí hay que notar, sin embargo, que los datos mensuales son dependientes en el tiempo. Esta característica de nuestros datos la ignoraremos por el momento.

Regresión Poisson

El modelo más simple que podemos estimar cuando tenemos una variable dependiente de este tipo (conteo) es un modelo Poisson. Para ello tenemos que escribir lo siguiente:

m_poisson <- glm(Energy ~ rmn1173 + grf0175 + grf575 + jec477 + jec1177 + jec479 + embargo + hostages + oilc +
                  Approval + Unemploy, family = poisson(link = log),
                  data = energy_data) ## estima modelo Poisson

En este modelo la cobertura televisiva de política energética es una función de seis términos de discursos presidenciales (dummy que toma valor de 1 en el mes del discurso y 0 para otros meses), un indicador de embargo de petróleo Árabe (dummy), un indicador para la crisis de rehenes en Irán, el precio del petróleo, la aprobación presidencial y la tasa de desempleo.

Ahora vemos los resultados de la regresión.

summary(m_poisson) ## objeto con resultados de regresión
## 
## Call:
## glm(formula = Energy ~ rmn1173 + grf0175 + grf575 + jec477 + 
##     jec1177 + jec479 + embargo + hostages + oilc + Approval + 
##     Unemploy, family = poisson(link = log), data = energy_data)
## 
## Deviance Residuals: 
##    Min      1Q  Median      3Q     Max  
## -8.383  -2.994  -1.054   1.536  11.399  
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept) 13.250093   0.329121  40.259  < 2e-16 ***
## rmn1173      0.694714   0.077009   9.021  < 2e-16 ***
## grf0175      0.468294   0.096169   4.870 1.12e-06 ***
## grf575      -0.130568   0.162191  -0.805 0.420806    
## jec477       1.108520   0.122211   9.071  < 2e-16 ***
## jec1177      0.576779   0.155511   3.709 0.000208 ***
## jec479       1.076455   0.095066  11.323  < 2e-16 ***
## embargo      0.937796   0.051110  18.349  < 2e-16 ***
## hostages    -0.094507   0.046166  -2.047 0.040647 *  
## oilc        -0.213498   0.008052 -26.515  < 2e-16 ***
## Approval    -0.034096   0.001386 -24.599  < 2e-16 ***
## Unemploy    -0.090204   0.009678  -9.321  < 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: 6009.0  on 179  degrees of freedom
## Residual deviance: 2598.8  on 168  degrees of freedom
## AIC: 3488.3
## 
## Number of Fisher Scoring iterations: 5

Nótese que la función en este caso es lograrítmica (link = log) y aunque los coeficientes no son muy informativos, la interpretación es relativamente sencilla.

Iniciamos por estimar la razón de conteo:

exp(m_poisson$coefficients[-1]) ## estima razón de conteo para todos coefs. excepto intercepto
##   rmn1173   grf0175    grf575    jec477   jec1177    jec479   embargo  hostages 
## 2.0031352 1.5972665 0.8775970 3.0298703 1.7802943 2.9342587 2.5543461 0.9098211 
##      oilc  Approval  Unemploy 
## 0.8077536 0.9664787 0.9137448

Digamos que queremos interpretar el efecto de aprobación presidencial con estos resultados. La razón de conteo de la cobertura televisiva de política energética es 0.9664. Y esto lo podemos interpretar de la siguiente forma: un incremento porcentual de un punto en aprobación presidencial resulta en una disminución de 3.4% en promedio en cobertura televisiva de política energética manteniendo todos las otras variables constantes.

Una forma fácil de estimar el cambio porcentual para cada coeficiente es la siguiente:

100*(exp(m_poisson$coefficients[-1]) - 1) ## estima cambio porcentual para todos coefs. excepto intercepto
##    rmn1173    grf0175     grf575     jec477    jec1177     jec479    embargo 
## 100.313518  59.726654 -12.240295 202.987027  78.029428 193.425875 155.434606 
##   hostages       oilc   Approval   Unemploy 
##  -9.017887 -19.224639  -3.352127  -8.625516

Regresión Binomial Negativa

Una característica interesante de la distribución Poisson es que la varianza es igual a la media. Por lo que cuendo modelamos la media, nuestro modelo esta simultaneamente modelando la varianza.

Muchas veces nos podemos encontrar con variables de conteo en donde la varianza es más amplia de lo que esperamos -un fenómeno llamado sobredispersión. La regresión binomial negativa nos ofrece una solución a este fenómeno estimando un parámetro de dispersión extra que permite que la varianza condicional sea diferente a la mdia condicional.

Estimamos un modelo de regresión binomial negativa:

library(MASS) ## llamamos paquete MASS
m_bn <- glm.nb(Energy ~ rmn1173 + grf0175 + grf575 + jec477 + jec1177 + jec479 + embargo + hostages + oilc +
                  Approval + Unemploy,
                  data = energy_data) ## estima modelo Binomal Negativo

Imprimimos los resultados:

summary(m_bn)
## 
## Call:
## glm.nb(formula = Energy ~ rmn1173 + grf0175 + grf575 + jec477 + 
##     jec1177 + jec479 + embargo + hostages + oilc + Approval + 
##     Unemploy, data = energy_data, init.theta = 2.149960724, link = log)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.7702  -0.9635  -0.2624   0.3569   2.2034  
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept) 15.299318   1.291013  11.851  < 2e-16 ***
## rmn1173      0.722292   0.752005   0.960  0.33681    
## grf0175      0.288242   0.700429   0.412  0.68069    
## grf575      -0.227584   0.707969  -0.321  0.74786    
## jec477       0.965964   0.703611   1.373  0.16979    
## jec1177      0.573210   0.702534   0.816  0.41455    
## jec479       1.141528   0.694927   1.643  0.10045    
## embargo      1.140854   0.350077   3.259  0.00112 ** 
## hostages     0.089438   0.197520   0.453  0.65069    
## oilc        -0.276592   0.030104  -9.188  < 2e-16 ***
## Approval    -0.032082   0.005796  -5.536  3.1e-08 ***
## Unemploy    -0.077013   0.037630  -2.047  0.04070 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Negative Binomial(2.15) family taken to be 1)
## 
##     Null deviance: 393.02  on 179  degrees of freedom
## Residual deviance: 194.74  on 168  degrees of freedom
## AIC: 1526.4
## 
## Number of Fisher Scoring iterations: 1
## 
## 
##               Theta:  2.150 
##           Std. Err.:  0.242 
## 
##  2 x log-likelihood:  -1500.427

Los resultados de un modelo binomial negativo pueden ser interpretados de la misma manera que los resultados de una regresión Poisson.

100*(exp(m_bn$coefficients[-1]) - 1) ## estima cambio porcentual para todos coefs. excepto intercepto
##    rmn1173    grf0175     grf575     jec477    jec1177     jec479    embargo 
## 105.914750  33.408010 -20.354455 162.731901  77.395238 213.154943 212.943976 
##   hostages       oilc   Approval   Unemploy 
##   9.355996 -24.163611  -3.157261  -7.412208

Una diferencia importante, reportada en nuestro summary, es el parámetro de dispersión estimado (\(\theta\)). En este caso \(\theta\) = 2.15 con un error estándar de 0.242. Esto indica que sobredispersión está presente en nuestra variable dependiente. En efecto, hay diferencias importantes entre los resultados del modelo Poisson y el Binomial Negativo: varias variables que son significantes en el modelo Poisson no lo son el el Binomial Negativo. Además el AIC es substancialmente más bajo en el modelo negativo binomial, lo que indica un mejor fit aún penalizando con el paramétro de sobredispersión.

Gráfica Predicción de Eventos

Aunque la razón de conteo es una forma fácil de interpretar los coeficientes de los modelos, también podemos graficar los resultados. Hay que recordar que en este caso estamos modelando el logaritmo de la media parametral, así que tenemos que exponenciar nuestra predicción lineal para predecir la cuenta de eventos esperados dadas nuestras variables explicativas. Como en los modelos logit y probit, usamos la función predict.

Vamos a suponer que queremos graficar el efecto de aprobación presidencial en el número de historias de televisión sobre energía.

La variable approve tiene un rango de 24% de aprobación a 72.3%. Por lo tanto, primero construimos un vector que incluya el rango completo de aprobación, así como valores plausibles para todas las otras variables.

approval <- seq(24, 72.3, by = .1)
inputs.4 <- cbind(1, 0, 0, 0, 0, 0, 0, 0, 0, mean(energy_data$oilc), approval, mean(energy_data$Unemploy))
colnames(inputs.4) <- c("constant", "rmn1173", "grf0175", "grf575", "jec477", "jec1177", "jec479", "embargo", 
                        "hostages", "oilc", "Approval", "Unemploy")
inputs.4 <- as.data.frame(inputs.4)

Una vez que tenemos el data frame de las variables predictivas, podemos usar predict para predecir el número de historias televisivas sobre energía para los modelos Poisson y Binomial Negativo.

pred.poisson <- predict(m_poisson, newdata = inputs.4, type = "response")
pred.bn <- predict(m_bn, newdata = inputs.4, type = "response")

Y para graficar la predicciones para cada modelo:

plot(y = pred.poisson, x = approval, type = "l", lwd = 2, ylim = c(0, 60), 
       xlab = "Aprobación Presidencial", ylab = "Predicción de # Historias Sobre Energía")
lines(y = pred.bn, x = approval, lty = 2, col = "blue", lwd = 2)
legend(x = 50, y = 50, legend = c("Poisson", "Binomial Negativa"), lty = c(1, 2), col = c("black", "blue"), lwd = 2)