Un Modelo Lineal Generalizado (GLM) asume que la variable respuesta \(Y_i\) pertenece a la familia exponencial de distribuciones.
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:
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.
Un GLM se compone de:
Componente aleatoria:
\(Y_i \sim\) familia
exponencial
Componente sistemática:
\[
\eta_i = X_i \beta
\]
Función de enlace:
\[
g(\mu_i) = \eta_i
\]
| 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 \]
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:
El método usa la recta tangente en \(x_t\) para aproximar la raíz.
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.
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:
En Modelos Lineales Generalizados, el método se usa para maximizar la log-verosimilitud.
En la práctica, se utiliza una variante llamada:
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.
# 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
# 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
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:
Se calcula:
\[ \hat{\phi} = \frac{\sum r_i^2}{n - p} \]
donde:
También puede usarse un test Chi-cuadrado:
\[ X^2 = \sum r_i^2 \]
y comparar con una \(\chi^2_{n-p}\).
# 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
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
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:
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.
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:
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.
Un GLM Gamma tiene:
\[ Y_i \sim \text{Gamma}(\mu_i, \phi) \]
\[ \eta_i = X_i \beta \]
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) \]
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).
Los parámetros se estiman por Máxima Verosimilitud usando:
glm())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
# Residuos deviance
plot(fitted(modelo_gamma),
residuals(modelo_gamma, type="deviance"),
main="Residuos vs Ajustados")
abline(h=0, col="red")
✔ Variable continua
✔ Valores estrictamente positivos
✔ Distribución asimétrica a la derecha
✔ Varianza creciente con la media
| 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.