5 Modelos lineales generalizados

El término “Modelo Lineal Generalizado” (MLG o GLM, por su nombre en inglés) se refiere a una clase más amplia de modelos popularizados por McCullagh y Nelder (1982). Esta clase de models incluyen casos como: la regresión lineal, el análisis de la varianza, los modelos logit y probit, los modelos de respuesta multinomial, entre otros.

Se agrupan debido a la cantidad de propiedades que tienen en común, como la linealidad o los métodos de estimación de sus parámetros. Por tanto, tiene sentido examinar esta diversidad de modelos dentro de una categoría única, la cual se denomina Modelos Lineales Generalizados.

Dado que es una extensión de los modelos lineales, algunas hipótesis de los modelos originales se conservan. Una de ellas es la necesidad de observaciones independientes o sin correlación, es decir, que el valor de cada uno de los datos no está influenciado de ninguna manera con el resto de los datos de la muestra.

No obstante, en los MLG, no es necesario imponer otras condiciones más restrictivas como la normalidad. Este hecho amplía considerablemente la aplicabilidad de estos modelos, lo que habilita el estudio de conjuntos de datos con distribuciones no gaussianas.

5.1 Estructura de los modelos lineales generalizados

Un modelo lineal generalizado tiene tres componentes básicos:

Componente aleatoria: Identifica la variable respuesta y su distribución de probabilidad.
Componente sistemática: Especifica las variables explicativas (independientes o predictoras) utilizadas en la función predictora lineal.
Función liga (o enlace -link-): Es una función del valor esperado de $Y$ como una combinación lineal de las variables predictoras.

5.1.1 Componente aleatoria

La componente aleatoria de un MLG consiste en una variable aleatoria $\boldsymbol{Y}$ con observaciones independientes $\left(Y_1, \ldots, Y_n\right)$.

En estos modelos, se asume que la variable respuesta sigue una distribución de la familia exponencial.

5.1.1.1 Familia Exponencial

La familia exponencial contiene como casos especiales la mayoría de las distribuciones discretas y continuas estándar que utilizamos para modelado práctico, como

Discretas	Continuas
Bernoulli	Normal
Binomial (con $n$ conocida)	Normal multivariada
Multinomial (con $n$ conocida)	Exponencial
Binomial negativa (con $n$ conocida)	Gamma
Geometrica	Ji-cuadrada
Poisson	Beta
	Dirichlet
	Wishart
	Wishart inversa

La razón del estatus especial de la familia exponencial es que varias cálculos importantes y útiles en estadística pueden realizarse de manera uniforme dentro del marco de esta familia.

Definición 5.1 (Familia exponencial) Sea $\boldsymbol{X}={X_1,\ldots,X_n}$ una m.a. con distribución $F_\theta$, donde $\boldsymbol{\theta} \in \Theta \subset \mathbb{R}^k$. La familia $\left\{f(\cdot\ ;\ \theta)\right\}$ pertenece a la familia exponencia si su f.d. (f.d.p.) se puede escribir como \[ f_{\boldsymbol{X}}(\boldsymbol{x} ;\ \boldsymbol{\theta} )= h(\boldsymbol{x}) \exp\left\{\sum_{i=1}^k \eta_i(\boldsymbol{\theta} )\ T_i(\boldsymbol{x})-\psi(\boldsymbol{\theta} )\right\}, \]

donde $\eta_i:\mathbb{R}^k \rightarrow \mathbb{R}$, $T_i:\mathbb{R}^n \rightarrow \mathbb{R}$, $i=1,\ldots,k$, $\psi: \mathbb{R}^k \rightarrow \mathbb{R}$ y $h: \mathbb{R}^n \rightarrow \mathbb{R}^+$ y además se cumple que la dimensión de $\Theta$ y la de la imagen de $\Theta$ bajo el mapeo \[ \left(\theta_1, \theta_2, \cdots, \theta_k\right) \longrightarrow\left(\eta_1\left(\theta_1, \theta_2, \cdots, \theta_k\right), \eta_2\left(\theta_1, \theta_2, \cdots, \theta_k\right), \cdots, \eta_k\left(\theta_1, \theta_2, \cdots, \theta_k\right)\right) \] son iguales a $k$.

En la familia exponencial, el vector de estadísticos $(T_1,\ldots, T_k)$ se llama el estadístico suficiente natural.

Ejemplos: Distribuciones que pertenecen a la familia exponencial.

Distribución Binomial. Sea $X \sim \operatorname{Bin}(n, p)$, con $n$ conocida. \[ \begin{gathered} f_X(x ;\ p)=\binom{n}{x} p^x(1-p)^{n-x}\ \mathbb{I}_{\{x \in\{0,1, \cdots, n\}\}}=\binom{n}{x}\left(\frac{p}{1-p}\right)^x(1-p)^n \ \mathbb{I}_{\{x \in\{0,1, \cdots, n\}\}} \\ =\binom{n}{x} \exp\left\{x \log \frac{p}{1-p}+n \log (1-p)\right\}\ \mathbb{I}_{\{x \in\{0,1, \cdots, n\}\}}. \end{gathered} \] Definiendo \[ \eta(p)=\log \frac{p}{1-p},\quad T(x)=x,\quad \psi(p)=-n \log (1-p), \quad \text{y}\quad h(x)=\binom{n}{x}\, \mathbb{I}_{\{x \in\{0,1, \cdots, n\}\}}, \] tenemos la representación de la familia exponencial para $p \in(0,1)$ para una sola observación. Claramente para $p=0$ o 1, la distribución se degenera.

Basta con obtener la expresión para una sola observación ya que si $\boldsymbol{X}={X_1,\ldots,X_n}$ es una m.a. entonces

\[ f_{\boldsymbol{X}}(\boldsymbol{x} ;\ \boldsymbol{\theta} ) = h(x)^n\exp\left\{ \sum_{i=1}^k\eta_i(\boldsymbol{\theta} ) \sum_{j=1}^n T_i(\boldsymbol{x_j})-n\psi(\boldsymbol{\theta} ) \right\}. \]

Distribución Gaussiana. Sea $X \sim N\left(\mu, \sigma^2\right)$, denotamos $(\mu, \sigma)=\left(\theta_1, \theta_2\right)=\boldsymbol{\theta}$, entonces \[ \begin{align} f(x ;\ \boldsymbol{\theta}) &= \frac{1}{\sqrt{2 \pi} \theta_2} \exp\left\{-\frac{\left(x-\theta_1\right)^2}{2 \theta_2^2}\right\} \mathbb{I}_{x \in \mathbb{R}}\\ &=\frac{1}{\sqrt{2 \pi} \theta_2} \exp\left\{-\frac{x^2}{2 \theta_2^2}+\frac{\theta_1 x}{\theta_2^2}-\frac{\theta_1^2}{2 \theta_2^2}\right\} \mathbb{I}_{x \in \mathbb{R}} \end{align} \] Que es miembro de la familia exponencial con \[ \begin{align} \eta_1(\boldsymbol{\theta})=-\frac{1}{2 \theta_2^2},\qquad \eta_2(\boldsymbol{\theta})=\frac{\theta_1}{\theta_2^2}&,\qquad T_1(x)=x^2,\qquad T_2(x)=x \\ \psi(\boldsymbol{\theta})=\frac{\theta_1^2}{2 \theta_2^2}+\log \theta_2&,\qquad h(x)=\frac{1}{\sqrt{2 \pi}} \mathbb{I}_{x \in \mathbb{R}}. \end{align} \] El número de parámetros $\boldsymbol{\theta}$ que varían libremente es el mismo que el número de elementos que varían libremente en $\boldsymbol{\eta}$, y ambos son iguales a $k=2$.
Distribución Gamma. Sea $X \sim \operatorname{Gam}(\alpha,\lambda)$ con f.d. \[ f_X(x ;\ \alpha,\lambda)=\frac{x^{\alpha-1}\exp\left\{-\frac{x}{\lambda} \right\}}{\lambda^\alpha \Gamma(\alpha)} \mathbb{I}_{x>0}. \] Sea $\boldsymbol{\theta}=(\alpha, \lambda)=\left(\theta_1, \theta_2\right)$, así \[ f_X(x ;\ \boldsymbol{\theta})=\exp\left\{-\frac{x}{\theta_2}+\theta_1 \log x-\theta_1 \log \theta_2-\log \Gamma\left(\theta_1\right)\right\} \frac{1}{x} \mathbb{I}_{x>0}, \] que es la forma de la familia exponencial de dos parámetros con \[ \begin{align} \eta_1(\boldsymbol{\theta})=-\frac{1}{\theta_2},\qquad \eta_2(\boldsymbol{\theta})=\theta_1&,\qquad T_1(x)=x,\qquad T_2(x)=\log x, \\ \psi(\boldsymbol{\theta})=\theta_1 \log \theta_2+\log \Gamma\left(\theta_1\right)&,\qquad h(x)=\frac{1}{x} \mathbb{I}_{x>0}. \end{align} \]

Los términos forma canónica, parámetro natural y espacio de parámetros naturales tendrán el mismo significado que en el caso de un solo parámetro. Así, si parametrizamos las distribuciones utilizando $\eta_1, \eta_2, \cdots, \eta_k$ como los $k$ parámetros, entonces el vector $\eta=\left(\eta_1, \eta_2, \cdots, \eta_k\right)$ se denomina vector de parámetros naturales. La parametrización \[ f(x ;\ \boldsymbol{\eta})=h(x)\exp\left\{\sum_{i=1}^k \eta_i T_i(x)-\psi(\eta)\right\} \]

se llama forma canónica, y el conjunto de todos los vectores $\boldsymbol{\eta}$ para los cuales $f(x ;\ \boldsymbol{\eta})$ es una densidad válida (o función de masa de probabilidad, pmf) se denomina espacio de parámetros naturales.

5.1.1.2 Variable respuesta

Como hemos comentado, en los MLG la variable respuesta $Y_i$ sigue una distribución de la familia exponencial.

En muchas aplicaciones, las observaciones

son binarias y se identifican como éxito y fracaso.
indican el número de éxitos de entre un número fijo de ensayos, y se modeliza como una distribución binomial.
cuentan el número de fracasos antes de obtener una cierta cantidad de éxitos.
cuentan el número de eventos que ocurren en un lapso de tiempo o espacio
cuentan el número de eventos de una categoría, de entre $m$ diferente categorías, donde el total es conocido.
Finalmente, si las observaciones son continuas se puede asumir para $Y$ una distribución normal. Así que el modelo de regresión ordinario es también un MLG.

5.1.2 Componente Sistemática

La componente sistemática de un MLG especifica las variables explicativas, que entran en forma de efectos fijos en un modelo lineal, es decir, las covariables variables $x_j$ que se relacionan mediante \[ \beta_0+\beta_1 x_1+\cdots+\beta_k x_k \]

La combinación lineal de variables explicativas se denomina predictor lineal.

Como hemos comentado, en los MLG la variable respuesta $y_i$ sigue una distribución de la familia exponencial. Además suponemos que su media $\mu_i$, es alguna función de las covariables. Aunque este modelo puede ser lineal como el de arriba algunos llamarían a estos modelos “no lineales” porque $\mu_i$ puede no se iguala directamente con la expresión lineal anteriorir, sino con una una función no lineal de ésta. Sin embargo, McCullagh y Nelder los siguen denominando modelo lineales porque las covariables afectan la distribución de $y_i$ sólo a través del predictor lineal $\boldsymbol{x}_i^{\top} \boldsymbol{\beta}$.

Considerando $n$ observaciones, denotamos el vector de predictores lineales como $\left(\eta_1, \ldots, \eta_n\right)$, donde \[ \eta_i=\sum_j \beta_j x_{i j} \]

El término independiente $\beta_0$ se obtienen con esta notación haciendo que todos los $x_{i j}$ sean igual a 1 para todos los $i$.

En cualquier caso, se pueden considerar variables que estén basadas en otras variables. Por ejemplo, $x_3=x_1 x_2$ o $x_3=x_2^2$, para modelizar interacciones entre variables o efectos curvilíneos de $x_2$.

5.1.3 Función liga

Se denota el valor esperado de $Y$ como $\mu=\mathbb{E}(Y)$, entonces la función link especifica una función $g(\cdot)$ que relaciona $\mu$ con el predictor lineal como \[ g(\mu)=\beta_0+\beta_1 x_1+\cdots+\beta_k x_k \]

Así, la función link $g(\cdot)$ relaciona las componentes aleatoria y sistemática. De este modo, para $i=1, \ldots, n$, \[ g\left(\mathbb{E}\left[Y_i\right]\right)=g\left(\mu_i\right)=\eta_i=\sum_j \beta_j x_{j i} \]

La función $g$ más simple es la identidad, $g(\mu)=\mu$, que da lugar al modelo de regresión lineal ordinario.

Los MLG permiten la unificación de una amplia variedad de métodos estadísticos como la regresión, los modelos ANOVA y los modelos de datos categóricos. En realidad se usa el mismo algoritmo para obtener los estimadores de máxima verosimilitud en todos los casos. Este algoritmo es la base del procedimiento GENMOD de SAS y de la función glm de R.

Para ilustrar la utilizad de las funciones liga, consideremos el caso en el que $Y$ produce datos binarios.

En muchos casos las respuestas tienen solo dos categorías del tipo sí/no (éxito/fracaso 1/0) que naturalmente lleva a considerar un modelo donde $Y \sim \operatorname{Bin}(1, \pi)$. En este caso \[ \begin{aligned} f(y ;\ \pi) & =\pi^y(1-\pi)^{1-y} \ \mathbb{I}_{y=\{0,1\}}\\ & =(1-\pi)\left(\frac{\pi}{1-\pi}\right)^y \ \mathbb{I}_{y=\{0,1\}}\\ & =\exp \left\{y \log \left(\frac{\pi}{1-\pi}\right)+\log(1-\pi)\right\} \ \mathbb{I}_{y=\{0,1\}} \end{aligned} \tag{5.1}\]

Vemos que el parámetro natural es \[ Q(\pi)=\log \left(\frac{\pi}{1-\pi}\right)=:\operatorname{logit}(\pi) \] y la media es \[ \mathbb{E}(Y)=\mathbb{P}(Y=1)=\pi \] que lleva al modelo \[ \pi_i=\pi(\boldsymbol{x}_i)=\boldsymbol{x}_i^{\top}\boldsymbol{\beta}. \] En respuestas binarias, un modelo análogo al de regresión lineal es \[ \pi(\boldsymbol{x}_i)=\boldsymbol{x}_i^{\top}\boldsymbol{\beta}, \] que se denomina modelo de probabilidad lineal, ya que la probabilidad de éxito cambia linealmente con respecto a $\boldsymbol{x}$.

