Modelo de regresión binomial negativa

Wed, Oct 7, 2020 9-minute read

Es muy común que la variable dependiente sea un número entero positivo o un conteo. Por ejemplo, el número de visitas al médico, el número de personas que tienen un infarto, el número fumado de cigarrillos, o el número de veces que asistes a tus clases virtuales. Estas variables, además de tener una función de distribución con masa de probabilidad concentrada en los enteros positivos, suelen estar sesgadas, presentar heteroscedasticidad y tener varianzas que se incrementan con su valor promedio.

En la literatura existen algunos modelos que nos permiten considerar variables de conteo. Entre estos se encuentran la regresión de Poisson y el modelo de regresión binomial negativa. En la regresión de Poisson, se supone que la variable dependiente, \(y_i\), sigue una distribución de Poisson con intensidad \(\mu_i\), donde:

\[\mu_i=\exp\{\mathbf{x}\boldsymbol{\beta}\},\]

y \(\boldsymbol{\beta}\) es, como siempre, un vector de coeficientes de dimensión \(k\times 1\). Dado que \(y_i\) sigue una distribución de Poisson tenemos que:

\[\begin{aligned} \mathbb{E}(y_i|\mathbf{x}) = \mu_i = \exp\{\mathbf{x}\boldsymbol{\beta}\} \\ \mathbb{V}(y_i|\mathbf{x}) = \mu_i = \exp\{\mathbf{x}\boldsymbol{\beta}\}. \end{aligned}\]

La limitación en la aplicabilidad del modelo de regresión de Poisson proviene de las restricciones impuestas por sus propios supuestos. La violación de tales supuestos implica la ausencia de equidispersión. A su vez, la violación del supuesto de equidispersión resulta en estimaciones ineficientes y en errores estándar sesgados.

El modelo de regresión binomial negativa es una extensión al modelo de regresión de Poisson que nos permite incorporar explícitamente la posibilidad de overdispersión. Para esto, consideramos una variable aleatoria, \(\xi_i\), tal que ahora:

\[\mu_i=\exp\{\mathbf{x}\boldsymbol{\beta} + \xi_i\}\]

Este nuevo elemento se justifica por la heterogeneidad de las observaciones. Esta heterogeneidad, por lo regular, no se observa (por ejemplo, diferencias en gustos y preferencias). Pero, estadísticamente hablando, ¿de dónde proviene \(\xi_i\)?

Regresión binomial negativa: teoría

Recuerda que, si una variable aleatoria \(y_i\) sigue una distribución binomial negativa con parámetros \(\mu\) y \(\alpha\), su función de densidad es:

\[ f(y_i;\mu, \alpha)=\frac{\Gamma(\alpha^{-1}+y_i)}{\Gamma(\alpha^{-1})\Gamma(y_i+1)}\left(\frac{\alpha^{-1}}{\alpha^{-1}+\mu}\right)^{\alpha^{-1}}\left(\frac{\mu}{\alpha^{-1}+\mu}\right)^{y_i}; \ \alpha>0 \]

donde \(\Gamma(\cdot)\) se refiere a la función gamma (la cual se simplifica al factorial en tanto su argumento sea un número entero). Con respecto a la distribución de Poisson, \(\alpha\) es un parámetro adicional y \(\mu\) tiene la misma interpretación.

Se puede demostrar que:

\[\begin{aligned} \mathbb{E}(y_i; \mu, \alpha) &= \mu \\ \mathbb{V}(y_i; \mu, \alpha) &= \mu + \alpha\mu^2 \end{aligned}\]

Nota que esta parametrización implica que la media condicional de una distribución binomial negativa es exactamente la misma que para una distribución de Poisson. Sin embargo, la varianza es más grande que la media (es decir, tenemos sobredispersión). Esto implica que la distribución binomial negativa introduce una mayor proporción de ceros y tiene una cola derecha más delgada.

Al igual que en modelo de regresión de Poisson, para derivar el modelo de regresión binomial negativa permitimos que el parámetro de intensidad varíe entre observaciones (por lo que tenemos \(\mu_i\) en lugar de la constante \(\mu\)). Este se usa para parametrizar la relación entre la variable dependiente \(y_i\) y el vector de regresores, \(\mathbf{x}\). En particular, suponemos que:

\[\mu_i=\exp\{\mathbf{x}\boldsymbol{\beta}\},\]

