Anteriormente se introdujo el modelo de regresión lineal, que asume que la varianza es constante (homocedasticidad) y, a menudo, que los errores siguen una distribución normal. Sin embargo, existen muchos tipos de datos en los que la variabilidad no es constante, por lo que el método de mínimos cuadrados no es adecuado. En estos casos es preferible usar la estimación de máxima verosimilitud para ajustar el modelo.
El componente aleatorio de los modelos de regresión en capítulos anteriores asume varianza constante y normalidad.
Existen tres situaciones comunes en las que la variación no es constante y, por ello, los modelos lineales no son adecuados:
En todas estas situaciones, la relación entre \(y\) y las variables explicativas suele ser no lineal, ya que \(y\) está acotado.
Anteriormente se desarrolló el principio de los mínimos cuadrados como criterio para estimar los parámetros en el predictor lineal de los modelos de regresión lineal. Este método es adecuado cuando los datos de respuesta se distribuyen aproximadamente de forma normal. Sin embargo, en situaciones en las que los datos no son normales (por ejemplo, modelos basados en distribuciones binomial, Poisson o gamma), se requiere una metodología de estimación más general: la estimación de máxima verosimilitud.
La máxima verosimilitud consiste en elegir los estimadores de los parámetros desconocidos que maximizan la función de verosimilitud, es decir, la función de densidad o probabilidad de los datos observados, considerada como función de los parámetros. Esta metodología incluye los mínimos cuadrados como un caso especial (cuando los errores son normales).
Supongamos que \(y_1,\dots,y_n\) son observaciones independientes de una distribución exponencial con parámetro de escala \(\theta\). La función de densidad de la distribución exponencial es: \[ \mathcal{P}(y;\theta)=\theta\exp(-\theta y). \]
La función de verosimilitud conjunta para los \(n\) datos es: \[ \mathcal{L}(\theta;y)=\prod_{i=1}^{n}\mathcal{P}(y_i;\theta)=\theta^n\exp\Big(-\theta\sum_{i=1}^{n}y_i\Big)=\theta^n\exp(-n\bar{y}\theta), \] donde \(\bar{y}\) es la media aritmética de los \(y_i\).
El principio de máxima verosimilitud consiste en elegir el valor \(\hat{\theta}\) que maximice \(\mathcal{L}(\theta;y)\). Es fácil demostrar que en este caso la función se maximiza cuando: \[ \hat{\theta}=\frac{1}{\bar{y}}. \]
En la práctica es más conveniente trabajar con la función logarítmica de la verosimilitud, ya que la transformación logarítmica convierte productos en sumas: \[ \ell(\theta;y)=\log \mathcal{L}(\theta;y)=\sum_{i=1}^{n}\log \mathcal{P}(y_i;\theta). \] Para el ejemplo exponencial, el log-verosimilitud es: \[ \ell(\theta;y)=n(\log \theta-\bar{y}\theta). \] Maximizar \(\ell(\theta;y)\) es equivalente a maximizar \(\mathcal{L}(\theta;y)\).
Se puede demostrar que los mínimos cuadrados son un caso especial de la máxima verosimilitud. Por ejemplo, considere un modelo de regresión normal: \[ y_i\sim N(\mu_i,\sigma^2),\quad \mu_i=\beta_0+\beta_1 x_{1i}+\cdots+\beta_p x_{pi}. \] La función de densidad para \(y_i\) es: \[ \mathcal{P}(y_i;\mu_i,\sigma^2)=\frac{1}{\sqrt{2\pi\sigma^2}}\exp\Biggl\{-\frac{(y_i-\mu_i)^2}{2\sigma^2}\Biggr\}. \] El logaritmo de esta función es: \[ \log \mathcal{P}(y_i;\mu_i,\sigma^2)=-\frac{1}{2}\log(2\pi\sigma^2)-\frac{(y_i-\mu_i)^2}{2\sigma^2}. \] Así, la log-verosimilitud para los \(n\) datos es: \[ \ell\bigl(\beta_0,\dots,\beta_p,\sigma^2;y\bigr) =-\frac{n}{2}\log(2\pi\sigma^2)-\frac{1}{2\sigma^2}\sum_{i=1}^{n}(y_i-\mu_i)^2, \] donde el término \(\sum_{i=1}^{n}(y_i-\mu_i)^2\) es la suma de cuadrados (RSS). Para un valor fijo de \(\sigma^2\), maximizar la log-verosimilitud equivale a minimizar la RSS, lo que demuestra que los mínimos cuadrados son equivalentes a la estimación de máxima verosimilitud cuando se asume normalidad.
En este ejemplo se registra la cantidad total de lluvia de julio (en milímetros) en Quilpie, Australia, junto con el valor del índice de oscilación del sur (SOI) mensual. El SOI es la diferencia estandarizada entre las presiones de aire en Darwin y Tahiti y se ha relacionado con la lluvia en partes de Australia. Algunos agricultores pueden retrasar la siembra hasta que se supere un cierto umbral de lluvia dentro de un periodo (una “ventana de lluvia”).
Para modelar este fenómeno, se define la variable respuesta \(y\) de la siguiente forma: \[ y = \begin{cases} 1 & \text{si la lluvia total de julio excede 10 mm}, \\ 0 & \text{en caso contrario.} \end{cases} \]
El parámetro desconocido es la probabilidad de que la lluvia supere los 10 mm, que denotamos por \(\mu\) (ya que \(E[y] = \mu = \Pr(y=1)\)). Por el momento se ignora el efecto del SOI y se considera que todas las observaciones son equivalentes.
La función de probabilidad de \(y\) se define como: \[
\mathcal{P}(y;\mu) = \mu^y (1-\mu)^{1-y},
\] para \(y = 0\) o \(1\). Esta es la distribución Bernoulli con probabilidad \(\mu\). En R, la función dbinom() evalúa esta función de probabilidad; fijando size=1, la distribución binomial se reduce a la distribución Bernoulli.
Se evalúa el log-verosimilitud para algunos valores candidatos de \(\mu\), lo que muestra que el estimador de máxima verosimilitud (MLE) de \(\mu\) se encuentra cerca de 0.5. El código empleado es el siguiente:
library(GLMsData)
data(quilpie); names(quilpie)[1] "Year" "Rain" "SOI" "Phase" "Exceed" "y" mu <- c(0.2, 0.4, 0.5, 0.6, 0.8) # Valores candidatos para mu
ll <- rep(0, 5) # Vector para guardar los log-verosimilitudes
for (i in 1:5){
ll[i] <- sum(dbinom(quilpie$y, size=1, prob=mu[i], log=TRUE))
}
data.frame(Mu = mu, LogLikelihood = ll) Mu LogLikelihood
1 0.2 -63.69406
2 0.4 -48.92742
3 0.5 -47.13401
4 0.6 -48.11649
5 0.8 -60.92148# Generar una secuencia de valores candidatos para mu (por ejemplo, de 0.1 a 0.9)
mu_vals <- seq(0.1, 0.9, length.out = 200)
# Calcular la función de log-verosimilitud para cada valor de mu
ll_vals <- sapply(mu_vals, function(mu) sum(dbinom(quilpie$y, size = 1, prob = mu, log = TRUE)))
# Convertir el log-verosimilitud a verosimilitud (exponenciando)
likelihood_vals <- exp(ll_vals)
# Graficar la función de verosimilitud y el log-verosimilitud en paneles separados
par(mfrow = c(2, 1), mar = c(4, 4, 2, 1)) # Configurar dos paneles verticales
# Panel superior: Gráfico de la función de verosimilitud
plot(mu_vals, likelihood_vals, type = "l", lwd = 2, col = "blue",
xlab = expression(mu), ylab = "Verosimilitud",
main = "Función de Verosimilitud")
abline(v = 0.5147, lty = 2, col = "red") # Línea vertical en el MLE
# Panel inferior: Gráfico de la función log-verosimilitud
plot(mu_vals, ll_vals, type = "l", lwd = 2, col = "blue",
xlab = expression(mu), ylab = "Log-Verosimilitud",
main = "Función Log-Verosimilitud")
abline(v = 0.5147, lty = 2, col = "red") # Línea vertical en el MLE
# Restablecer la distribución de los gráficos a una sola figura (opcional)
par(mfrow = c(1, 1))La figura anterior muestra el gráfico de la función de verosimilitud (panel superior) y el log-verosimilitud (panel inferior) para un rango mayor de valores de \(\mu\). Visualmente, el MLE de \(\mu\) aparece justo por encima de 0.5, y se señala una línea vertical en \(\hat{\mu} = 0.5147\).
En esta sección se desarrolla un enfoque sistemático para maximizar el log-verosimilitud usando cálculo, y se introduce el concepto de score.
n <- length(quilpie$y)
muhat <- mean(quilpie$y)
Info <- n / (muhat * (1 - muhat))
c(muhat = muhat, FisherInfo = Info) muhat FisherInfo
0.5147059 272.2354978 1/sqrt(Info)[1] 0.06060767Esto indica que el estimador \(\hat{\mu}\) tiene un error estándar aproximado de 0.0606.
En esta sección se extiende la metodología de máxima verosimilitud al caso en que se estiman múltiples parámetros en un modelo de regresión. El proceso involucra la construcción de las ecuaciones score para cada parámetro, el cálculo de la información (observada y esperada), y la obtención de los errores estándar de los estimadores.
Modelo General:
Se asume que cada observación \(y_i\) sigue una distribución parametrizada por un parámetro de localización (media) \(\mu_i\) y un parámetro de dispersión \(\phi\). La media se define en función de las covariables \(x_{ij}\) y los parámetros de regresión \(\beta_j\) a través del predictor lineal
\[
\eta_i=\beta_0+\beta_1x_{i1}+\cdots+\beta_px_{ip}
\] y mediante una función de enlace que relaciona \(\mu_i\) y \(\eta_i\), es decir, \(g(\mu_i)=\eta_i\).
Función Log-Verosimilitud y Funciones Score:
La log-verosimilitud para \(n\) observaciones es \[
\ell(\beta_0,\beta_1,\dots,\beta_p;y) = \sum_{i=1}^{n} \log \mathcal{P}(y_i;\mu_i,\phi).
\] Las funciones score para cada parámetro se definen como: \[
U(\beta_j)=\frac{\partial \ell(\beta_0,\beta_1,\dots,\beta_p;y)}{\partial \beta_j}
= \sum_{i=1}^{n} \frac{\partial \log \mathcal{P}(y_i;\mu_i,\phi)}{\partial \mu_i}\frac{\partial \mu_i}{\partial \beta_j}.
\] Así, hay una ecuación score por cada parámetro desconocido en el modelo, y se resuelven conjuntamente para obtener los estimadores de máxima verosimilitud.
Información Observada:
La segunda derivada negativa de la log-verosimilitud respecto a los parámetros mide la “curvatura” del log-verosimilitud alrededor del máximo y se define para el par de parámetros \((\beta_j, \beta_k)\) como: \[
\mathcal{J}_{jk}(\beta) = -\frac{\partial^2 \ell(\beta;y)}{\partial \beta_j \partial \beta_k}
= -\frac{\partial U(\beta_j)}{\partial \beta_k}.
\]
Información Esperada o Fisher Information:
Se define como \[
\mathcal{I}_{jk}(\beta) = E[\mathcal{J}_{jk}(\beta)].
\] En particular, la cantidad \(\mathcal{I}_{jj}(\beta)\) es una medida del monto de información disponible sobre el parámetro \(\beta_j\) y, de forma general, la varianza de un estimador de máxima verosimilitud se aproxima por: \[
\operatorname{var}[\hat{\beta}_j] \approx \frac{1}{\mathcal{I}_{jj}(\beta)}.
\]
Utilizando la información esperada se puede obtener una aproximación de la precisión de los estimadores. Específicamente, el error estándar (la desviación típica) del estimador de \(\beta_j\) se estima como: \[ \operatorname{se}(\hat{\beta}_j) \approx \frac{1}{\sqrt{\mathcal{I}_{jj}(\hat{\beta})}}. \]
Este resultado muestra que la precisión del MLE está inversamente relacionada con la cantidad de información disponible; mientras mayor sea \(\mathcal{I}_{jj}\), menor será el error estándar del estimador.
En esta sección se extiende el método de máxima verosimilitud al contexto de varios parámetros y se lo expresa usando álgebra matricial. Esto facilita el manejo de modelos complejos y permite utilizar métodos iterativos para resolver las ecuaciones score.
Se asume que las respuestas provienen de una distribución con función de probabilidad \[ \mathcal{P}(\mathbf{y};\boldsymbol{\zeta}), \] donde \(\boldsymbol{\zeta} = [\zeta_1, \dots, \zeta_q]\) es el vector de parámetros desconocidos. La función de verosimilitud es la función de probabilidad vista como función de los parámetros: \[ \mathcal{L}(\boldsymbol{\zeta};\mathbf{y}) = \mathcal{P}(\mathbf{y};\boldsymbol{\zeta}). \] La log-verosimilitud se define como: \[ \ell(\boldsymbol{\zeta};\mathbf{y}) = \log \mathcal{L}(\boldsymbol{\zeta};\mathbf{y}). \] Los valores de los parámetros que maximizan la verosimilitud son los estimadores de máxima verosimilitud (MLE), denotados por \(\hat{\boldsymbol{\zeta}}\).
La derivada primera de la log-verosimilitud con respecto a \(\boldsymbol{\zeta}\) se conoce como el vector score: \[ U(\boldsymbol{\zeta}) = \frac{\partial \ell(\boldsymbol{\zeta};\mathbf{y})}{\partial \boldsymbol{\zeta}} = \sum_{i=1}^{n} \frac{\partial \log \mathcal{P}(y_i; \boldsymbol{\zeta})}{\partial \boldsymbol{\zeta}}. \] Cada componente \(U(\zeta_j)\) de este vector es la derivada parcial con respecto a \(\zeta_j\). El MLE se obtiene resolviendo la ecuación de score: \[ U(\hat{\boldsymbol{\zeta}}) = \mathbf{0}. \] En el caso de modelos de regresión, a menudo el interés se centra en los parámetros de la regresión. Por ejemplo, si el predictor lineal está dado por \[ \boldsymbol{\eta} = \beta_0 + \beta_1 x_{i1} + \cdots + \beta_p x_{ip} = X\boldsymbol{\beta}, \] y se define el enlace mediante \(\mu = g^{-1}(\eta)\), entonces se pueden escribir las ecuaciones score para cada \(\beta_j\) como: \[ U(\beta_j) = \frac{d\ell(\boldsymbol{\beta};\mathbf{y})}{d\mu} \frac{\partial \mu}{\partial \beta_j}. \]
La información se mide a través de la segunda derivada de la log-verosimilitud. Se define la información observada como: \[ \mathcal{J}(\boldsymbol{\zeta}) = -\frac{\partial^2 \ell(\boldsymbol{\zeta}; \mathbf{y})}{\partial \boldsymbol{\zeta} \partial \boldsymbol{\zeta}^T}. \] La información esperada (o información de Fisher) es: \[ \mathcal{I}(\boldsymbol{\zeta}) = E[\mathcal{J}(\boldsymbol{\zeta})]. \] En el contexto de modelos de regresión, con el predictor lineal \(\eta = X\boldsymbol{\beta}\), la información se puede expresar en términos de derivadas parciales con respecto a los parámetros y usualmente resulta en una matriz de dimensión \(p' \times p'\). Por ejemplo, la entrada \((j,k)\) se expresa como: \[ \mathcal{J}_{jk}(\boldsymbol{\beta}) = -\frac{\partial^2 \ell(\boldsymbol{\beta};\mathbf{y})}{\partial \beta_j \partial \beta_k}. \] Las propiedades importantes son: 1. La esperanza del vector score es cero: \(E[U(\boldsymbol{\zeta})] = \mathbf{0}\). 2. La varianza del vector score es igual a la información esperada: \[ \operatorname{var}[U(\boldsymbol{\zeta})] = \mathcal{I}(\boldsymbol{\zeta}) = E[U(\boldsymbol{\zeta})U(\boldsymbol{\zeta})^T]. \]
A partir de la matriz de información, se puede aproximar la varianza de cada estimador mediante el inverso de la matriz: \[ \operatorname{var}[\hat{\beta}_j] \approx (\mathcal{I}^{-1}(\boldsymbol{\beta}))_{jj}. \] Por lo tanto, el error estándar de \(\hat{\beta}_j\) se estima como: \[ \operatorname{se}(\hat{\beta}_j) \approx \frac{1}{\sqrt{(\mathcal{I}^{-1}(\hat{\boldsymbol{\beta}}))_{jj}}}. \]
La estimación de máxima verosimilitud se obtiene al resolver \(U(\hat{\boldsymbol{\zeta}}) = \mathbf{0}\). Una técnica iterativa para hacer esto es el método de Newton-Raphson, que en forma matricial es: \[ \hat{\boldsymbol{\zeta}}^{(r+1)} = \hat{\boldsymbol{\zeta}}^{(r)} + \mathcal{J}(\hat{\boldsymbol{\zeta}}^{(r)})^{-1}U(\hat{\boldsymbol{\zeta}}^{(r)}). \] Debido a que la matriz de información observada puede ser complicada de calcular, se utiliza en su lugar la matriz de información esperada (Fisher scoring): \[ \hat{\boldsymbol{\zeta}}^{(r+1)} = \hat{\boldsymbol{\zeta}}^{(r)} + \mathcal{I}(\hat{\boldsymbol{\zeta}}^{(r)})^{-1}U(\hat{\boldsymbol{\zeta}}^{(r)}). \]
Esta es una forma más estable y a menudo más sencilla de implementar la estimación de máxima verosimilitud.