El parámetro $\beta_j$ representa el cambio en probabilidad por unidad de $x_{j i}$. Este modelo es un MLG con un componente aleatorio binomial y con función link igual a la identidad.

Sin embargo, este modelo tiene el problema de que aunque las probabilidades deben estar entre 0 y 1, el modelo puede predecir a veces valores $\pi(\boldsymbol{x})>1$ o $\pi(\boldsymbol{x})<0$.

Para evitar este problema se puede utiizar la función $\text{logit}:[0,1]\rightarrow \mathbb{R}$ como función liga.

Así, \[ \operatorname{logit}\left(\pi(\boldsymbol{x}_i)\right)=\boldsymbol{x}_i^{\top}\boldsymbol{\beta}. \] que equivale a \[ \pi(\boldsymbol{x}_i)=\text{logit}^{-1}(\boldsymbol{x}_i^{\top}\boldsymbol{\beta})=\frac{\exp\{\boldsymbol{x}_i^{\top}\boldsymbol{\beta}\}}{1+\exp\{\boldsymbol{x}_i^{\top}\boldsymbol{\beta}\}}. \] A la funcion $\text{logit}^{-1}(\cdot)$ se le conoce como función expit.

Algunos otros ejemplos de MLG se presentan en la siguiente tabla, antes de proceder a abordar algunos de ellos con mayor detalle.

Distribución	Función de Enlace Común	Propósito/Uso	Fórmula
Normal	Identidad	$Y$ toma valores continuos. Relaciona directamente el predictor lineal con la variable respuesta.	$\eta=\mu$
Binomial	Logit	$Y$ dicotómica. Transforma la probabilidad de éxito a una escala no acotada.	$\eta=\log \left(\frac{\mu}{1-\mu}\right)$
Poisson	Log	$Y$ de conteo. Relaciona el logaritmo de la cuenta media con los predictores lineales.	$\eta=\log (\mu)$
Gamma	Inversa	$Y$ toma valores en $\mathbb{R}^+$ como el tiempo hasta un evento. Modela la inversa de la variable respuesta.	$\eta=\frac{1}{\mu}$
Multinomial	Softmax (Logit para múltiples categorías)	$Y$ toma valores categoricos (más de dos tipos). Generaliza la función logit para resultados multiclase.	$\eta_i=\log \left(\frac{e^{x_i}}{\sum e^{X_j}}\right)$, para la $i$-ésima categoría
Exponencial	Log	$Y>0$ como tiempo hasta un evento. Aplicable en casos de procesos de Poisson.	$\eta=\log (\mu)$
Inv. Gaussiana	Inv. al cuadrado	$Y$ continua. Usada para modelar datos con sesgo positivo.	$\eta=\frac{1}{\mu^2}$

Las funciones liga no son únicas para cada distribución y en particular, en el caso de los datos binarios la función logit no es la única liga popular. A continuación se presentan otras dos funciones liga para comparar:

liga probit: $\Phi^{-1}(\cdot):[0,1]\rightarrow \mathbb{R}$, donde $\Phi$ es la función de distribución estándar.
liga cloglog: $\log(-\log(1-\pi))$

Código

library(ggplot2)

x <- seq(-5, 5, by = 0.1)

# valores evaluados de las funciones logit y probit 
logit_vals <- plogis(x)
probit_vals <- pnorm(x)

# calculamos los valores de la función cloglog
cloglog_vals <- 1 - exp(-exp(x))

# Creamos data frame
data <- data.frame(x, logit_vals, probit_vals, cloglog_vals)

ggplot(data, aes(x = x)) +
  geom_line(aes(y = logit_vals, color = "Logit"), size = 1) +
  geom_line(aes(y = probit_vals, color = "Probit"), size = 1) +
  geom_line(aes(y = cloglog_vals, color = "CLogLog"), size = 1) +
  labs(title = "Inversa de Funciones liga: Logit, Probit y CLogLog",
       x = "x", y = "Probability") +
  scale_color_manual(values = c("Logit" = "red", "Probit" = "blue", "CLogLog" = "olivedrab3")) +
  theme_minimal(base_size = 15)

Figura 5.1: Visualización de tres funciones liga

La elección de la función liga es parte del modelo y otra función liga que puede tomarse en cuenta es en realidad una famiia de funciones, denominada familia de enlaces exponencial, la cual se define como: \[ {\eta}= \begin{cases}{\mu}^\lambda, & \text { si } \lambda \neq 0 \\ \log ({\mu}), & \text { si } \lambda=0\end{cases} \] La selección de una función liga (selección de $\lambda$) puede compararse a la elección de una transformación en la variable respuesta. Sin embargo, es importante destacar que la función liga afecta exclusivamente a la media de la variable respuesta, a diferencia de las transformaciones usuales que influyen directamente en la función de distribución. Al igual que con las transformaciones, optar por una función liga inadecuada puede generar problemas en la estimación de los parámetros del modelo.

5.2 Casos específicos

5.2.1 Regresión con variable binomial

Consideramos variables respuesta de tipo binario, es decir, para cada individuo, la variable respuesta para cada uno de los sujetos puede tomar únicamente dos posibles valores y \[ \mathbb{P}(Y_i=0)=1−\pi_i,\qquad \mathbb{P}(Y_i=1)=\pi_i \] Si además se han observado una serie de covariables continuas o de tipo factor (categoricas), el objetivo del análisis con el modelo lineal generalizado será predecir la probabilidad asociada al éxito, $\widehat{\pi}_i$ en función de dichas covariables.

Existen dos formas habituales en las que se presenta la información experimental sobre variables que miden éxito o fracaso:

Información individualizada por sujeto. En este caso la base de datos recoge la información de cada sujeto de la muestra con una variable que indica el éxito o el fracaso. Dicha variable se suele codificar 1 (éxito) y 0 (fracaso). Asociada con esta variable binaria se pueden recoger variables predictoras, en cada uno de los sujetos, para tratar de explicar el comportamiento de la respuesta.
Información agrupada por sujetos. La base de datos de este tipo suelen identificar por fila el “tratamiento” al que se ven sometidos un grupo de sujetos, registrándose el número total de sujetos en esa combinación y el número de éxitos (en algunos casos se recogen los éxitos y fracasos para dicha combinación). El “tratamiento” puede ser la combinación de una o más predictoras. Estos bancos de datos son muy habituales en ensayos de dosis-repuesta donde únicamente queremos valorar eficacia, ya que no es necesario recoger una información exhaustiva sobre los sujetos bajo estudio. En estos casos no estimamos la probabilidad de éxito individual sino del grupo de sujetos que se ven sometidos al mismo tratamiento. La información muestral debe recoger la probabilidad de éxito asociada a cada combinación.

Las hipótesis que debe verificar los modelos de este tipo son:

independencia entre las observaciones
linealidad entre transformaciones de la proporción de éxitos y de las variables explicativas continuas (función link)

Además de la verificación de las hipótesis el modelo también se evalua considerando la consistencia entre los resultados del modelo y la interpretación física.

Si $\pi_i$ es la probabilidad de éxito asociado al i-ésimo experimento binomial con respuesta $Y_i$, y condiciones de experimentación observadas en las variables predictoras $x_1, x_2, \ldots, x_k$, el MLG se establece como: \[ g\left(\pi_i\right)=\eta_i=\beta_0+\beta_1 x_{i 1}+\beta_2 x_{i2}+\ldots+\beta_k x_{i k}, \qquad i=1,\ldots,n. \] donde $g(\cdot)$ es la función liga que en anteriormente hemos planteado puede ser la función logit, pero no es la única que puede utilizarse. A continuación, presentamos algunas de las más comunes funciones liga para este tipo de datos.

5.2.1.1 Regresión logit

La función liga logit proporciona los comúnmente conocidos como modelos de regresión logística: \[ g\left(\pi_i\right)=\log \left(\frac{\pi_i}{1-\pi_i}\right)=\beta_0+\beta_1 x_{i 1}+\beta_2 x_{i2}+\ldots+\beta_k x_{i k} \] Bajo este modelo la probabilidad de éxito se puede escribir en términos de las variables predictoras (despejando de la ecuación anterior) como: \[ \pi_i=\frac{\exp \left(\beta_0+\beta_1 x_{i 1}+\beta_2 x_{i2}+\ldots+\beta_k x_{i k}\right)}{1+\exp \left(\beta_0+\beta_1 x_{i 1}+\beta_2 x_{i2}+\ldots+\beta_k x_{i k}\right)}=\frac{\exp(\boldsymbol{x}_i^{\top}\boldsymbol{\beta})}{1+\exp(\boldsymbol{x}_i^{\top}\boldsymbol{\beta})}=\frac{1}{1+\exp(-\boldsymbol{x}_i^{\top}\boldsymbol{\beta})} \] Si $\pi\left(\boldsymbol{x}_1\right)$ y $\pi\left(\boldsymbol{x}_2\right)$ denotan las probabilidades de éxito para dos valores de la predictora ( $\boldsymbol{x}_1$ y $\boldsymbol{x}_2$ ), pero consideremos la $i$-ésima variable explicativa y $\pi\left(x_{1r}\right)$ y $\pi\left(x_{2r}\right)$ las probabilidades con $\boldsymbol{x}$ fija, exceptuando en el elemento de la $r$-ésima variable para la cual se tienen dos valores diferentes. Entonce el conciente \[ \operatorname{logit}(\pi\left(x_{1r}\right))-\operatorname{logit}(\pi\left(x_{2r}\right))=\log\left[\frac{\pi\left(x_{1r}\right)}{ 1 - \pi\left(x_{1r}\right)}\div\frac{\pi\left(x_{2r}\right)}{ 1-\pi\left(x_{2r}\right)}\right]=(x_{1r}-x_{2r})\beta_r \] se denomina tasa de momios y compara los momios bajo dos escenarios de valores de la covariable $r$, de forma que:

un valor mayor que 0 indica un aumento de la probabilidad de éxito en $x_{1r}$ con respecto a $x_{2r}$ y naturalmente,
un valor menor que 0 indica un decremento de la probabilidad de éxito en $x_{1r}$ con respecto a $x_{2r}$.

En términos del modelo ajustado se puede valorar el incremento o decremento de los momios con el valor de $\exp \left(\beta_r\right)$, para cada uno de los efectos (variables numéricas o niveles de un factor) presentes en el modelo. Si $x_{1r}>x_{2r}$, valores negativos de $\beta_r$ dan lugar a odds ratios menores que 0, mientras que valores de $\beta_r$ positivos dan odds ratios positivos.

De hecho, esta propiedad es la que motiva que esta función de enlace sea la más utilizada en este tipo de modelos.

Los modelos de regresión logística contribuyen entonces a:

Identificar factores de riesgo.
Evaluar el riego (probabilidad de ocurrencia del evento respuesta) para individuos específicos.
Clasificar / discriminar a grupos de individuos como alto riesgo o bajo riesgo.

Clasificamos a un individuo como de alto riesgo si la probabilidad de ocurrencia es mayor o igual a cierto umbral; entonces, diremos que un individuo con atributos $\boldsymbol{x}$ es de alto riesgo si \[ \frac{1}{1+\exp\left\{-\boldsymbol{x}^{\top}\boldsymbol{\beta}\right\}} > p, \] esto es, si \[ \boldsymbol{x}^{\top} \boldsymbol{\beta} \; > \; \text{log}\left(\frac{p}{1-p}\right). \] En el caso particular de $p=0.5$ tenemos que los individuos de alto riesgo satisfacen $\boldsymbol{x}^{\top}\boldsymbol{\beta} \; > \; 0$.

Verosimilitud, logverosimilitud y función Score

La verosimilitud de los datos independientes $(x_1^{\top},y_1)$, $\cdots,(x_n^{\top},y_n)$ es \[ L(\boldsymbol{\beta}) = \prod_{i=1}^n \pi_i^{y_i} (1-\pi_i)^{1-y_i}, \qquad \text{con} \quad \pi_i = \frac{1}{1+\exp(-\boldsymbol{x}_i^{\top}\boldsymbol{\beta})} \] Estimamos $\boldsymbol{\beta}$ como aquel valor que maximiza la logverosimilitud \[ \ell(\boldsymbol{\beta}) =\log L(\boldsymbol{\beta})= \sum_{i=1}^n \left[ \, y_i \log(\pi_i) + (1-y_i) \log(1-\pi_i) \, \right]. \] Derivando la logverosimilitud, tenemos $S(\boldsymbol{\beta})$: \[ \begin{align} S(\boldsymbol{\beta})=\frac{\partial \ell(\boldsymbol{\beta})}{\partial \boldsymbol{\beta}} &= \sum_{i=1}^n \left[ \; y_i \frac{1}{\pi_i} \frac{\partial \pi_i}{\partial \boldsymbol{\beta}} - (1-y_i) \frac{1}{1-\pi_i} \frac{\partial \pi_i}{\partial \boldsymbol{\beta}} \; \right] \\ &= \sum_{i=1}^n \left[ \; \frac{ y_i}{\pi_i} - \frac{1-y_i}{1-\pi_i} \; \right] \frac{\partial \pi_i}{\partial \boldsymbol{\beta}} \\ &= \sum_{i=1}^n \left[ \; \frac{ y_i - \pi_i }{\pi_i(1-\pi_i)} \; \right] \frac{\partial \pi_i}{\partial \boldsymbol{\beta}}. \end{align} \] Por otro lado, \[ \begin{align} \frac{\partial \pi_i}{\partial \boldsymbol{\beta}} = & \; -\frac{\exp\left\{-\boldsymbol{x}_i^{\top}\boldsymbol{\beta}\right\} (-\boldsymbol{x}_i) }{\left( \; 1+\exp\left\{-\boldsymbol{x}_i^{\top}\boldsymbol{\beta}\right\} \; \right)^{2}} \\ = & \; \frac{1}{1+\exp\left\{-\boldsymbol{x}_i^{\top}\boldsymbol{\beta}\right\}} \; \frac{\exp\left\{-\boldsymbol{x}_i^{\top}\boldsymbol{\beta}\right\}}{1+\exp\left\{-\boldsymbol{x}_i^{\top}\boldsymbol{\beta}\right\}} \; \boldsymbol{x}_i \; = \; \pi_i(1-\pi_i) \boldsymbol{x}_i, \end{align} \] entonces \[ S(\boldsymbol{\beta})=\frac{\partial \ell(\boldsymbol{\beta})}{\partial \boldsymbol{\beta}} \; = \; \sum_{i=1}^n \left[ \; \frac{ y_i - \pi_i }{\pi_i(1-\pi_i)} \; \right] \pi_i(1-\pi_i) \boldsymbol{x}_i \; = \; \sum_{i=1}^n ( y_i - \pi_i ) \; \boldsymbol{x}_i. \]