donde \(\boldsymbol{\beta}\) es, como siempre, un vector de coeficientes de dimensión \(k\times 1\). La diferencia entre el modelo de regresión de Poisson y el modelo de regresión binomial negativa se encuentra en la especificación de la varianza condicional de \(y_i\). En el modelo de regresión binomial negativa suponemos que:

\[ \mathbb{V}(y_i; \mu, \alpha) = \mu_i + \alpha\mu_i^2 \] Existe una segunda opción (menos natural y, por lo tanto, menos común en la práctica). Esta consiste en suponer que la varianza condicional de \(y_i\) no es una función \(\mu_i^2\), sino que estas guardan una relación lineal. Para esto, suponemos que \(\alpha=\delta/\mu_i\), por lo que:

\[ \mathbb{V}(y_i; \mu, \alpha) = \mu_i + \delta\mu_i \] Esta especificación lineal se conoce en la literatura como NB1 (de negative binomial 1). En cambio, la especificación cuadrática se conoce como NB2 (de negative binomial 2).

Para terminar esta sección teórica, es importante destacar que la distribución binomial negativa no es miembro de la familia exponencial, por lo que es muy sensible a los errores de especificación. Sin embargo, para valores fijos de \(\alpha\) (o \(\delta\)) es posible demostrar (Cameron y Trivedi, 2013) que los estimadores de un modelo de regresión binomial negativa son robustos a errores de especificación de momentos de orden mayor a 2.

Regresión binomial negativa: práctica

Para el siguiente ejemplo, estaremos utilizando los siguientes paquetes de R:

# Instalar paquetes si es necesario
# install.packages("tidyverse")
# install.packages("MASS")
# install.packages("margins")
# install.packages("countreg", repos="http://R-Forge.R-project.org")

# Leer paquetes
library(tidyverse)

## Warning: package 'tidyverse' was built under R version 4.0.2

## ── Attaching packages ─────────────────────────── tidyverse 1.3.0 ──

## ✓ ggplot2 3.3.2     ✓ purrr   0.3.4
## ✓ tibble  3.0.1     ✓ dplyr   1.0.0
## ✓ tidyr   1.1.0     ✓ stringr 1.4.0
## ✓ readr   1.3.1     ✓ forcats 0.5.0

## ── Conflicts ────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()

library(MASS)

## 
## Attaching package: 'MASS'

## The following object is masked from 'package:dplyr':
## 
##     select

library(margins)

## Warning: package 'margins' was built under R version 4.0.2

library(countreg)

Consideremos una muestra del Medical Expenditure Panel Survey 2004 de los Estados Unidos. Tómate un tiempo para analizar la base de datos, así como para realizar un análisis de estadística descriptiva.

db <- read_csv("https://raw.githubusercontent.com/amosino/courses--econometria/master/econometria_salud/econometria_salud--datos/MEPSData.csv")

## Parsed with column specification:
## cols(
##   .default = col_double(),
##   famidyr = col_character(),
##   hieuidx = col_character(),
##   female = col_character(),
##   race_bl = col_character(),
##   race_oth = col_character(),
##   eth_hisp = col_character(),
##   ed_hs = col_character(),
##   ed_hsplus = col_character(),
##   ed_col = col_character(),
##   reg_midw = col_character(),
##   reg_south = col_character(),
##   reg_west = col_character(),
##   anylim = col_character(),
##   ins_mcare = col_character(),
##   ins_mcaid = col_character(),
##   ins_unins = col_character(),
##   ins_dent = col_character()
## )

## See spec(...) for full column specifications.

Deseamos modelar el número de visitas al consultorio médico (use_off), la cual juega el papel de la variable dependiente y es, claramente, una variable de conteo, en función de la edad de la persona (age) y de su género (female). Esta es una variable dummy que vale Female si la observación es una mujer y Malesi es hombre. Aunque R no tiene ningún problema manejando variables dummy cuyas características estén codificadas como texto, podemos recodificar, si así lo deseamos, la variable female, tal que esta valga 1 si la observación es una mujer y 0 si es hombre.

db$female <- as.factor(as.numeric(db$female=="Female"))

Un histograma de la variable dependiente revela que esta está muy sesgada hacia la izquierda. Esto, y el hecho de que la variable dependiente sea discreta, justifican la aplicación de un modelo de conteo.

hist(db$use_off)

Ajustamos un modelo de regresión binomial negativa. Esto se hace con la función glm.nb(), la cual pertenece al paquete MASS. Esta función nos permite estimar el modelo en su versión NB2. (Nota. El lector interesado puede investigar sobre las funciones ml.b1() y ml.b2() del paquete COUNT.)

