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 Modelos Lineales Generalizados, 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.
6.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.
6.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.
6.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)
exponential
binomial negativa (con \(n\) conocida)
gamma
geometrica
chi-squared
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 6.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
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 miembreo 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.
6.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.
6.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\) ó \(x_3=x_2^2\), para modelizar interacciones entre variables o efectos curvilíneos de \(x_2\).
6.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{6.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 \(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)=\frac{\exp\{\boldsymbol{x}_i^{\top}\boldsymbol{\beta}\}}{1+\exp\{\boldsymbol{x}_i^{\top}\boldsymbol{\beta}\}}.
\]
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\) dicotomica. 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
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 cloglogcloglog_vals<-1-exp(-exp(x))# Creamos data framedata<-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 ="Funciones liga: Logit, Probit y CLogLog", x ="x", y ="Probability")+scale_color_manual(values =c("Logit"="red", "Probit"="blue", "CLogLog"="green"))+theme_minimal()
Figura 6.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.
6.2 Casos específicos
6.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()\) 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.
6.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\).
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.
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 - p_i) \\ \sum ( y_i - p_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)}).
\] En esta expresión, tenemos \[
\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)}).
\]
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” \[
\tilde{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 \tilde{y}.
\end{align}
\]
Esto es, \(\boldsymbol{\beta}^{(k+1)} = (X^TWX)^{-1}X^TW \tilde{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
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 verosimilitudy<-coron<-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 coronariatolm<-1e-6# tolerancia (norma minima de delta)iterm<-100# numero maximo de iteracionestolera<-1# inicializar toleraitera<-0# inicializar iterahisto<-b# inicializar historial de iteracioneswhile((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+deltatolera<-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
#> [,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 originalxx<-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} 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}})\).
#>
#> 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
6.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
6.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:
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{6.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 Ecuación 6.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{6.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 Ecuación 6.3, podemos obtener el intervalo de confianza (IC) del \((1-\alpha)\times 100\%\) para el coeficiente \(\beta_j\), \(j=0, \ldots, k\):
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\}.
\]
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.
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”
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{6.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.
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
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\) en (5.20). Si se utiliza la función de enlace canónica, esto equivale a establecer (_i = g(Y_i)) (ver (5.12)). La desviación se define entonces 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}}\).
La log-verosimilitud \(\ell_{\text {modelo }}(\hat{\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 la estadística \(R^2\), que 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 termino \(R^2\) en los modelos lineales generalizados 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”.
#>
#> 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ós está expresada como la diferencia de los logaritmos de verosimilitude, 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.
6.2.2 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
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.
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)
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 interaccionesbio<-seq(0, 10, l =100)z<-mod_especies1$coefficients[1]+mod_especies1$coefficients[4]*biolines(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 interaccionesbio<-seq(0, 10, l =100)z<-mod_especies2$coefficients[1]+mod_especies2$coefficients[2]*biolines(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)
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 (Ecuación 6.4) al cuadrado dividido entre sus grados de libertad: \[
\cal{X}=\sum_{i=1}^n\frac{(r_{i\text{P}})^2}{n-p}.
\tag{6.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 simulacionesset.seed(163)#variables explicativasx1<-runif(nobs,min=-3,max=3)x2<-runif(nobs)x3<-runif(nobs)#simulación de varible de interéspy<-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
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.
6.2.3.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 6.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.
\]
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, \(\hat{\boldsymbol{\beta}}\) y \(\hat{\gamma}\) son asintóticamente independientes, y el error estándar en muestras grandes (\(SE\)) de \(\hat{\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 \(\hat{\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 (ver ?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 \(\hat{k}\) de \(k\) tiende a ser no correlacionado con los coeficientes estimados \(\hat{\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):
Dilution
Pock
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
Los datos de En un experimento se colectaron para evaluar la actividad viral, se contabilizaron las marcas de pock en la membrana a diversas diluciones del medio viral
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).
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.
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.