Así, las ecuaciones de verosimilitud son: \[ \frac{\partial \ell(\boldsymbol{\beta})}{\partial \boldsymbol{\beta}} \; = \; \sum_{i=1}^n ( y_i - \pi_i ) \; \boldsymbol{x}_i = 0. \] Este es un sistema de $p$ ecuaciones nolineales en $n$ incógnitas que puede resolverse con métodos numéricos como el método de Newton, que es un procedimiento iterativo.

Método de Newton

Las ecuaciones de verosimilitud son: \[ \frac{\partial \ell(\boldsymbol{\beta})}{\partial \boldsymbol{\beta}} \; = \; \sum_{i=1}^n ( y_i - \pi_i ) \; \boldsymbol{x}_i \; = \; 0. \]

Para el caso de una sola covariable, se ven como: \[ \begin{align} \frac{\partial \ell(\beta)}{\partial \boldsymbol{\beta}} = & \; \sum_{i=1}^n ( y_i - \pi_i ) \; \left[ \begin{array}{c} 1 \\ x_i \end{array} \right] = \left[ \begin{array}{c} \sum (y_i - \pi_i) \\ \sum ( y_i - \pi_i ) x_i \end{array} \right] \\[5pt] = \; : &\; \; \left[ \begin{array}{c} h_1(\boldsymbol{\beta}) \\ h_2(\boldsymbol{\beta}) \end{array} \right] = S(\boldsymbol{\beta}) = 0, \quad (\text{donde } \; S: \mathbb{R}^2 \rightarrow \mathbb{R}^2) \end{align} \]

Para resolver $S(\boldsymbol{\beta})=0$, hacemos una aproximación de primer orden para $S$ alrededor de algún valor inicial razonable $\boldsymbol{\beta}^{(0)}$: \[ S(\boldsymbol{\beta}) \approx S(\boldsymbol{\beta}^{(0)}) + \frac{\partial S(\boldsymbol{\beta}^{(0)})}{\partial \boldsymbol{\beta}} (\boldsymbol{\beta} - \boldsymbol{\beta}^{(0)}). \]

Para el caso de una sola covariable, \[ \frac{\partial S(\boldsymbol{\beta}^{(0)})}{\partial \boldsymbol{\beta}} = \left[\begin{array}{c} \frac{\partial h_1(\boldsymbol{\beta}^{(0)}) }{ \partial \boldsymbol{\beta}^{\top}} \\ \frac{\partial h_2(\boldsymbol{\beta}^{(0)}) }{ \partial \boldsymbol{\beta}^{\top}} \end{array} \right] \qquad \text{es una matriz} \; 2 \times 2. \]

Entonces, en vez de resolver $S(\boldsymbol{\beta})=\boldsymbol{0}$, resolvemos el problema más fácil \[ S(\boldsymbol{\beta}^{(0)}) + \frac{\partial S(\boldsymbol{\beta}^{(0)})}{\partial \boldsymbol{\beta}} (\boldsymbol{\beta} - \boldsymbol{\beta}^{(0)})=0, \]

así, despejando para $\boldsymbol{\beta}$, obtenemos \[ \boldsymbol{\beta}^{(1)} = \boldsymbol{\beta}^{(0)} - \left[ \; \frac{\partial S(\boldsymbol{\beta}^{(0)})}{\partial \boldsymbol{\beta}} \; \right]^{-1} S(\boldsymbol{\beta}^{(0)}). \] Esta expresión la iteramos hasta convergencia \[ \boldsymbol{\beta}^{(k+1)} = \boldsymbol{\beta}^{(k)} - \left[ \; \frac{\partial S(\boldsymbol{\beta}^{(k)})}{\partial \boldsymbol{\beta}} \; \right]^{-1} S(\boldsymbol{\beta}^{(k)}). \]

Las diferentes partes que intervienen en la expresión anterior son: \[ \begin{align} S(\boldsymbol{\beta}) &= \sum ( y_i - \pi_i ) \boldsymbol{x}_i \\[7pt] \frac{\partial S(\boldsymbol{\beta})}{\partial \boldsymbol{\beta}^{\top}} &= \frac{\partial^2 \ell(\boldsymbol{\beta})}{\partial \boldsymbol{\beta}\partial \boldsymbol{\beta^{\top}}} -\sum_{i=1}^n \frac{\partial \pi_i}{d\beta} \boldsymbol{x}_i^{\top} -\sum_{i=1}^n \pi_i(1-\pi_i) \boldsymbol{x}_i \boldsymbol{x}_i^{\top} \end{align} \]

Reescribimos en notación matricial. Definamos $X$ y $W$ de tamaños $n \times p$ y $n \times n$, respectivamente, mientras $y$ y $p$ son vectores de $n \times 1$: \[ X = \left[ \begin{array}{c} \boldsymbol{x}_1^{\top} \\ \vdots \\ \boldsymbol{x}_n^{\top} \end{array} \right] \quad W = \left[ \begin{array}{ccc} \pi_1(1-\pi_1) & \cdots & 0 \\ \vdots & \ddots & \vdots \\ 0 & \cdots & \pi_n(1-\pi_n) \end{array} \right] \] \[ \boldsymbol{Y} = \left[ \begin{array}{c} y_1 \\ \vdots \\ y_n \end{array} \right] \quad \text{y} \quad \boldsymbol{\pi} = \left[ \begin{array}{c} \pi_1 \\ \vdots \\ \pi_n \end{array} \right] \] con esto, es fácil ver que \[ \boldsymbol{\beta}^{(k+1)} \; = \; \boldsymbol{\beta}^{(k)} - \left[ \; \frac{\partial S(\boldsymbol{\beta}^{(k)} )}{\partial \boldsymbol{\beta}} \; \right]^{-1} S(\boldsymbol{\beta}^{(k)} ) \; = \; \boldsymbol{\beta}^{(k)} + (X^{\top}WX)^{-1}X^{\top}(\boldsymbol{Y}-\boldsymbol{\pi}). \]

Mínimos Cuadrados Reponderados Iterativamente (IRWLS)

En la literatura de Modelos Lineales Generalizados encontramos que la estimación se hace comunmente mediante “Mínimos Cuadrados Reponderados Iterativamente” (Iteratively Reweighted Least Squares). Este método es precisamente el que acabamos de ver.

Si definimos el vector de “observaciones de trabajo” \[ \widetilde{y\ } = X\boldsymbol{\beta}^{(k)} + W^{-1}(\boldsymbol{Y}-\boldsymbol{\pi}) \] entonces, el Método de Newton es \[ \begin{align} \boldsymbol{\beta}^{(k+1)} = & \; \boldsymbol{\beta}^{(k)} + (X^TWX)^{-1}X^{\top}(\boldsymbol{Y}-\boldsymbol{\pi}) \\ = & \; (X^{\top}WX)^{-1}X^{\top}WX \boldsymbol{\beta}^{(k)} +(X^{\top}WX)^{-1}X^{\top}WW^{-1}(\boldsymbol{Y}-\boldsymbol{\pi}) \\ = & \; (X^{\top}WX)^{-1}X^{\top}W [\; X \boldsymbol{\beta}^{(k)} + W^{-1}(\boldsymbol{Y}-\boldsymbol{\pi}) \; ] \\ = & \; (X^{\top}WX)^{-1}X^{\top}W \widetilde{y\ }. \end{align} \]

Esto es, $\boldsymbol{\beta}^{(k+1)} = (X^TWX)^{-1}X^TW \widetilde{y\ }$, y de aquí es de donde le viene el nombre de “mínimos cuadrados ponderados”.

Implementación con R

El siguiente ejemplo considera datos de Edad y enfermedad coronaria incluido en Hosmer, D.W. & Lemeshow, S. (1989) Applied logistic regression. Wiley.

Ejemplo: Hosmer & Lemeshow, S. (1989)

Código

edad <- c(
 20, 23, 24, 25, 25, 26, 26, 28, 28, 29, 30, 30, 30, 30, 30, 30, 32, 32, 33, 33,
 34, 34, 34, 34, 34, 35, 35, 36, 36, 36, 37, 37, 37, 38, 38, 39, 39, 40, 40, 41, 
 41, 42, 42, 42, 42, 43, 43, 43, 44, 44, 44, 44, 45, 45, 46, 46, 47, 47, 47, 48,
 48, 48, 49, 49, 49, 50, 50, 51, 52, 52, 53, 53, 54, 55, 55, 55, 56, 56, 56, 57,
 57, 57, 57, 57, 57, 58, 58, 58, 59, 59, 60, 60, 61, 62, 62, 63, 64, 64, 65, 69)

coro <- c(
 0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,1,0,0,0,0,0,1,0,0,1,0,0,0,0,1,0,1,0,
 0,0,0,0,1,0,0,1,0,0,1,1,0,1,0,1,0,0,1,0,1,1,0,0,1,0,1,0,0,1,1,1,1,0,1,1,1,1,1,0,
 0,1,1,1,1,0,1,1,1,1,0,1,1,1,1,1,0,1,1,1)
 

# Gráfica de los datos
edaj <- jitter(edad)    # solo con fines de graficación

plot(edaj, coro, xlab="Edad", ylab="Indicador CHD", ylim=c(-.1, 1.1),
 mgp=c(1.5,.5,0), cex.axis=.8, cex.lab=.8, cex.main=1, xlim=c(15,75), cex=1.2,
 main="Regresión lineal", pch=19,col=rgb(.5,0,.5,.5))
rug(edaj)
out <- lm(coro ~ edad)
abline(out,lwd=2,col="blue")
grid()

Código

plot(edaj, coro, xlab="Edad", ylab="Indicador CHD", ylim=c(-.1, 1.1),
 mgp=c(1.5,.5,0), cex.axis=.8, cex.lab=.8, cex.main=1, xlim=c(15,75), cex=1.2,
 main="Regresión lineal", pch=19,col=rgb(.5,0,.5,.5))
rug(edaj)
# Resolviendo ecuaciones de verosimilitud
y       <- coro
n       <- length(y)
X       <- cbind(rep(1,n),edad)
b       <- c(-10,.2) # valores iniciales
# Las 4 líneas anteriores son específicas para los datos de coronaria
tolm    <- 1e-6      # tolerancia (norma minima de delta)
iterm   <- 100       # numero maximo de iteraciones
tolera  <- 1         # inicializar tolera
itera   <- 0         # inicializar itera
histo   <- b         # inicializar historial de iteraciones

while((tolera > tolm) & (itera < iterm)){
  p      <- 1/( 1+exp( -as.vector(X%*%b) ) )
  W      <- p*(1-p)
  delta  <- as.vector( solve(t(X*W)%*%X, t(X)%*%(y-p)) )
  b      <- b + delta
  tolera <- sqrt( sum(delta*delta) )
  histo  <- rbind(histo,b)
  itera  <- itera + 1 
} 

head(histo)

#>             [,1]       [,2]
#> histo -10.000000 0.20000000
#> b      -1.497206 0.03767488
#> b      -4.380358 0.09253679
#> b      -5.221685 0.10918224
#> b      -5.308597 0.11090419
#> b      -5.309453 0.11092114

Código

tail(histo)

#>        [,1]       [,2]
#> b -1.497206 0.03767488
#> b -4.380358 0.09253679
#> b -5.221685 0.10918224
#> b -5.308597 0.11090419
#> b -5.309453 0.11092114
#> b -5.309453 0.11092114

Código

# Agregamos curva logística a la gráfica original
xx   <- seq(15,75,length=200)
X    <- cbind(rep(1,n),xx)
p    <- 1/( 1+exp( -as.vector(X%*%b) ) )
lines(xx,p,lwd=2,col="blue")
grid()

Usamos la función glm:

Código

dat <- data.frame(edad=edad, cor=coro)
fit_Bin <- glm(coro ~ edad, data = dat, family = "binomial")
coefficients(fit_Bin)

#> (Intercept)        edad 
#>  -5.3094534   0.1109211

Errores estándar

Una forma aproximada de calcular la varianza del estimador de máxima verosimilitud es mediante la función de información observada: \[ \operatorname{Var}(\widehat{\boldsymbol{\beta}}) \doteq \left[ -\frac{\partial^2 \ell(\widehat{\boldsymbol{\beta}}) }{\partial \boldsymbol{\beta} \partial \boldsymbol{\beta}^{\top} } \right]^{-1}, \] esto es, la varianza (asintótica) de $\widehat{\boldsymbol{\beta}}$ se aproxima por el inverso de la Matriz observada de Información.

Interpretación heurística: La logverosimilitud es globalmente cóncava (para los modelos lineales generalizados), por lo tanto su segunda derivada es negativa, además, la segunda derivada mide el grado de curvatura de la logverosimilitud. Mientras más grande sea la segunda derivada, más “picuda” es la logverosimilitud y mejor definido está su máximo.

Entonces, si tomamos el recíproco del negativo de la curvatura tenemos una medida de que tan mal está nuestro estimador, (i.e. que tanta varianza tiene), a mayor curvatura menor varianza.

Para el caso de regresión logística vimos que \[ \begin{align} \frac{\partial \ell(\boldsymbol{\beta})}{\partial \boldsymbol{\beta}} \; = & \; S(\boldsymbol{\beta}) \; = \; \sum_{i=1}^n (y_i-\pi_i)\boldsymbol{x}_i \\ \frac{\partial^2 \ell(\boldsymbol{\beta})}{\partial \boldsymbol{\beta} \partial \boldsymbol{\beta}^{\top}} \; = & \; \frac{\partial S(\boldsymbol{\beta}) }{\partial \boldsymbol{\beta}^{\top}} \; = \; -\sum_{i=1}^n \pi_i (1-\pi_i) \boldsymbol{x}_i \boldsymbol{x}_i^{\top} \; = \; -X^{\top}WX \end{align} \] así que, la varianza se puede aproximar por: \[ \operatorname{Var}(\widehat{\boldsymbol{\beta}}) = \left[\, X^{\top} \widehat{W} X \, \right]^{-1}. \]

Los errores estándar de los estimadores se obtienen sacando raíz cuadrada a los elementos diagonales de $\operatorname{Var}(\widehat{\boldsymbol{\beta}})$.

Ejemplo: Hosmer & Lemeshow, S.(1989) - Cont.

Código

summary(fit_Bin)