nbm.fit <- glm.nb(use_off~ age + female, data=db)
summary(nbm.fit)

## 
## Call:
## glm.nb(formula = use_off ~ age + female, data = db, init.theta = 0.5208532183, 
##     link = log)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.9452  -1.3161  -0.5186   0.0705   6.3859  
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept) 0.0802363  0.0317258   2.529   0.0114 *  
## age         0.0276920  0.0006026  45.952   <2e-16 ***
## female1     0.5192054  0.0212025  24.488   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Negative Binomial(0.5209) family taken to be 1)
## 
##     Null deviance: 23622  on 19385  degrees of freedom
## Residual deviance: 21108  on 19383  degrees of freedom
## AIC: 102672
## 
## Number of Fisher Scoring iterations: 1
## 
## 
##               Theta:  0.52085 
##           Std. Err.:  0.00632 
## 
##  2 x log-likelihood:  -102664.19500

La lectura de los resultados es idéntica que en el modelo de Poisson; esto es, los coeficientes pueden interpretarse como semielasticidades. En particular, podemos ver que número de visitas al médico aumenta en un 2.8% con cada año adicional de vida [\(100{\exp(0.0277)-1}\)] y las mujeres visitan al médico 52% [\(100{\exp(0.5192)-1}\)] más que los hombres.

También, resulta interesante calcular los efectos marginales e incrementales con respecto a \(y_i\) en lugar de los efectos marginales e incrementales con respecto a \(\ln(y_i)\) (que es lo que acabamos de hacer). Para esto, en R podemos utilizar la función margins() del paquete del mismo nombre.

summary(margins(nbm.fit))

##   factor    AME     SE       z      p  lower  upper
##      age 0.1645 0.0047 34.6647 0.0000 0.1552 0.1738
##  female1 2.9253 0.1226 23.8694 0.0000 2.6851 3.1655

La tabla anterior nos indica que el efecto marginal promedio para la variable age es 0.16. Es decir, en promedio, cada año adicional de vida implica visitar al médico 0.16 veces más. El efecto incremental de la variable female es 2.92. Esto quiere decir que las mujeres, en promedio, visitan al médico 2.92 veces más que los hombres.

También es posible evaluar los efectos marginales e incrementales en los valores promedio de las variables. Para esto, basta con agregar el argumento at a la función margins(). El siguiente bloque de código calcula el efecto marginal evaluando en la edad promedio y fijando el valor de la variable dummy en 1.

summary(margins(nbm.fit, at=list(age=mean(db$age)), variables="female"))

##   factor     age    AME     SE       z      p  lower  upper
##  female1 45.3609 2.5902 0.1082 23.9431 0.0000 2.3781 2.8022

Entonces, las mujeres de 45.36 años visitan al médico 2.59 veces en promedio.

Finalmente, nota que de la tabla de regresión de un modelo de regresión binomial negativo también podemos encontrar el valor del parámetro \(\alpha\). Este se calcula como \(\alpha= 1/\theta = 1/0.52085 = 1.9199\). (Nota. En muchos libros de texto, así como paquetes estadísticos, la \(\mathbb{V}(y_i; \mu, \theta)\) se expresa en función de \(\theta\) en lugar de \(\alpha\). Esto resulta intuitivo, puesto que queda claro que, conforme \(\theta \rightarrow \infty\), el modelo binomial negativo converge a uno de Poisson).

Grado de ajuste

Finalmente, podemos analizar el grado de ajuste de un modelo de regresión binomial negativo mediante un rootograma. Este representa la raiz cuadrada tanto de las frecuencias pronosticadas por el modelo como de las frecuencias empíricas (observadas) Las primeras se representan por puntos conectados por segmentos de línea y las segundas como un gráfico de barras.

Existen varias versiones de un rootograma tiene varias versiones. En su versión colgante, la barra cuelga de los puntos. Si la barra no llega al eje de las abscisas tenemos que la predicción de la frecuencia es superior a la frecuencia empírica. Por el contrario, si la barra corta el eje de las abscisas, el modelo predice una frecuencia inferior a la empírica. Naturalmente, el ajuste es mejor en tanto las barras justo toquen el eje de las abscisas.

En R podemos utilizar la función rootogram() que pertenece al paquete countreg. El uso de esta función es muy simple como se muestra a continuación.

rootogram(nbm.fit)

Como vemos, el ajuste de las frecuencias es bastante bueno.

Última actualización: 07-10-2020.

Econometría & Co.