Familia Exponencial en los Modelos Lineales Generalizados (GLM)

Un Modelo Lineal Generalizado (GLM) asume que la variable respuesta \(Y_i\) pertenece a la familia exponencial de distribuciones.

1. Forma general

Una distribución pertenece a la familia exponencial si su función de densidad (o masa) puede escribirse como:

\[ f(y_i; \theta_i, \phi) = \exp\left( \frac{y_i \theta_i - b(\theta_i)}{\phi} + c(y_i, \phi) \right) \]

donde:

  • \(\theta_i\) = parámetro canónico
  • \(\phi\) = parámetro de dispersión
  • \(b(\theta_i)\) = función cumulante
  • \(c(y_i, \phi)\) = función conocida

2. Media y Varianza

Para distribuciones de la familia exponencial:

\[ \mathbb{E}(Y_i) = \mu_i = b'(\theta_i) \]

\[ \text{Var}(Y_i) = \phi \, b''(\theta_i) \]

También se expresa como:

\[ \text{Var}(Y_i) = \phi \, V(\mu_i) \]

donde \(V(\mu_i)\) es la función de varianza.


3. Estructura de un GLM

Un GLM se compone de:

  1. Componente aleatoria:
    \(Y_i \sim\) familia exponencial

  2. Componente sistemática:
    \[ \eta_i = X_i \beta \]

  3. Función de enlace:
    \[ g(\mu_i) = \eta_i \]


4. Ejemplos importantes

Distribución Media Varianza Enlace Canónico
Normal \(\mu\) \(\phi\) Identidad
Binomial \(n p\) \(n p(1-p)\) Logit
Poisson \(\mu\) \(\mu\) Log

El método de Newton-Raphson es un procedimiento iterativo para encontrar raíces de una función \(f(x)\), es decir, resolver:

\[ f(x) = 0 \]


1. Idea básica

A partir de una aproximación inicial \(x_0\), se construye una sucesión:

\[ x_{t+1} = x_t - \frac{f(x_t)}{f'(x_t)} \]

donde:

  • \(f'(x_t)\) es la derivada de \(f(x)\)
  • \(t\) indica la iteración

El método usa la recta tangente en \(x_t\) para aproximar la raíz.


2. Interpretación geométrica

La recta tangente en el punto \((x_t, f(x_t))\) es:

\[ y = f(x_t) + f'(x_t)(x - x_t) \]

Al hacer \(y = 0\) y despejar \(x\), se obtiene la fórmula iterativa.


3. Extensión multivariada

Para un vector de parámetros \(\boldsymbol{\beta}\):

\[ \boldsymbol{\beta}^{(t+1)} = \boldsymbol{\beta}^{(t)} - H^{-1}(\boldsymbol{\beta}^{(t)}) \nabla \ell(\boldsymbol{\beta}^{(t)}) \]

donde:

  • \(\nabla \ell(\boldsymbol{\beta})\) es el gradiente
  • \(H(\boldsymbol{\beta})\) es la matriz Hessiana
  • \(\ell(\boldsymbol{\beta})\) es la log-verosimilitud

4. Newton-Raphson en GLM

En Modelos Lineales Generalizados, el método se usa para maximizar la log-verosimilitud.

En la práctica, se utiliza una variante llamada:

  • Fisher Scoring

que reemplaza la Hessiana por su valor esperado:

\[ \boldsymbol{\beta}^{(t+1)} = \boldsymbol{\beta}^{(t)} + I^{-1}(\boldsymbol{\beta}^{(t)}) \nabla \ell(\boldsymbol{\beta}^{(t)}) \]

donde \(I(\boldsymbol{\beta})\) es la matriz de información de Fisher.


5. Ejemplo en R

# Ejemplo simple: encontrar raíz de f(x) = x^2 - 2

f <- function(x) x^2 - 2
df <- function(x) 2*x

x <- 1  # valor inicial

for(i in 1:10){
  x <- x - f(x)/df(x)
}

x  # aproximación de sqrt(2)
## [1] 1.414214

6. Condiciones de convergencia

  • Requiere que \(f'(x_t) \neq 0\)
  • Buena aproximación inicial
  • Convergencia cuadrática cerca de la solución
  • Puede divergir si el punto inicial está lejos
# Regresión Poisson con datos por defecto de R (warpbreaks)

# Cargar datos
data(warpbreaks)

# Ajustar modelo Poisson
modelo_poisson <- glm(breaks ~ wool + tension,
                      data = warpbreaks,
                      family = poisson(link = "log"))

# Resumen del modelo
summary(modelo_poisson)
## 
## Call:
## glm(formula = breaks ~ wool + tension, family = poisson(link = "log"), 
##     data = warpbreaks)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  3.69196    0.04541  81.302  < 2e-16 ***
## woolB       -0.20599    0.05157  -3.994 6.49e-05 ***
## tensionM    -0.32132    0.06027  -5.332 9.73e-08 ***
## tensionH    -0.51849    0.06396  -8.107 5.21e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 297.37  on 53  degrees of freedom
## Residual deviance: 210.39  on 50  degrees of freedom
## AIC: 493.06
## 
## Number of Fisher Scoring iterations: 4
# Predicciones en escala original
predict(modelo_poisson, type = "response")
##        1        2        3        4        5        6        7        8 
## 40.12354 40.12354 40.12354 40.12354 40.12354 40.12354 40.12354 40.12354 
##        9       10       11       12       13       14       15       16 
## 40.12354 29.09722 29.09722 29.09722 29.09722 29.09722 29.09722 29.09722 
##       17       18       19       20       21       22       23       24 
## 29.09722 29.09722 23.89035 23.89035 23.89035 23.89035 23.89035 23.89035 
##       25       26       27       28       29       30       31       32 
## 23.89035 23.89035 23.89035 32.65424 32.65424 32.65424 32.65424 32.65424 
##       33       34       35       36       37       38       39       40 
## 32.65424 32.65424 32.65424 32.65424 23.68056 23.68056 23.68056 23.68056 
##       41       42       43       44       45       46       47       48 
## 23.68056 23.68056 23.68056 23.68056 23.68056 19.44298 19.44298 19.44298 
##       49       50       51       52       53       54 
## 19.44298 19.44298 19.44298 19.44298 19.44298 19.44298

1. Teoría de Sobredispersión

En un modelo Poisson se asume que:

\[ \mathbb{E}(Y_i) = \mu_i \]

\[ \text{Var}(Y_i) = \mu_i \]

Es decir, la media es igual a la varianza.

Existe sobredispersión cuando:

\[ \text{Var}(Y_i) > \mu_i \]

Esto provoca:


2. Test basado en estadístico de Pearson

Se calcula:

\[ \hat{\phi} = \frac{\sum r_i^2}{n - p} \]

donde:

  • \(r_i\) = residuos de Pearson
  • \(n\) = número de observaciones
  • \(p\) = número de parámetros

Regla práctica

  • Si \(\hat{\phi} \approx 1\) → No hay sobredispersión
  • Si \(\hat{\phi} > 1.5\) → Posible sobredispersión
  • Si \(\hat{\phi} >> 1\) → Fuerte sobredispersión

También puede usarse un test Chi-cuadrado:

\[ X^2 = \sum r_i^2 \]

y comparar con una \(\chi^2_{n-p}\).


3. Ejemplo en R (datos warpbreaks)

# Datos por defecto
data(warpbreaks)

# Modelo Poisson
modelo_poisson <- glm(breaks ~ wool + tension,
                      family = poisson(link = "log"),
                      data = warpbreaks)

summary(modelo_poisson)
## 
## Call:
## glm(formula = breaks ~ wool + tension, family = poisson(link = "log"), 
##     data = warpbreaks)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  3.69196    0.04541  81.302  < 2e-16 ***
## woolB       -0.20599    0.05157  -3.994 6.49e-05 ***
## tensionM    -0.32132    0.06027  -5.332 9.73e-08 ***
## tensionH    -0.51849    0.06396  -8.107 5.21e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 297.37  on 53  degrees of freedom
## Residual deviance: 210.39  on 50  degrees of freedom
## AIC: 493.06
## 
## Number of Fisher Scoring iterations: 4
# ---- Test de Sobredispersión ----

# Estadístico de Pearson
pearson_chi2 <- sum(residuals(modelo_poisson, type = "pearson")^2)

# Grados de libertad
gl <- modelo_poisson$df.residual

# Parámetro de dispersión
phi_hat <- pearson_chi2 / gl

phi_hat
## [1] 4.261522
# p-valor del test Chi-cuadrado
p_value <- 1 - pchisq(pearson_chi2, df = gl)

p_value
## [1] 0

4. ¿Qué hacer si hay sobredispersión?

Opción 1: Quasi-Poisson

modelo_quasi <- glm(breaks ~ wool + tension,
                    family = quasipoisson(link = "log"),
                    data = warpbreaks)

summary(modelo_quasi)
## 
## Call:
## glm(formula = breaks ~ wool + tension, family = quasipoisson(link = "log"), 
##     data = warpbreaks)
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  3.69196    0.09374  39.384  < 2e-16 ***
## woolB       -0.20599    0.10646  -1.935 0.058673 .  
## tensionM    -0.32132    0.12441  -2.583 0.012775 *  
## tensionH    -0.51849    0.13203  -3.927 0.000264 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for quasipoisson family taken to be 4.261537)
## 
##     Null deviance: 297.37  on 53  degrees of freedom
## Residual deviance: 210.39  on 50  degrees of freedom
## AIC: NA
## 
## Number of Fisher Scoring iterations: 4

Opción 2: Binomial Negativa

library(MASS)

modelo_nb <- glm.nb(breaks ~ wool + tension,
                    data = warpbreaks)

summary(modelo_nb)
## 
## Call:
## glm.nb(formula = breaks ~ wool + tension, data = warpbreaks, 
##     init.theta = 9.944385436, link = log)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   3.6734     0.0979  37.520  < 2e-16 ***
## woolB        -0.1862     0.1010  -1.844   0.0651 .  
## tensionM     -0.2992     0.1217  -2.458   0.0140 *  
## tensionH     -0.5114     0.1237  -4.133 3.58e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Negative Binomial(9.9444) family taken to be 1)
## 
##     Null deviance: 75.464  on 53  degrees of freedom
## Residual deviance: 53.723  on 50  degrees of freedom
## AIC: 408.76
## 
## Number of Fisher Scoring iterations: 1
## 
## 
##               Theta:  9.94 
##           Std. Err.:  2.56 
## 
##  2 x log-likelihood:  -398.764

Si \(\hat{\phi} > 1\), el modelo Poisson viola el supuesto de equidispersión y se recomienda usar:

  • Quasi-Poisson (corrige errores estándar)
  • Binomial Negativa (modela explícitamente la sobredispersión)

Regresión Gamma

La regresión Gamma es un caso particular de los Modelos Lineales Generalizados (GLM) que se utiliza cuando:

Ejemplos típicos: costos médicos, tiempos de espera, ingresos, montos de seguros.


1. Distribución Gamma

Una variable aleatoria \(Y\) sigue una distribución Gamma si:

\[ f(y; \alpha, \lambda) = \frac{\lambda^\alpha}{\Gamma(\alpha)} y^{\alpha - 1} e^{-\lambda y}, \quad y > 0 \]

donde:


2. Gamma como Familia Exponencial

En forma de familia exponencial:

\[ f(y; \theta, \phi) = \exp \left( \frac{y\theta - b(\theta)}{\phi} + c(y, \phi) \right) \]

Para la distribución Gamma:

La función de varianza es:

\[ V(\mu_i) = \mu_i^2 \]

Esto implica que la varianza crece cuadráticamente con la media.


3. Estructura del Modelo Gamma

Un GLM Gamma tiene:

Componente aleatoria:

\[ Y_i \sim \text{Gamma}(\mu_i, \phi) \]

Componente sistemática:

\[ \eta_i = X_i \beta \]

Función de enlace:

La más común es el enlace logarítmico:

\[ g(\mu_i) = \log(\mu_i) \]

Entonces:

\[ \log(\mu_i) = X_i \beta \]

y por lo tanto:

\[ \mu_i = \exp(X_i \beta) \]


4. Interpretación de coeficientes (enlace log)

Si:

\[ \log(\mu_i) = \beta_0 + \beta_1 x_i \]

entonces:

\[ \mu_i = e^{\beta_0} e^{\beta_1 x_i} \]

Un incremento de 1 unidad en \(x\) multiplica la media por:

\[ e^{\beta_1} \]

Interpretación: efecto multiplicativo (similar a Poisson).


5. Estimación

Los parámetros se estiman por Máxima Verosimilitud usando:


6. Ejemplo en R

Usaremos datos simulados positivos.

set.seed(123)

# Simulación
n <- 100
x <- runif(n, 0, 5)
mu <- exp(1 + 0.5*x)

y <- rgamma(n, shape = 2, scale = mu/2)

datos <- data.frame(y, x)

# Modelo Gamma con enlace log
modelo_gamma <- glm(y ~ x,
                    family = Gamma(link = "log"),
                    data = datos)

summary(modelo_gamma)
## 
## Call:
## glm(formula = y ~ x, family = Gamma(link = "log"), data = datos)
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.96549    0.12644   7.636  1.5e-11 ***
## x            0.48975    0.04409  11.108  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Gamma family taken to be 0.3907621)
## 
##     Null deviance: 89.111  on 99  degrees of freedom
## Residual deviance: 41.455  on 98  degrees of freedom
## AIC: 604.42
## 
## Number of Fisher Scoring iterations: 5

7. Diagnóstico básico

# Residuos deviance
plot(fitted(modelo_gamma),
     residuals(modelo_gamma, type="deviance"),
     main="Residuos vs Ajustados")

abline(h=0, col="red")


8. ¿Cuándo usar Regresión Gamma?

✔ Variable continua
✔ Valores estrictamente positivos
✔ Distribución asimétrica a la derecha
✔ Varianza creciente con la media


9. Comparación con otros modelos

Modelo Tipo de variable Varianza
Normal Continua Constante
Poisson Conteo μ
Gamma Continua positiva μ²

La regresión Gamma es ideal para modelar datos positivos y asimétricos donde la varianza aumenta proporcionalmente al cuadrado de la media, y su interpretación con enlace log es multiplicativa.