#> 
#> Call:
#> glm(formula = coro ~ edad, family = "binomial", data = dat)
#> 
#> Coefficients:
#>             Estimate Std. Error z value Pr(>|z|)    
#> (Intercept) -5.30945    1.13365  -4.683 2.82e-06 ***
#> edad         0.11092    0.02406   4.610 4.02e-06 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> (Dispersion parameter for binomial family taken to be 1)
#> 
#>     Null deviance: 136.66  on 99  degrees of freedom
#> Residual deviance: 107.35  on 98  degrees of freedom
#> AIC: 111.35
#> 
#> Number of Fisher Scoring iterations: 4

5.2.1.2 Regresión probit

El link probit proporciona los comúnmente conocidos como modelos de regresión probit y su expresión viene dada por:

\[ g\left(\pi_i\right)=\Phi^{-1}\left(\pi_i\right)=\beta_0+\beta_1 x_{i 1}+\beta_2 x_{i2}+\ldots+\beta_k x_{i k} \] donde $\Phi^{-1}$ es la función inversa de la función de distribución Normal estándar. En esta situación la probabilidad de éxito se puede escribir en términos de las variables predictoras (despejando de la ecuación anterior) como: \[ \pi_i=\Phi\left(\beta_0+\beta_1 x_{i 1}+\beta_2 x_{i2}+\ldots+\beta_k x_{i k}\right) \]

donde $\Phi$ es la función de distribución Normal estándar.

Ejemplo: probit

Código

fit_BinProbit <- glm(coro ~ edad, data = dat, family = binomial(link = "probit"))

summary(fit_BinProbit)

#> 
#> Call:
#> glm(formula = coro ~ edad, family = binomial(link = "probit"), 
#>     data = dat)
#> 
#> Coefficients:
#>             Estimate Std. Error z value Pr(>|z|)    
#> (Intercept) -3.14573    0.62460  -5.036 4.74e-07 ***
#> edad         0.06580    0.01335   4.930 8.20e-07 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> (Dispersion parameter for binomial family taken to be 1)
#> 
#>     Null deviance: 136.66  on 99  degrees of freedom
#> Residual deviance: 107.50  on 98  degrees of freedom
#> AIC: 111.5
#> 
#> Number of Fisher Scoring iterations: 4

5.2.1.3 Regresión cloglog

La liga cloglog (‘complementary log log’) es el menos habitual en la práctica y su expresión viene dada por:

\[ g\left(\pi_i\right)=\log \left(-\log \left(1-\pi_i\right)\right)=\beta_0+\beta_1 x_{i 1}+\beta_2 x_{i2}+\ldots+\beta_k x_{i k}. \]

Ejemplo: cloglog

Código

fit_BinCloglog <- glm(coro ~ edad, data = dat, family = binomial(link = "cloglog"))

summary(fit_BinCloglog)

#> 
#> Call:
#> glm(formula = coro ~ edad, family = binomial(link = "cloglog"), 
#>     data = dat)
#> 
#> Coefficients:
#>             Estimate Std. Error z value Pr(>|z|)    
#> (Intercept) -4.24889    0.83108  -5.112 3.18e-07 ***
#> edad         0.07918    0.01630   4.858 1.18e-06 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> (Dispersion parameter for binomial family taken to be 1)
#> 
#>     Null deviance: 136.66  on 99  degrees of freedom
#> Residual deviance: 107.55  on 98  degrees of freedom
#> AIC: 111.55
#> 
#> Number of Fisher Scoring iterations: 5

¿Cuál es la probabilidad de éxito en este caso?

La diferencia entre una función liga y otra hace referencia a como modelamos las probabilidades más extremas, es decir, las muy bajas o muy altas, pero en la práctica proporcionan resultados muy similares. Por eso casi siempre el modelo utilizado es el de regresión logística.

Distribuciones de los coeficientes ajustados

La distribución de $\widehat{\boldsymbol{\beta}}$ se obtiene a partir de la teoría asintótica de los EMV: \[ \widehat{\boldsymbol{\beta}} \stackrel{D}{\rightarrow} N_{p}\left(\boldsymbol{\beta}, \mathcal{I}(\boldsymbol{\beta})^{-1}\right) \quad \text{cuando} \quad n \rightarrow \infty, \tag{5.2}\]

donde $\mathcal{I}(\boldsymbol{\beta})$ es la función de Información de Fisher: \[ \mathcal{I}(\boldsymbol{\beta}):=-\mathbb{E}\left[\frac{\partial^2 \ell(\boldsymbol{\beta})}{\partial \boldsymbol{\beta} \partial \boldsymbol{\beta}^{\top}}\right] \] Su nombre proviene del hecho de que mide la información disponible en la muestra para estimar $\boldsymbol{\beta}$. Cuanto “más grande” sea la matriz, más precisa será la estimación de $\boldsymbol{\beta}$, ya que esto resulta en varianzas más pequeñas en (5.2).

El inverso de la matriz de información de Fisher, para el caso binomial con liga logit, es: \[ \mathcal{I}(\boldsymbol{\beta})^{-1}=\left[\, X^{\top} W X \, \right]^{-1} \] En el caso de la regresión lineal múltiple, $\mathcal{I}(\boldsymbol{\beta})^{-1}=\sigma^2\left(\mathbf{X}^{\prime} \mathbf{X}\right)^{-1}$, por lo que la presencia de $W$ aquí es una consecuencia de la heterocedasticidad del modelo logístico.

De estos resultados vemos algunos elementos que inciden en la calidad de la estimación:

Sesgo: Las estimaciones son asintóticamente insesgadas.
Varianza: Depende de:
- Tamaño de muestra $n$: Oculto dentro de $X^{\top} W X$. A medida que $n$ crece, aumenta la precisión de los estimadores.
- Dispersión ponderada de los predictores $\left(X^{\top} W X\right)^{-1}$: Cuanto más “dispersos” sean los predictores, más precisa será $\widehat{\boldsymbol{\beta}}$. Cuando $p=1$, $X^{\top} W X$ es una versión ponderada de $s_x^2$.

La precisión de $\widehat{\boldsymbol{\beta}}$ está influenciada por el valor de $\boldsymbol{\beta}$, que está oculto dentro de $W$. Esto contrasta fuertemente con el modelo lineal, donde la precisión del estimador de mínimos cuadrados no se veía afectada por $\boldsymbol{\beta}$. Esto se debe en parte a la heterocedasticidad de la regresión logística, lo que implica una dependencia de la varianza de $Y$ en la curva logística y, por ende, en $\boldsymbol{\beta}$.

Al igual que en la regresión lineal, el problema con las expresiones de varianza asintótica es que $W$ es desconocida en la práctica, ya que depende de $\boldsymbol{\beta}$. Al sustituir la estimación $\widehat{\boldsymbol{\beta}}$ en $W$, se obtiene el estimador $\widehat{W}$. Ahora podemos usar $\widehat{W}$ para obtener: \[ \frac{\widehat{\beta}_j-\beta_j}{\widehat{\mathrm{SE}}\left(\widehat{\beta}_j\right)} \stackrel{D}{\sim} N(0,1), \qquad \text{con} \quad \widehat{\mathrm{SE}}\left(\widehat{\beta}_j\right)^2:=v_j \tag{5.3}\]

donde $v_j$ es el $j$-ésimo elemento de la diagonal de $\left(X^{\top} \widehat{W} X\right)^{-1}$.

Intervalos de confianza para $\boldsymbol{\beta}$

Gracias a la ecuación (5.3), podemos obtener el intervalo de confianza (IC) del $(1-\alpha)\times 100\%$ para el coeficiente $\beta_j$, $j=0, \ldots, k$:

\[ \left(\widehat{\beta}_j \ \pm \ \widehat{\mathrm{SE}}\left(\widehat{\beta}_j\right) z_{\alpha / 2}\right) \]

donde $z_{\alpha / 2}$ es el cuantil superior $\alpha / 2$ de la distribución $N(0,1)$.

Si estamos interesados en el IC de $\exp{\beta_j}$, simplemente podemos tomar la exponencial del intervalo de confianza anterior. Así, el IC del$(1-\alpha)\times 100\%$ para $\exp^{\beta_j}$, $j=0, \ldots, p$ es: \[ \exp\left\{\widehat{\beta}_j \ \pm \ \widehat{\mathrm{SE}}\left(\widehat{\beta}_j\right) z_{\alpha / 2}\right\}. \]

5.2.1.4 Más sobre Regresión logit

5.2.1.4.1 Prueba de hipótesis para los coeficientes

Las distribuciones asintótica del EMV $\widehat{\boldsymbol{\beta}}$ también nos permiten realizar una prueba formal de hipótesis sobre los coeficientes $\beta_j$, con $j=0, \ldots, p$. Por ejemplo, la prueba de significancia: \[ H_0: \beta_j=0. \]

Como hemos visto, esta prueba particularmente interesante, ya que nos permite determinar si la variable $X_j$ tiene un efecto significativo en $Y$.

La estadística utilizada para probar la significancia es la estadística de Wald: \[ \frac{\widehat{\beta}_j-0}{\widehat{\operatorname{SE}}\left(\widehat{\beta}_j\right)} \]

la cual se distribuye asintóticamente como $N(0,1)$ bajo la veracidad de la hipótesis nula.

$H_0$ se contrasta con la hipótesis alternativa de dos colas $H_1: \beta_j \neq 0$.

A partir de la relación entre intervalos de confianza y pruebas de hipótesis, tenemos que si el intervalo de confianza (IC) de $\beta_j$ está por debajo (por encima) de 0 con nivel $\alpha$, entonces se puede rechazar $H_0$ con el mismo nivel. Las pruebas de significancia están incorporadas en la función summary(object). Sin embargo, debido a discrepancias entre esta función y confint(object, parm, level = 0.95, …), se recomienda precaución al aplicar la regla de decisión anterior para rechazar $H_0$ en términos del IC.

Ejemplo: Hosmer & Lemeshow, S.(1989). Cont.

Código

summary(fit_Bin)

#> 
#> Call:
#> glm(formula = coro ~ edad, family = "binomial", data = dat)
#> 
#> Coefficients:
#>             Estimate Std. Error z value Pr(>|z|)    
#> (Intercept) -5.30945    1.13365  -4.683 2.82e-06 ***
#> edad         0.11092    0.02406   4.610 4.02e-06 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> (Dispersion parameter for binomial family taken to be 1)
#> 
#>     Null deviance: 136.66  on 99  degrees of freedom
#> Residual deviance: 107.35  on 98  degrees of freedom
#> AIC: 111.35
#> 
#> Number of Fisher Scoring iterations: 4

Código

cbind(coefficients(summary(fit_Bin))[,1]-1.96*coefficients(summary(fit_Bin))[,2],coefficients(summary(fit_Bin))[,1]+1.96*coefficients(summary(fit_Bin))[,2])

#>                   [,1]       [,2]
#> (Intercept) -7.5314145 -3.0874922
#> edad         0.0637639  0.1580784

Código

confint.default(fit_Bin)  # based on asymptotic normality

#>                   2.5 %     97.5 %
#> (Intercept) -7.53137370 -3.0875331
#> edad         0.06376477  0.1580775

Código

#based on profiling the least-squares estimation surface
confint(fit_Bin)          # needs MASS to be installed

#>                   2.5 %     97.5 %
#> (Intercept) -7.72587162 -3.2461547
#> edad         0.06693158  0.1620067

Los errores estándar y las pruebas de Wald enfrentan dificultades particulares cuando

los valores ajustados en modelos lineales generalizados binomiales (GLMs) están muy cerca de cero o uno.
Cuando el predictor lineal incluye factores, en la práctica puede ocurrir que haya un nivel del factor en el que los valores ( $y_i$ ) sean todos cero o todos uno. En esta situación, los valores ajustados estimados por el modelo también serán cero o uno para ese nivel del factor.

A pesar de los problemas con las pruebas de Wald, la prueba de razón de verosimilitud (y la prueba score) suele seguir siendo bastante útiles en estas situaciones, incluso cuando los valores ajustados son cero o uno.

Sobre los modelos y ligas en glm La función glm de R ofrece estas opciones para la familia con las siguientes funciones de enlace predeterminadas:

Familia	Funciones ligas:
binomial	“logit”, “probit”, “cloglog”
gaussian	“identity”
Gamma	“inverse”, “identity”,“log”
inverse.gaussian	“1/mu^2”, “inverse”, “identity”, “log”
poisson	“log”, “identity”, “sqrt”

5.2.1.4.2 Residuales de Pearson

Las distancias $Y_i-\mu_i$ se denominan los residuales de respuesta, pero no son adecuados para usarse directamente para evaluar el modelo ya que la varianza de estos puede depender de su media.

En forma similar a los residuales Studentizados, los modificaremos para inducir una varianza constante. Una opción son los residuales de Pearson definidos como \[ r_{i\text{P}}=\frac{Y_i-\widehat{\mu}}{\sqrt{\widehat{\text{Var}}({Y_i})}}. \tag{5.4}\]

Para la regresión logística tenemos \[ r_{i\text{P}}=\frac{y_i-\widehat{\pi}_i}{\sqrt{\widehat{\pi}_i\left(1-\widehat{\pi}_i\right)}}, \qquad \text{donde} \qquad \widehat{\pi_i}= \frac{1}{1+\exp(-\boldsymbol{x}_i^{\top}\boldsymbol{\beta})}. \] En R, puedes obtener distintos tipos de residuos utilizando la función ‘residuals’ y el argumento ‘type’. Por ejemplo: residuals(objeto_glm, type = “pearson”).

Los gráficos de los residuos frente a los valores ajustados $\widehat{\mu}$ y los residuos frente a $x_j$ son las principales herramientas para el análisis diagnóstico. Se recomienda utilizar residuos de desviación estandarizados o residuos cuantílicos en estos gráficos, ya que presentan una varianza aproximadamente constante.

Si el modelo de regresión logística ajustado fuera correcto, esperaríamos ver una banda horizontal con la mayoría de los residuos cayendo dentro de ±3.
Si el número de observaciones en cada categoría, $n_i$ son pequeños, entonces el gráfico podría no tener un patrón tan claro, incluso si el modelo es correcto.

Ejemplo: Hosmer & Lemeshow, S.(1989). Cont.

Código

residuals(fit_Bin, type = "pearson")

#>          1          2          3          4          5          6          7 
#> -0.2132020 -0.2517966 -0.2661559 -0.2813341  3.5544929 -0.2973778 -0.2973778 
#>          8          9         10         11         12         13         14 
#> -0.3322623 -0.3322623 -0.3512103 -0.3712389 -0.3712389 -0.3712389 -0.3712389 
#>         15         16         17         18         19         20         21 
#> -0.3712389  2.6936834 -0.4147877 -0.4147877 -0.4384420 -0.4384420 -0.4634451 
#>         22         23         24         25         26         27         28 
#> -0.4634451  2.1577527 -0.4634451 -0.4634451 -0.4898742 -0.4898742 -0.5178104 
#>         29         30         31         32         33         34         35 
#>  1.9312088 -0.5178104 -0.5473397  1.8270188 -0.5473397 -0.5785531 -0.5785531 
#>         36         37         38         39         40         41         42 
#> -0.6115464  1.6351988 -0.6464213  1.5469788 -0.6832850 -0.6832850 -0.7222509 
#>         43         44         45         46         47         48         49 
#> -0.7222509 -0.7222509  1.3845604 -0.7634389 -0.7634389  1.3098624 -0.8069758 
#>         50         51         52         53         54         55         56 
#> -0.8069758  1.2391945  1.2391945 -0.8529955  1.1723391 -0.9016396  1.1090906 
#>         57         58         59         60         61         62         63 
#> -0.9530577 -0.9530577  1.0492544 -1.0074080  0.9926464  0.9926464 -1.0648578 
#>         64         65         66         67         68         69         70 
#> -1.0648578  0.9390925 -1.1255838  0.8884278 -1.1897729 -1.2576225  0.7951512 
#>         71         72         73         74         75         76         77 
#>  0.7522522  0.7522522  0.7116677 -1.4852821  0.6732728  0.6732728  0.6369493 
#>         78         79         80         81         82         83         84 
#>  0.6369493  0.6369493 -1.6595158 -1.6595158  0.6025854  0.6025854  0.6025854 
#>         85         86         87         88         89         90         91 
#>  0.6025854 -1.7541536  0.5700755  0.5700755  0.5393196  0.5393196 -1.9599277 
#>         92         93         94         95         96         97         98 
#>  0.5102229  0.4826960  0.4566543  0.4566543  0.4320175 -2.4467237  0.4087098 
#>         99        100 
#>  0.3866596  0.3097305

Código

hist(residuals(fit_Bin,type = "pearson"))

Código

#Grafica de residuales contra mu gorro
plot(fit_Bin$fitted.values, residuals(fit_Bin, type = "pearson"))

5.2.1.4.3 Desviaciones (Deviances)

Hay otra forma de pensar en los residuos en este contexto, recurriendo a la idea de desviación (como la de prueba de razón de verosimilitudes).

La desviación es un concepto clave en los modelos lineales generalizados. Intuitivamente, mide la discrepancia del modelo lineal generalizado ajustado con respecto a un modelo “perfecto” para la muestra $\left\{\left(\mathbf{x}_i, Y_i\right)\right\}_{i=1}^n$. Este modelo perfecto, conocido como modelo saturado, es aquel que ajusta exactamente los datos, en el sentido de que las respuestas ajustadas $\left(\widehat{Y}_i\right)$ son iguales a las respuestas observadas $\left(Y_i\right)$. Por ejemplo, en regresión logística, este modelo sería aquel en el que

\[ \widehat{\mathbb{P}}\left[Y=1 \mid X_1=X_{i 1}, \ldots, X_k=X_{i k}\right]=Y_i, \quad i=1, \ldots, n. \]

En este sentido el modelo saturado es una idealización de un modelo perfecto con verosimilitud igual a \[ L(\boldsymbol{\beta})=\prod_{i=1}^n\hat{\pi_i}^{y_i}(1-\hat{\pi}_i)^{1-y_i}, \] donde si $y_i=1$ el modelo estima $\hat{\pi}_i=1$ y si $y_i=1$, el modelo estima $\hat{\pi}_i=0$. Por lo tanto la verosimiliutd del modelo saturado, $L_{\text{sat}}(\boldsymbol{\hat{\beta}})$, es igual a 1.

Formalmente, la desviación se define mediante la diferencia entre los log-verosímiles del modelo ajustado, $\ell(\widehat{\boldsymbol{\beta}})$, y el modelo saturado, $\ell_{\text{sat}}$. Calcular $\ell_{\text{sat}}$ implica sustituir $\mu_i$ por $Y_i$. La desviación se define como: \[ D := 2 \log \left(\frac{L_{\text {sat }}(\widehat{\boldsymbol{\beta}})}{L_{\text {modelo}}(\widehat{\boldsymbol{\beta}})}\right)=2\left(\ell_{\text {sat }}(\widehat{\boldsymbol{\beta}})-\ell_{\text {modelo }}(\widehat{\boldsymbol{\beta}})\right). \] con grados de libertad $\text{df} = \text{df}_{\text{sat}} - \text{df}_{\text{modelo}}$.

Entonces la deviación para el modelo saturado es \[ −2\log(1/1)=0, \] con 0 grados de libertad.

La log-verosimilitud $\ell_{\text {modelo }}(\widehat{\boldsymbol{\beta}})$ es siempre menor que $\ell_{\text {sat }}$ (ya que este coincide perfectamente con la muestra). Como consecuencia, la desviación es siempre mayor o igual a cero, siendo cero únicamente cuando el ajuste del modelo es perfecto.

En el otro lado del espectro podemos considerar el modelo que tiene como único parámetro el correspondiente a la ordenada al origen en la función lineal, $\beta_0$. A este modelo nos referimos como el modelo nulo. Ahora definimos a la desviación nula (Null deviance) como: \[ D_0 := 2 \log \left(\frac{L_{\text {sat }}(\widehat{\boldsymbol{\beta}})}{L_{\text {modelo nulo}}(\widehat{\beta_0})}\right)=2\left(\ell_{\text {sat }}(\widehat{\boldsymbol{\beta}})-\ell_{\text {modelo nulo}}(\widehat{\beta_0})\right). \] y con grados de libertad $\text{df} = \text{df}_{\text{sat}} - \text{df}_{\text{Null}}$.

El modelo saturado es un modelo que asume que cada punto de datos tiene sus propios parámetros (lo que significa que tienes que estimar $n$ parámetros).
El modelo nulo asume lo contrario: considera un solo parámetro para todos los puntos de datos, lo que significa que solo estimas un parámetro.
El modelo propuesto asume que puedes explicar los puntos de datos con $p$ parámetros (incluyendo el intercepto).

En R, se puede obtener usando la función summary() evaluado en el objeto del modelo ajustado.

Verosimilitudes componentes de desviaciones

Utilizando la desviación y la desviación nula, podemos comparar cuánto ha mejorado el modelo al agregar los predictores $X_1, \ldots, X_k$ y cuantificar el porcentaje de desviación explicada. Esto se puede hacer mediante un elemento similar a la estadística $R^2$

5.2.1.4.4 Pseudo $R^2$

Es una generalización del coeficiente de determinación en la regresión lineal: \[ R^2:=1-\frac{D}{D_0}, \qquad \text{que en el modelo lineal es igual a} \quad 1-\frac{\mathrm{SS_{\text{Res}}}}{\mathrm{SS_{\text{T}}}}. \]

El término $R^2$ en los modelos lineales generalizados (Pseudo $R^2$) sigue la misma filosofía que el coeficiente de determinación en la regresión lineal: es una medida proporcional de qué tan bueno es el ajuste del modelo. Si el ajuste es perfecto, $D=0$ y $R^2=1$. Si los predictores no tienen relación con $Y$, entonces $D=D_0$ y $R^2=0$.

Sin embargo, este $R^2$ tiene una interpretación diferente a la de la regresión lineal. En particular:

x- No representa el porcentaje de varianza explicada por el modelo, sino más bien una proporción que indica qué tan cercano es el ajuste al mejor o al peor escenario. - No está relacionado con ningún coeficiente de correlación.

La desviación se obtiene mediante la función summary. Es importante recordar que en R la desviación se denomina “Residual deviance”, mientras que la desviación nula se conoce como “Null deviance”.

Ejemplo: Hosmer & Lemeshow, S.(1989). Cont.

Código

summary(fit_Bin)

#> 
#> Call:
#> glm(formula = coro ~ edad, family = "binomial", data = dat)
#> 
#> Coefficients:
#>             Estimate Std. Error z value Pr(>|z|)    
#> (Intercept) -5.30945    1.13365  -4.683 2.82e-06 ***
#> edad         0.11092    0.02406   4.610 4.02e-06 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> (Dispersion parameter for binomial family taken to be 1)
#> 
#>     Null deviance: 136.66  on 99  degrees of freedom
#> Residual deviance: 107.35  on 98  degrees of freedom
#> AIC: 111.35
#> 
#> Number of Fisher Scoring iterations: 4

Código

1-fit_Bin$deviance/fit_Bin$null.deviance

#> [1] 0.2144684

Entonces el término $R^2$ si muestra una mejora en el modelo ajustado en comparación al modelo nulo.

Como la desviación está expresada como la diferencia de los logaritmos de verosimilitudes, se puede descomponerse en una suma de valores provenientes de cada punto de datos, y se puede escribir como: \[ D=\sum_{i=1}^n d_i, \qquad \text{con }\qquad d_i=2\left(\ell_{\text {sat }}(\widehat{\beta};Y_i)-\ell_{\text {modelo }}(\widehat{\beta};Y_i)\right) \] donde $d_i$ representa la desviación debida al punto de datos $i$.

El residuo de desviación para un punto individual se define como: \[ \operatorname{dev}_i=\operatorname{signo}\left(y_i-\widehat{\mu}_i\right) \sqrt{d_i} \]

Observa que los residuos de desviación utilizan el signo de “valor observado menos predicho”. Así que, al igual que antes, si subestimaste la respuesta, obtendrás un residuo positivo, mientras que si la sobrestimaste, obtendrás un residuo negativo. Pero aquí, la magnitud del residuo se define utilizando esa desviación individual.

Un punto con un residuo grande es uno que puede tener una gran influencia en la calidad de la inferencia.

En R podemos usar summary(objeto_glm)$deviance.resid.

5.2.2 Regresión Multinomial

Podemos extender el caso binomial ( con $r=2$ categorías) al caso $r > 2$, es decir cuando tenemos una variable de respuesta multicategórica.

La Regresión Logística Multinomial modela cómo una variable de respuesta multinomial $Y$ depende de un conjunto de $k$ variables explicativas, $\boldsymbol{x} = (x_1, x_2, \ldots, x_k).$

En este GLM (Modelo Lineal Generalizado) el componente aleatorio asume que la distribución de $Y$ es multinomial($n, \boldsymbol{\pi}$), donde $\boldsymbol{\pi}$ es un vector con las probabilidades de “éxito” para las categorías.

Al igual que en la regresión logística binaria, el componente sistemático consiste en variables explicativas (pueden ser continuas, discretas o ambas) y son lineales en los parámetros. La función de enlace es el logit generalizado, es decir, el enlace logit para cada par de logits no redundantes.

En este caso la regresión logit se extiende a $r(r-1)/2$ logits (razones de probabilidad) que podemos formar, por ejemplo:

\[ \log\left(\frac{\pi_{\text{cat A}}}{\pi_{\text{cat B}}}\right) = \beta_0 + \beta_1 X_1 + \ldots + \beta_k X_k, \]

pero solo $(r-1)$ son no redundantes.

Ilustramos este último punto a continuación.

Con $r=3$ categorías, hay $r \times(r-1)=3 \times 2=6$ pares (A vs B, B vs A, A vs C, C vs A, B vs C, C vs B). Como el logit de A vs B es el negativo del logit de B vs A \[ \log \left(\frac{\pi_A}{\pi_B}\right) = -\log \left(\frac{\pi_B}{\pi_A}\right), \] los pares no ordenados son ${r(r-1)}/{2}=3$:

$L_{A B}=\log \left(\pi_A / \pi_B\right)$
$L_{A C}=\log \left(\pi_A / \pi_C\right)$
$L_{B C}=\log \left(\pi_B / \pi_C\right)$

Por otro lado, estos tres logits no son independientes. Existe una relación exacta que los conecta: \[ L_{B C}=\log \left(\frac{\pi_B}{\pi_C}\right)=\log \left(\frac{\pi_B}{\pi_A} \cdot \frac{\pi_A}{\pi_C}\right)=\log \left(\frac{\pi_B}{\pi_A}\right)+\log \left(\frac{\pi_A}{\pi_C}\right). \] Así, \[ L_{B C}=L_{B A}+L_{A C}=-L_{A B}+L_{A C} \] y el tercer logit ($L_{B C}$) se puede calcular exactamente como una combinación lineal de los otros dos $\left(L_{A C}\right.$ y $\left.L_{A B}\right)$.

Existen diferentes formas de construir un conjunto de $(r-1)$ logits no redundantes, y estas conducirán a diferentes modelos de regresión logística multinomial.

Al analizar una respuesta multicategórica, es importante notar si la respuesta es ordinal (categorías ordenadas) o nominal (categorías no ordenadas). En el modelo logístico binario, esta cuestión no surge. Algunos tipos de modelos son apropiados solo para respuestas ordinales (por ejemplo, modelo de logits acumulativos, modelo de categorías adyacentes).

Si la respuesta es ordinal, no es obligatorio tener en cuenta el orden, pero no hacerlo puede llevar a modelos subóptimos. Usar el orden natural puede:

conducir a un modelo más simple y parsimonioso, y
aumentar la potencia para detectar relaciones con otras variables.

Ejemplo: High School and Beyond

Doscientas observaciones fueron muestreadas aleatoriamente de la encuesta “High School and Beyond”, una encuesta realizada a estudiantes de último año de secundaria por el Centro Nacional de Estadísticas Educativas.

La base tiene 200 observaciones y 11 variables:

id. Student ID.
gender. Student’s gender, with levels female and male.
race. Student’s race, with levels african american, asian, hispanic, and white.
ses. Socio economic status of student’s family, with levels low, middle, and high.
schtyp. Type of school, with levels public and private.
prog. Type of program, with levels general, academic, and vocational.
read. Standardized reading score.
write. Standardized writing score.
math. Standardized math score.
science. Standardized science score.
socst. Standardized social studies score.

Código

library(foreign)
library(nnet)
library(stargazer)


mydata = read.csv("https://raw.githubusercontent.com/leticiaram/data/refs/heads/main/hsb2.csv",head=TRUE)

head(mydata)

#>    id gender  race    ses schtyp       prog read write math science socst
#> 1  70   male white    low public    general   57    52   41      47    57
#> 2 121 female white middle public vocational   68    59   53      63    61
#> 3  86   male white   high public    general   44    33   54      58    31
#> 4 141   male white   high public vocational   63    44   47      53    56
#> 5 172   male white middle public   academic   47    52   57      53    61
#> 6 113   male white middle public   academic   44    52   51      63    61

Código

table(mydata$ses)

#> 
#>   high    low middle 
#>     58     47     95

Código

mydata$ses<-as.factor(mydata$ses)
mydata$gender<-as.factor(mydata$gender)

# Por devfault la primera categoria es la referencia.
# para cambiarla a low hacemos:
mydata$ses2 = relevel(mydata$ses, ref = "low")

multi1 = multinom(ses2 ~ read + science + socst, data=mydata)

#> # weights:  15 (8 variable)
#> initial  value 219.722458 
#> iter  10 value 196.207110
#> final  value 196.180071 
#> converged

Código

multi1

#> Call:
#> multinom(formula = ses2 ~ read + science + socst, data = mydata)
#> 
#> Coefficients:
#>        (Intercept)         read    science      socst
#> high      -6.50201  0.018211643 0.04339229 0.06714246
#> middle    -2.26803 -0.003267788 0.02956575 0.03351177
#> 
#> Residual Deviance: 392.3601 
#> AIC: 408.3601

La función multinom() no reporta los p-valores, pero se pueden calcular.

Código

z <- summary(multi1)$coefficients / summary(multi1)$standard.errors

p <- (1 - pnorm(abs(z), 0, 1)) * 2

p

#>         (Intercept)      read   science       socst
#> high   5.106042e-06 0.5506418 0.1222936 0.009825357
#> middle 4.457373e-02 0.9028890 0.2053895 0.114949253

5.2.3 Regresión Poisson

Con fecuencia se tiene que modelar el número de eventos que ocurren (o se cuentan) durante un cierto periodo de tiempo. Uno de los primeros candidatos que podemos considerar, si

el número de eventos son independientes en intervalos ajenos de tiempo y
la probabilidad de ocurrencia de dos eventos o más se vuelve despreciable en intervalos pequeños es un proceso Poisson pero en el que la tasa de ocurrencia no es constante. Es decir, la tasas de incidencia puede no solo variar en el tiempo sino que puede depender de ciertas variables predictoras. Para modelar este tipo de datos se utiliza la distribución de Poisson, $Y \sim \operatorname{Pois}(\mu)$, donde $\mu=\mu(\boldsymbol{x})$.

En general, el parámetro describe el número promedio de observaciones por unidad de tiempo, pero el proceso puede considerar también otra unidades como unidades de espacio o unidades en exposición.

En este tipo de modelos es muy habitual trabajar con variables categóricas de tipo ordinal. Dichas variables pueden ser incorporadas en el modelo como factores pero en ocasiones, cuando son ordinales, interesa introducirlas como numéricas para captar las posibles tendencias de la respuesta conforme aumenta la relevancia del factor. Recordemos que una variable categórica ordinal se puede obtener siempre a partir de una variable de tipo numérico. Para poder hacer esto introducimos un código numérico artificial asociado con cada categoría. Disponemos de dos alternativas:

Asignar un código continuo: $1,2,3, \ldots$ donde los valores más bajos se asocian con los niveles más bajos de la variable categórica.
Cuando la variable categórica viene dada en términos de un intervalo se puede usar el punto medio de dicho intervalo para establecer el código numérico.

De hecho, en el análisis descriptivo inicial de los datos (sobre todo en los gráficos) esta forma de proceder permite analizar con más detalle la evolución del factor ordinal en la respuesta.

El modelo de regresión de Poisson se utiliza para modelar datos de recuento y tablas de contingencia. Asume el logaritmo de los valores esperados (media) que pueden modelarse en forma lineal mediante algunos parámetros desconocidos. Esto es la función liga es $\log()$ y así el MLG es \[ \log\left(\mathbb{E}\left[Y\right]\right)=\log\left(\mu\right)=\eta_i=\boldsymbol{x}^{\top}\boldsymbol{\beta}=\sum_{j=0}^k \beta_j x_{j}, \] donde $Y\sim\operatorname{Pois}(\mu)$.

La componente sistemáatica del modelo se puede escribir como:

\[ \mu = \exp\{\beta_0\}\exp\{\beta_1\}^{x_{1}}\exp\{\beta_2\}^{x_{2}}\cdots \exp\{\beta_k\}^{x_{k}}, \] de donde es evidente que el impacto de cada variable explicativa en la media, es de tipo multiplicativo. Esto es, incrementando $x_{j}$ en una unidad, modificamos $\mu$ por un factor $\exp\{\beta_j\}$. Si $\beta_j=0$ entonces $\exp\{\beta_j\}=1$ que permanece constante independientemente del valor de $x_{j}$, mientras que si $\beta_j>0$, $\mu$ incrementa con $x_{j}$ y tenemos el caso contrario cuando $\beta_j<0$.

Aunque la función liga $\log()$ es muy popular, otras opciones son $\eta=\mu$ y $\eta=\sqrt{\mu}$.

Cuando todas las covariables son categóricas, los datos pueden resumirse como tablas de contingencia y el modelo se suele denominar como “modelo log-linear”. En caso de que alguna de las covariables son numéricas, el modelo se denomina “regresión Poisson”.

Ejemplo: Maron, M. (2007)

El número de mielero ruidosos (Manorina melanocephala) detectados en varios transectos de 2 hectáreas en parches de bosques de buloke dentro de las llanuras de Wimmera en el oeste de Victoria, Australia.

Región de Wimmera y foto de Minero ruidoso

La investigación tenía como objetivo determinar si la presencia de eucaliptos en los bosques de buloke facilita la invasión del mielero ruidoso (Manorina melanocephala), un competidor agresivo ausente en los bosques de buloke puros.

Se realizó un estudio de aves en remanentes de bosque de buloke donde los eucaliptos eran una especie subdominante con densidades de 0 a 16 por hectárea. El estudio indaga sobre la probabilidad de presencia del mielero ruidoso en los bosques de buloke y en relación a la precencia de eucaliptos.

La base de datos contiene 9 variables entre las cuales hay:

Miners Variable que indica la presencia o ausencia de mieleros ruidosos. Es un vector numérico con niveles 1 (presencia) o 0 (ausencia).
Eucs El número de ecualiptos en cada 2 hectáreas. whether the area was grazed or not; a numeric vector with levels 0 (grazed) or 1 (not grazed)
Bulokes El número de arboles de buloke en cada 2 hectáreas.
Timber El número de piezas de madera caída en el transecto.
Minerab El número de mieleros ruidosos (abundancia) observados en tres encuentas de 20 minutos.

Código

library(GLMsData)
data(nminer)

nm.m1 <- glm(Minerab ~ Eucs, data = nminer, family = poisson)
nm.m2 <- glm(Minerab ~ Bulokes, data = nminer, family = poisson)
nm.m3 <- glm(Minerab ~ Eucs + Bulokes + Timber, data = nminer, family = poisson)
summary(nm.m1)

#> 
#> Call:
#> glm(formula = Minerab ~ Eucs, family = poisson, data = nminer)
#> 
#> Coefficients:
#>             Estimate Std. Error z value Pr(>|z|)    
#> (Intercept) -0.87621    0.28279  -3.098  0.00195 ** 
#> Eucs         0.11398    0.01243   9.169  < 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: 150.545  on 30  degrees of freedom
#> Residual deviance:  63.318  on 29  degrees of freedom
#> AIC: 121.47
#> 
#> Number of Fisher Scoring iterations: 5

Código

summary(nm.m2)

#> 
#> Call:
#> glm(formula = Minerab ~ Bulokes, family = poisson, data = nminer)
#> 
#> Coefficients:
#>              Estimate Std. Error z value Pr(>|z|)    
#> (Intercept)  1.776648   0.215383   8.249  < 2e-16 ***
#> Bulokes     -0.005867   0.001587  -3.696 0.000219 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> (Dispersion parameter for poisson family taken to be 1)
#> 
#>     Null deviance: 150.55  on 30  degrees of freedom
#> Residual deviance: 134.85  on 29  degrees of freedom
#> AIC: 193
#> 
#> Number of Fisher Scoring iterations: 6

Código

summary(nm.m3)

#> 
#> Call:
#> glm(formula = Minerab ~ Eucs + Bulokes + Timber, family = poisson, 
#>     data = nminer)
#> 
#> Coefficients:
#>              Estimate Std. Error z value Pr(>|z|)    
#> (Intercept) -1.531667   0.517195  -2.961  0.00306 ** 
#> Eucs         0.129355   0.016794   7.702 1.34e-14 ***
#> Bulokes      0.002680   0.001616   1.659  0.09709 .  
#> Timber       0.001246   0.006863   0.182  0.85595    
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> (Dispersion parameter for poisson family taken to be 1)
#> 
#>     Null deviance: 150.545  on 30  degrees of freedom
#> Residual deviance:  60.645  on 27  degrees of freedom
#> AIC: 122.8
#> 
#> Number of Fisher Scoring iterations: 6

Calculamos las $R^2$ de las regresiones Poisson

Código

1-nm.m1$deviance/nm.m1$null.deviance

#> [1] 0.5794092

Código

1-nm.m2$deviance/nm.m2$null.deviance

#> [1] 0.1042727

Código

1-nm.m3$deviance/nm.m3$null.deviance

#> [1] 0.5971666

Código

exp(nm.m1$coefficients)

#> (Intercept)        Eucs 
#>   0.4163573   1.1207312

Entonces

0.4163 es el número promedio de mieleros cuando el numero de Eucaliptos es 0
Por cada Eucalipto que se cuenta, el número promedio de mieleros aumenta un factor de 1.1207

Ahora graficamos las desviaciones residuales

Código

plot(nminer$Minerab,summary(nm.m1)$deviance.resid)

Ejemplo: Especies

Consideramos la base de datos species.txt que contiene datos de 90 parcelas de distintos países en las que se midieron la biomasa, el pH del terreno (variable categórica) y el número de especies.

Código

especies <- read.table("https://raw.githubusercontent.com/egarpor/handy/master/datasets/species.txt", sep="\t", head=TRUE)

La variable pH es categórica con valores “high”, “low”, “med”. La transformamos a factor para que en el modelo sea identificada como tal.

Código

str(especies)

#> 'data.frame':    90 obs. of  3 variables:
#>  $ pH     : chr  "high" "high" "high" "high" ...
#>  $ Biomass: num  0.469 1.731 2.09 3.926 4.367 ...
#>  $ Species: int  30 39 44 35 25 29 23 18 19 12 ...

Código

table(especies$pH)

#> 
#> high  low  med 
#>   30   30   30

Código

especies$pH <- as.factor(especies$pH)
str(especies$pH)

#>  Factor w/ 3 levels "high","low","med": 1 1 1 1 1 1 1 1 1 1 ...

Graficamos los datos y coloreamos de acuerdo a la variable pH

Código

plot(Species ~ Biomass, data = especies, col = as.numeric(pH), pch=19)
legend("topright", legend = c("High pH", "Medium pH", "Low pH"),
       col = c(1, 3, 2), lwd = 2)

Ajustamos el modelo Poisson

Código

mod_especies1 <- glm(Species ~ ., data = especies, family = poisson)
summary(mod_especies1)

#> 
#> Call:
#> glm(formula = Species ~ ., family = poisson, data = especies)
#> 
#> Coefficients:
#>             Estimate Std. Error z value Pr(>|z|)    
#> (Intercept)  3.84894    0.05281  72.885  < 2e-16 ***
#> pHlow       -1.13639    0.06720 -16.910  < 2e-16 ***
#> pHmed       -0.44516    0.05486  -8.114 4.88e-16 ***
#> Biomass     -0.12756    0.01014 -12.579  < 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: 452.346  on 89  degrees of freedom
#> Residual deviance:  99.242  on 86  degrees of freedom
#> AIC: 526.43
#> 
#> Number of Fisher Scoring iterations: 4

Interpretamos los coeficientes

Código

exp(mod_especies1$coefficients)

#> (Intercept)       pHlow       pHmed     Biomass 
#>  46.9433686   0.3209744   0.6407222   0.8802418

46.9433 es el promedio de número de especies cuando Biomass = 0 y the pH es ‘high’
Por cada incremento de una unidad en Biomass, el número de especies decrece por un factor de 0.88 (12%)
Si pH decrece a ‘med’ (‘low’), el número de especies decrece por un factor de 0.6407 (0.3209)

Modelo con interacciones

Código

mod_especies2 <- glm(Species ~ Biomass * pH, data = especies, family = poisson)
summary(mod_especies2)

#> 
#> Call:
#> glm(formula = Species ~ Biomass * pH, family = poisson, data = especies)
#> 
#> Coefficients:
#>               Estimate Std. Error z value Pr(>|z|)    
#> (Intercept)    3.76812    0.06153  61.240  < 2e-16 ***
#> Biomass       -0.10713    0.01249  -8.577  < 2e-16 ***
#> pHlow         -0.81557    0.10284  -7.931 2.18e-15 ***
#> pHmed         -0.33146    0.09217  -3.596 0.000323 ***
#> Biomass:pHlow -0.15503    0.04003  -3.873 0.000108 ***
#> Biomass:pHmed -0.03189    0.02308  -1.382 0.166954    
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> (Dispersion parameter for poisson family taken to be 1)
#> 
#>     Null deviance: 452.346  on 89  degrees of freedom
#> Residual deviance:  83.201  on 84  degrees of freedom
#> AIC: 514.39
#> 
#> Number of Fisher Scoring iterations: 4

Código

int<-confint.default(mod_especies2)
int

#>                     2.5 %      97.5 %
#> (Intercept)    3.64752555  3.88872163
#> Biomass       -0.13160981 -0.08264976
#> pHlow         -1.01712752 -0.61401495
#> pHmed         -0.51210872 -0.15081642
#> Biomass:pHlow -0.23348884 -0.07656751
#> Biomass:pHmed -0.07711965  0.01333562

Código

plot(0,xlim=c(min(int),max(int)), ylim=c(0,5),t="n",xlab="",ylab="parámetros",main="Intervalos de Confianza")
grid()
segments(int[,1],0:5,int[,2],0:5,col=rgb(.7,0,.2,1),lwd=6)
abline(v=0,lty=2)

Exponencial de coeficientes

Código

exp(mod_especies2$coefficients)

#>   (Intercept)       Biomass         pHlow         pHmed Biomass:pHlow 
#>    43.2987424     0.8984091     0.4423865     0.7178730     0.8563910 
#> Biomass:pHmed 
#>     0.9686112

Comparamos los modelos

Graficamos los ajuste de los dos modelos

Código

plot(Species ~ Biomass, data = especies, col = as.numeric(pH), pch=19)
legend("topright", legend = c("High pH", "Medium pH", "Low pH"),
       col = c(1, 3, 2), lwd = 2) 

# Sin interacciones
bio <- seq(0, 10, l = 100)
z <- mod_especies1$coefficients[1] + mod_especies1$coefficients[4] * bio
lines(bio, exp(z), col = 1)
lines(bio, exp(mod_especies1$coefficients[2] + z), col = 2)
lines(bio, exp(mod_especies1$coefficients[3] + z), col = 3)

# Con interacciones
bio <- seq(0, 10, l = 100)
z <- mod_especies2$coefficients[1] + mod_especies2$coefficients[2] * bio
lines(bio, exp(z), col = 1, lty = 2)
lines(bio, exp(mod_especies2$coefficients[3] + mod_especies2$coefficients[5] * bio + z),
      col = 2, lty = 2)
lines(bio, exp(mod_especies2$coefficients[4] + mod_especies2$coefficients[6] * bio + z),
      col = 3, lty = 2)

Obtenemos $R^2$

Código

1-mod_especies1$deviance/mod_especies1$null.deviance

#> [1] 0.7806071

Código

1-mod_especies2$deviance/mod_especies2$null.deviance

#> [1] 0.8160674

Residuales

Código

par(mfrow=c(1,2),oma=c(0,0,0,0))
hist(residuals(mod_especies1, type = "pearson"), main="Residuales mod especies1")
hist(residuals(mod_especies2, type = "pearson"), main="Residuales mod especies2")

Código

qqnorm(residuals(mod_especies1, type = "pearson"))
qqnorm(residuals(mod_especies2, type = "pearson"))

5.2.4 Regresión Binomial Negativa

5.2.4.1 Sobredispesión

La sobredispersión es un fenómeno estadístico que ocurre cuando la varianza observada en un conjunto de datos es mayor que la esperada según un modelo estadístico determinado, en particular en el contexto de datos de recuento.

En el contexto de la regresión de Poisson, que se utiliza comúnmente para modelar datos de recuento, la sobredispersión es común en muchos datos reales y plantea un desafío importante. La distribución de Poisson supone que la media y la varianza de los datos son iguales. Sin embargo, cuando hay sobredispersión, este supuesto se viola, lo que lleva a errores estándar subestimados y estadísticas de prueba infladas.

Detectar la sobredispersión es un paso crítico en el proceso de análisis de datos.

No toda sobredispersión es real; en algunos casos, la sobredispersión aparente puede ser identificada y corregida mediante ajustes al modelo. Primero, abordaremos la diferencia entre sobredispersión real y aparente, así como una de las estrategias para mitigar estas últimas.

La sobredispersión puede deberse

a una variabilidad excesiva en las probabilidades de respuesta o los conteos o
cuando se incumplen los supuestos de distribución de los datos, como en el caso de datos agrupados que violan la independencia de observaciones asumida en la verosimilitud.

¿Por qué es un problema la sobredispersión? Porque puede provocar que los errores estándar de las estimaciones se reduzcan artificialmente o sean subestimados, lo que puede llevar a que una variable parezca ser un predictor significativo cuando, en realidad, no lo es.

La cantidad de dispersión en el modelo se determina típicamente mediante el estadístico de bondad de ajuste Ji-cuadrada de Pearson (McCullagh & Nelder, 1989), que es una medida del ajuste general del modelo proporcionada por la salida computacional. El estadístico se basa en la suma de los residuales de Pearson (5.4) al cuadrado dividido entre sus grados de libertad: \[ \cal{X}=\sum_{i=1}^n\frac{(r_{i\text{P}})^2}{n-p}. \tag{5.5}\] Bajo ciertas condiciones $(r_{i\text{P}})^2$ tiene una distribución aproximadamente $\chi^2$. Sin embargo, un criterio de sobredispersión puede hacerse directamente a partir del valor de $\cal{X}$ que se conoce con el nombre de índice de dispersión.

Un modelo puede presentar sobredispersión si el valor de $\cal{X}$ es mayor que 1.0.

Pequeñas cantidades de sobredispersión no suelen ser preocupantes; sin embargo, si el índice de dispersión supera 1.25 en modelos de tamaño moderado, puede ser necesario aplicar una corrección. En modelos con un gran número de observaciones, la sobredispersión puede manifestarse incluso con un índice de dispersión de 1.05.

También debemos estar al tanto de la sobredispersión aparente, que ocurre cuando:

El modelo omite predictores explicativos importantes.
Los datos contienen valores atípicos.
El modelo no incluye un número suficiente de términos de interacción.
Un predictor requiere una transformación a otra escala.
La relación asumida entre la respuesta y la función de enlace, así como los predictores, es incorrecta.

Ejemplo: Fuentes de sobredispersion aparente. Basado en ejemplo en Hilbe J.M (2011).

Simulamos los datos Poisson($\mu(\boldsymbol{x})$), con $\log(\mu)=1+0.5x_1- 0.75x_2 + 0.25x_3$.

Código

library(MASS) 
nobs <- 50000          #número de simulaciones
set.seed(163)
#variables explicativas
x1   <- runif(nobs,min=-3,max=3)
x2   <- runif(nobs)
x3   <- runif(nobs)
#simulación de varible de interés
py   <- rpois(nobs, exp(1 + 0.5*x1 - 0.75*x2 + 0.25*x3))

cnt  <- table(py)
df   <- data.frame(prop.table(table(py)))
df$cumulative <- cumsum(df$Freq)
dfall <- data.frame(cnt, Porciento=df$Freq *100, F=df$cumulative)
dfall

#>    py  Freq Porciento       F
#> 1   0  9966    19.932 0.19932
#> 2   1 10212    20.424 0.40356
#> 3   2  7610    15.220 0.55576
#> 4   3  5520    11.040 0.66616
#> 5   4  4011     8.022 0.74638
#> 6   5  3010     6.020 0.80658
#> 7   6  2380     4.760 0.85418
#> 8   7  1954     3.908 0.89326
#> 9   8  1472     2.944 0.92270
#> 10  9  1135     2.270 0.94540
#> 11 10   832     1.664 0.96204
#> 12 11   622     1.244 0.97448
#> 13 12   439     0.878 0.98326
#> 14 13   298     0.596 0.98922
#> 15 14   231     0.462 0.99384
#> 16 15   131     0.262 0.99646
#> 17 16    73     0.146 0.99792
#> 18 17    49     0.098 0.99890
#> 19 18    28     0.056 0.99946
#> 20 19    12     0.024 0.99970
#> 21 20     7     0.014 0.99984
#> 22 21     2     0.004 0.99988
#> 23 22     3     0.006 0.99994
#> 24 23     3     0.006 1.00000

Ahora ajustamos el modelo real

Código

poy1 <- glm(py ~ x1 + x2 + x3, family = poisson)
summary(poy1)

#> 
#> Call:
#> glm(formula = py ~ x1 + x2 + x3, family = poisson)
#> 
#> Coefficients:
#>              Estimate Std. Error z value Pr(>|z|)    
#> (Intercept)  0.994864   0.006945  143.25   <2e-16 ***
#> x1           0.500739   0.001798  278.51   <2e-16 ***
#> x2          -0.740817   0.008898  -83.26   <2e-16 ***
#> x3           0.254381   0.008776   28.98   <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: 157721  on 49999  degrees of freedom
#> Residual deviance:  54262  on 49996  degrees of freedom
#> AIC: 173004
#> 
#> Number of Fisher Scoring iterations: 5

Código

confint(poy1)

#>                  2.5 %     97.5 %
#> (Intercept)  0.9812413  1.0084647
#> x1           0.4972180  0.5042658
#> x2          -0.7582599 -0.7233816
#> x3           0.2371810  0.2715834

Código

print("índice de dispersión")

#> [1] "índice de dispersión"

Código

sum((py-poy1$fit)^2/poy1$fit)/(poy1$df.residual)

#> [1] 0.9958592

A continuación ajustamos el modelo sin una covariable

Código

poy2 <- glm(py ~ x2 + x3, family=poisson)
summary (poy2)

#> 
#> Call:
#> glm(formula = py ~ x2 + x3, family = poisson)
#> 
#> Coefficients:
#>              Estimate Std. Error z value Pr(>|z|)    
#> (Intercept)  1.356715   0.006528  207.83   <2e-16 ***
#> x2          -0.773092   0.008901  -86.86   <2e-16 ***
#> x3           0.272632   0.008795   31.00   <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: 157721  on 49999  degrees of freedom
#> Residual deviance: 149076  on 49997  degrees of freedom
#> AIC: 267816
#> 
#> Number of Fisher Scoring iterations: 5

Código

confint(poy2)

#>                  2.5 %     97.5 %
#> (Intercept)  1.3439097  1.3694995
#> x2          -0.7905407 -0.7556510
#> x3           0.2553951  0.2898711

Código

print("índice de dispersión")

#> [1] "índice de dispersión"

Código

sum((py-poy2$fit)^2/poy2$fit)/(poy2$df.residual)

#> [1] 3.015921

Ahora introducimos unos valores atípicos en las observaciones originales y ajustamos el modelo verdadero

Código

py_atip <- py
mues<-sample(1:nobs,50)
py_atip[mues]<-py_atip[mues]+round(rnorm(50,30,2))
plot(py,py_atip,pch=19)
grid()

Código

poy3 <- glm(py_atip ~ x1 + x2 + x3, family = poisson)
summary(poy3)

#> 
#> Call:
#> glm(formula = py_atip ~ x1 + x2 + x3, family = poisson)
#> 
#> Coefficients:
#>              Estimate Std. Error z value Pr(>|z|)    
#> (Intercept)  1.012524   0.006899  146.77   <2e-16 ***
#> x1           0.492455   0.001779  276.84   <2e-16 ***
#> x2          -0.731210   0.008852  -82.61   <2e-16 ***
#> x3           0.251748   0.008734   28.82   <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: 162277  on 49999  degrees of freedom
#> Residual deviance:  60625  on 49996  degrees of freedom
#> AIC: 179523
#> 
#> Number of Fisher Scoring iterations: 5

Código

confint(poy3)

#>                  2.5 %     97.5 %
#> (Intercept)  0.9989911  1.0260338
#> x1           0.4889715  0.4959445
#> x2          -0.7485624 -0.7138645
#> x3           0.2346304  0.2688689

Código

#plot(py_atip,poy3$fit)

print("índice de dispersión")

#> [1] "índice de dispersión"

Código

sum((py_atip-poy3$fit)^2/poy3$fit)/(poy3$df.residual)

#> [1] 1.709428

La sobredispersión puede ser motivada también por la transformación de la variable de interés, derivados de los procedimientos o tecnologías para capturar los datos, tales como:

datos inflados en cero
datos truncados
datos censurados

Existen diversos métodos para afrontar la sobredispersión real (ver Hilbe, J. M., 2011) pero una opción es utilizar la distribución Binomial Negativa.

5.2.4.2 Binomial negativa como mezcla gamma de distribuciones Poisson

Un modelo de mezcla es una estrategia flexible para abordar la sobredispersión. Dado un conjunto fijo de variables explicativas observadas y una media $\lambda$, supongamos que la distribución de $Y$ es Poisson$(\lambda)$, pero que $\lambda$ varía debido a covariables no medidas. Definiendo $\mu = \mathbb{E}(\lambda)$, tenemos \[ \begin{gathered} \mathbb{E}(Y) = \mathbb{E}[\mathbb{E}(Y \mid \lambda)] = \mathbb{E}(\lambda) = \mu \\[7pt] \operatorname{Var}(Y) = \mathbb{E}[\operatorname{Var}(Y \mid \lambda)] + \operatorname{Var}[\mathbb{E}(Y \mid \lambda)] = \mathbb{E}(\lambda) + \operatorname{Var}(\lambda) = \mu + \operatorname{Var}(\lambda) > \mu. \end{gathered} \]

Un caso relevante de modelo de mezcla para datos de conteo es el siguiente.

Proposición 5.1 (Binomial Negativa como mezcla de distribuciones) Suponemos que:

Dado $\lambda$, la variable $Y$ sigue una distribución Poisson$(\lambda)$.
$\lambda$ sigue una distribución gamma con fd \[ \begin{aligned} &f(\,y \,; \, k, \mu\,)=\frac{(k / \mu)^k}{\Gamma(k)} y^{k-1} e^{-k y / \mu}, \qquad y \geq 0\\[5pt] &\text { para el que tenemos }\quad \mathbb{E}(y)=\mu \quad \text{y} \quad \operatorname{Var}(y)=\mu^2 / k . \end{aligned} \]

Marginalmente, la mezcla gamma de distribuciones de Poisson da lugar a la distribución binomial negativa para $Y$, cuya f.d.p. es:
\[ p(y \,; \mu, k) = \frac{\Gamma(y+k)}{\Gamma(k) \Gamma(y+1)}\left(\frac{\mu}{\mu+k}\right)^y\left(\frac{k}{\mu+k}\right)^k, \quad y=0,1,2, \ldots. \]

Demostración

\[ \begin{aligned} \mathbb{P}\left(Y=y\right) & =\int_0^{\infty} \mathbb{P}\left(Y=y \mid \lambda\right) f(\lambda; k,\mu) d\lambda \\ & =\int_0^{\infty} \frac{\lambda^y \exp\left(-\lambda\right)}{y!} \frac{(k / \mu)^k}{\Gamma(k)} {\lambda}^{k-1} \exp(-k {\lambda} / \mu) d\lambda\\ & =\frac{(k / \mu)^k}{\Gamma(y+1)\Gamma(k)} \int_0^{\infty} \lambda^y \exp \left(-\lambda\right) \lambda^{k-1} \exp(-k\lambda / \mu) d\lambda \\ & =\frac{(k / \mu)^k}{\Gamma(y+1)\Gamma(k)} \int_0^{\infty} \lambda^{y+k-1} \exp \left[-\left(\lambda+k\lambda / \mu\right)\right] d\lambda \\ & =\frac{(k / \mu)^k}{\Gamma(y+1)\Gamma(k)} \int_0^{\infty} \left[\left(\frac{\mu}{\mu+k}\right)u\right]^{y+k-1} \exp \left[-u\right] \left(\frac{\mu}{\mu+k}\right)du \\ & =\frac{(k / \mu)^k}{\Gamma(y+1)\Gamma(k)} \left(\frac{\mu}{\mu+k}\right)^{y+k}\int_0^{\infty} u^{y+k-1} \exp \left[-u\right] du \\ & =\frac{\Gamma(y+k)}{\Gamma(k) \Gamma(y+1)}\left(\frac{\mu}{\mu+k}\right)^y\left(\frac{k}{\mu+k}\right)^k. \end{aligned} \]

En la familia BN de dos parámetros, definamos $\gamma = 1 / k$. Entonces, la variable $y$ tiene:
\[ \mathbb{E}(Y) = \mu \qquad \text{y} \qquad \operatorname{Var}(Y) = \mu + \gamma \mu^2. \]
El índice $\gamma > 0$ actúa como un parámetro de dispersión. Cuanto mayor sea $\gamma$, mayor será la sobredispersión relativa a la distribución de Poisson. Cuando $\gamma \rightarrow 0$, la varianza de $y$ se aproxima a $\mu$, y la distribución binomial negativa converge a la de Poisson.

La distribución BN abarca un espectro más amplio que la distribución de Poisson. Por ejemplo:

En la distribución de Poisson, la moda es la parte entera de la media y solo es 0 cuando $\mu < 1$.
La BN también es unimodal, pero su moda es 0 cuando $\gamma \geq 1$. En otros casos, la moda corresponde a la parte entera de $\mu(1 - \gamma)$.
A diferencia de la Poisson, la moda puede ser 0 para cualquier valor de $\mu$.

Los modelos lineales generalizados (GLMs) binomial negativa suelen utilizar la función de enlace logarítmica, como en los modelos loglineales de Poisson. Para simplificar, consideramos el parámetro de dispersión $\gamma$ como una constante para todas las $n$ observaciones, tratándolo como un valor desconocido, de manera análoga a la varianza en los modelos normales.

A partir de f.d.p. de $Y$ en términos del parámetro de dispersión $\gamma$, y $\eta_i =\log(\mu_i)=\boldsymbol{x}_i^{\top}\boldsymbol{\beta}$, la logverosimilitud para un GLM binomial negativo con $n$ observaciones independientes es:
\[ \begin{aligned} \ell(\boldsymbol{\beta}, \gamma ; \boldsymbol{y}) = & \sum_{i=1}^n \left[\log \Gamma\left(y_i+\frac{1}{\gamma}\right) - \log \Gamma\left(\frac{1}{\gamma}\right) - \log \Gamma\left(y_i+1\right) \right] \\ & + \sum_{i=1}^n \left[y_i \log \left(\frac{\gamma \mu_i}{1+\gamma \mu_i}\right) - \left(\frac{1}{\gamma}\right) \log \left(1+\gamma \mu_i\right)\right]. \end{aligned} \]

Las ecuaciones de verosimilitud obtenidas al diferenciar $\ell(\boldsymbol{\beta}, \gamma ; \boldsymbol{y})$ con respecto a $\boldsymbol{\beta}$ tienen la forma estándar para un GLM:
\[ \sum_{i=1}^n \frac{\left(y_i-\mu_i\right) x_{i j}}{\operatorname{Var}\left(y_i\right)}\left(\frac{\partial \mu_i}{\partial \eta_i}\right) = \sum_i \frac{\left(y_i-\mu_i\right) x_{i j}}{\mu_i+\gamma \mu_i^2} \left(\frac{\partial \mu_i}{\partial \eta_i}\right) = 0, \quad j=1,2, \ldots, p \]
La logverosimilitud produce una matriz Hessiana con el siguiente término:
\[ \frac{\partial^2 \ell(\boldsymbol{\beta}, \gamma ; \boldsymbol{y})}{\partial \beta_j \partial \gamma} = -\sum_i \frac{\left(y_i-\mu_i\right) x_{i j}}{\left(1+\gamma \mu_i\right)^2} \left(\frac{\partial \mu_i}{\partial \eta_i}\right). \]

Así, se cumple que $E\left(\partial^2 \ell / \partial \beta_j \partial \gamma\right) = 0$ para cada $j$, lo que implica que $\beta$ y $\gamma$ son parámetros ortogonales. En consecuencia, $\widehat{\boldsymbol{\beta}}$ y $\widehat{\gamma}$ son asintóticamente independientes, y el error estándar en muestras grandes ($SE$) de $\widehat{\beta}_j$ es el mismo, independientemente de si $\gamma$ es conocido o estimado.

El algoritmo iterativo de mínimos cuadrados ponderados para Fisher scoring es aplicable para ajustar modelos por máxima verosimilitud. La matriz de covarianza estimada de $\widehat{\boldsymbol{\beta}}$ es:
\[ \widehat{\operatorname{Var}}(\widehat{\boldsymbol{\beta}}) = \left(X^{\top} \widehat{W} X\right)^{-1}, \]
donde $W$, para la función de enlace logarítmica, es la matriz diagonal con elementos:
\[ w_i = \left(\frac{\partial \mu_i}{\partial \eta_i}\right)^2 \frac{1}{\operatorname{Var}(Y_i)} = \frac{\mu_i}{1+\gamma \mu_i}. \]

La desviación para un MLG binomial negativo está dada por:
\[ D(\boldsymbol{y}, \widehat{\boldsymbol{\mu}}) = 2 \sum_i \left[Y_i \log \left(\frac{Y_i}{\widehat{\mu}_i}\right) - \left(Y_i + \frac{1}{\widehat{\gamma}}\right) \log \left(\frac{1+\widehat{\gamma} Y_i}{1+\widehat{\gamma} \widehat{\mu}_i}\right)\right], \]
que se aproxima a la desviación del MLG de Poisson cuando $\widehat{\gamma}$ es cercano a 0.

Una forma en la que podemos comparar los modelos lineales generalizados (MLGs) de Poisson y BN con las mismas variables explicativas para determinar si el modelo binomial negativo proporciona un mejor ajuste puede realizarse utilizando los valores de AIC.

Para una prueba de significancia formal, podemos evaluar la hipótesis nula $H_0: \gamma = 0$, dado que la distribución de Poisson es el caso límite de la binomial negativa cuando $\gamma \downarrow 0$.

Dado que $\gamma$ es positivo, el valor $\gamma = 0$ se encuentra en el límite del espacio paramétrico. Como resultado, el estadístico de razón de verosimilitudes no sigue una distribución nula asintótica Ji-cuadrada estándar. En su lugar, sigue una mezcla equitativa entre una distribución puntual en 0 (que ocurre cuando $\widehat{\gamma} = 0$) y una distribución ji-cuadrada con grados de libertad $\operatorname{df} = 1$.

Por lo tanto, el p-valor de la prueba es la mitad del que se obtendría al tratar el estadístico como una Ji-cuadrada con $\operatorname{df} = 1$ (Ver Self y Liang, 1987).

En R

El modelo BN es un modelo lineal generalizado cuando el parámetro de sobredispersión $k$ es conocido. En aplicaciones, usualmente no lo conocemos y debe ser estimado junto con los otros parámetros del modelo.

glm(., family = negative.binomial(k)) requiere que se tenga un valor de $k$ que se pueda proporcionar. Sin emabargo glm.nb(), del paquete MASS, ajusta el modelo BN tradicional donde $k$ es estimado. La función glm.nb() realiza la estimación por máxima verosimilitud tanto para el parámetro $k$ como para los coeficientes del modelo lineal generalizado (MLG), $\beta_j$, de manera simultánea (usa ?glm.nb).

La estimación de $k$ introduce una capa adicional de incertidumbre en un MLG binomial negativo. Sin embargo, como hemos visto arriba, el estimador de máxima verosimilitud $\widehat{k}$ de $k$ tiende a ser no correlacionado con los coeficientes estimados $\widehat{\beta}_j$, según las aproximaciones asintóticas habituales. Por lo tanto, el ajuste del GLM tiende a ser relativamente estable con respecto a la estimación de $k$.

Los MLGs BN producen errores estándar más grandes que sus equivalentes en Poisson, dependiendo del tamaño de $k=1 / \gamma$. Por otro lado, las estimaciones de los coeficientes $\widehat{\beta}_j$ en un GLM BN pueden ser similares a las obtenidas en un MLG de Poisson correspondiente.

A diferencia de glm(), donde la función de enlace predeterminada para cada familia es el enlace canónico, la función de enlace predeterminada en glm.nb() es el enlace logarítmico. De hecho, el enlace logarítmico es casi siempre utilizado en MLGs BN para garantizar que $\mu > 0$ para cualquier valor del predictor lineal. La función glm.nb() también permite el uso de las funciones de enlace “sqrt” y “identity”.

Para MLGs BN, se recomienda mucho el uso de residuos cuantílicos para su diagnóstico (Dunn, y Smyth, 1996).

Ejemplo: Dunn, P. K., & Smyth, G. K. (2018).

Consideramos los datos de pock muestran sobredispersión (Ejemplo 10.9 en Dunn & Smyth, 2018):

Se trata de un experimento donde se evaluó la actividad viral mediante el conteo de pocks (lesiones o placas) en una membrana, sometiendo el virus a una serie de diluciones.

El conjunto de datos pock es un data frame con 48 observaciones y 2 variables:

Count (variable numérica): Es el número de pocks (lesiones) contados en la membrana. Es la variable de respuesta.
Dilution (variable numérica): Es el factor de dilución del medio viral. Es la variable explicativa.

Dilution	Counts
1	116	151	171	194	196	198	208	259
2	71	74	79	93	94	115	121	123	135	142
4	27	33	34	44	49	51	52	59	67	92
8	8	10	15	22	26	27	30	41	44	48
16	5	6	7	7	8	9	9	9	11	20

Utilizamos el logaritmo en base 2 de la dilución como una covariable, dado que los niveles de dilución aumentan en potencias de 2, lo que sugiere que este aspecto fue considerado en el diseño del experimento.

El gráfico de los datos muestra una relación clara entre las variables (panel izquierdo) y evidencia que la varianza aumenta conforme lo hace la media (panel derecho).

Código

pockc <- c(116, 151, 171, 194, 196, 198, 208, 259, 71, 74, 79, 93, 94, 115, 121, 123, 135, 142, 27, 33, 34, 44, 49, 51, 52, 59, 67, 92, 8, 10, 15, 22, 26, 27, 30, 41, 44, 48, 5, 6, 7, 7, 8, 9, 9, 9, 11, 20 )
dilution <- c(rep(1,8),rep(c(2,4,8,16),each=10))
pock<-data.frame(Count=pockc,Dilution=dilution)

head(pock)

#>   Count Dilution
#> 1   116        1
#> 2   151        1
#> 3   171        1
#> 4   194        1
#> 5   196        1
#> 6   198        1

Código

plot(Count ~ jitter(log2(Dilution)), data=pock, las=1, xlab="Log (base 2) of dilution", ylab="Pock mark count")

Código

mn <- with(pock, tapply(Count, log2(Dilution), mean) ) # Media de grupos
vr <- with(pock, tapply(Count, log2(Dilution), var) ) # Varianzas de grupos
plot( log(vr) ~ log(mn), las=1, xlab="Group mean", ylab="Group variance")

Ajustamos un MLG BN, estimando $k$ mediante la función glm.nb() (esta función usa theta para denotar $k$).

Código

library(MASS) 

m.nb <- glm.nb(Count ~ log2(Dilution), data=pock )
m.nb$theta  # Este es el valor de k (llamado theta en MASS)

#> [1] 9.892894

El objeto de salida m.nb incluye información sobre la estimación de $k$.

La salida del modelo glm.nb() se convierte al formato de salida de glm() utilizando glm.convert().

Código

m.nb <- glm.convert(m.nb)
printCoefmat(coef(summary(m.nb, dispersion=1)))

#>                Estimate Std. Error z value  Pr(>|z|)    
#> (Intercept)     5.33284    0.08786  60.697 < 2.2e-16 ***
#> log2(Dilution) -0.72460    0.03886 -18.646 < 2.2e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Dado que $k \approx 10$, el modelo BN utiliza la función de varianza $\operatorname{Var}(\mu) \approx \mu+\mu^2 / 10$.

El coeficiente de variación de la distribución de mezcla ($\gamma = 1 / k$) se estima en aproximadamente 10%, lo cual representa un nivel razonable de variación entre réplicas. La comparación entre los modelos de Poisson y BN muestra que las estimaciones de los parámetros son bastante cercanas, pero los errores estándar son significativamente diferentes.

Código

m1 <- glm(Count ~ log2(Dilution), data=pock,family=poisson)
printCoefmat( coef( summary(m1)) )

#>                 Estimate Std. Error z value  Pr(>|z|)    
#> (Intercept)     5.267932   0.022552 233.596 < 2.2e-16 ***
#> log2(Dilution) -0.680944   0.015443 -44.093 < 2.2e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Comparamos los aic

Código

m1$aic

#> [1] 562.4242

Código

m.nb$aic

#> [1] 410.0057

Además obtenemos unas figuras para el dianostico

QQplot y Distancia de Cook para el modelo Poisson

Código

plot(m1,2)

Código

plot(m1,4)

QQplot y Distancia de Cook para el modelo BN

Código

plot(m.nb,2)

Código

plot(m.nb,4)

Bibliografía:

Dunn, P. K., & Smyth, G. K. (2018). Generalized linear models with examples in R (Vol. 53, p. 16). New York: Springer. Springer link
McCullagh, P. (2019). Generalized linear models. Routledge.

Referencias:

Dunn, P.K., Smyth, G.K. (1996). Randomized quantile residuals. Journal of Computational and Graphical Statistics 5(3), 236–244.
Hilbe, J. M. (2011). Negative binomial regression. Cambridge University Press.
McCullagh, P., & Nelder, J. A. (1989). Generalized linear models (2nd ed.). London: Chapman & Hall
Maron, M. (2007). Threshold effect of eucalypt density on an aggressive avian competitor. Biological Conservation 136, 100–107.
Self, S. G., and Liang, K.-Y. (1987). Asymptotic properties of maximum likelihood estimators and likelihood ratio tests under nonstandard conditions. J. Am. Stat. Assoc. 82: 605–610.

5.1 Estructura de los modelos lineales generalizados

5.1.1 Componente aleatoria

5.1.1.1 Familia Exponencial

5.1.1.2 Variable respuesta

5.1.2 Componente Sistemática

5.1.3 Función liga

5.2 Casos específicos

5.2.1 Regresión con variable binomial

5.2.1.1 Regresión logit

5.2.1.2 Regresión probit

5.2.1.3 Regresión cloglog

5.2.1.4 Más sobre Regresión logit

5.2.1.4.1 Prueba de hipótesis para los coeficientes

5.2.1.4.2 Residuales de Pearson

5.2.1.4.3 Desviaciones (Deviances)

5.2.1.4.4 Pseudo \(R^2\)

5.2.2 Regresión Multinomial

5.2.3 Regresión Poisson

5.2.4 Regresión Binomial Negativa

5.2.4.1 Sobredispesión

5.2.4.2 Binomial negativa como mezcla gamma de distribuciones Poisson