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.
El método de máxima verosimilitud no sólo ofrece un camino unificado para estimar parámetros en modelos normales y no‒normales: bajo supuestos regulares —básicamente, que la información aumenta al crecer el tamaño muestral \(n\)— sus estimadores disfrutan de atributos teóricos muy atractivos. A continuación se describen esas propiedades tanto para un único parámetro como para el caso multivariado, manteniendo la notación y las ecuaciones originales.
Supóngase que cada respuesta proviene de una distribución parametrizada por \(\zeta\) (o por el vector \(\boldsymbol{\zeta}\)). El estimador de máxima verosimilitud, denotado \(\hat{\zeta}\), se define como el valor que maximiza la función de verosimilitud (o su logaritmo) y, en consecuencia, coincide con la solución de la ecuación de “score” \(U(\hat\zeta)=0\). A diferencia de métodos basados en momentos, la verosimilitud utiliza la forma completa de la densidad y esto, junto con condiciones de regularidad suaves, confiere a \(\hat{\zeta}\) una serie de propiedades asintóticas que justifican su uso práctico.
Para un único parámetro \(\zeta\) las conclusiones clásicas son cinco.
- invarianza: cualquier transformación biyectiva \(s(\zeta)\) tiene por MLE el valor transformado \(s(\hat\zeta)\).
- insesgamiento asintótico: \(\mathrm{E}[\hat\zeta]=\zeta\) cuando \(n\to\infty\); con muestras pequeñas el sesgo puede no ser despreciable.
- eficiencia asintótica: ningún otro estimador insesgado asintótico posee varianza menor; si existiera, sería estadísticamente equivalente a \(\hat\zeta\).
- consistencia: \(\hat\zeta\rightarrow\zeta\) en probabilidad a medida que crece el tamaño muestral.
- normalidad asintótica: cuando \(\zeta_{0}\) es el valor verdadero,
\[ \hat{\zeta}\sim N\!\left(\zeta_{0},\,1/\mathcal{I}(\zeta_{0})\right),\qquad n\to\infty, \]
donde \(\mathcal{I}(\zeta_{0})\) es la información de Fisher. De aquí se deduce directamente
\[ \operatorname{var}[\hat{\zeta}]=1/\mathcal{I}(\zeta_{0}), \]
así que su error estándar se estima mediante \(\sqrt{\mathcal{I}(\zeta_{0})^{-1}}\).
Cuando el parámetro es un vector \(\boldsymbol{\zeta}=(\zeta_{1},\dots,\zeta_{q})\) las cinco propiedades se mantienen reinterpretadas en clave matricial.
La invarianza se refiere a transformaciones biyectivas de \(\boldsymbol{\zeta}\).
El insesgamiento asegura \(\mathrm{E}[\hat{\boldsymbol{\zeta}}]=\boldsymbol{\zeta}\) asintóticamente.
La eficiencia implica que ningún otro estimador vectorial insesgado posee matriz de varianzas–covarianzas más pequeña.
La consistencia equivale a \(\hat{\boldsymbol{\zeta}}\to\boldsymbol{\zeta}\) en probabilidad.
Por último, la normalidad asintótica adopta la forma
\[ \hat{\boldsymbol{\zeta}}\sim N_{q}\!\left(\boldsymbol{\zeta}_{0},\,\mathcal{I}(\boldsymbol{\zeta}_{0})^{-1}\right),\qquad n\to\infty. \]
En consecuencia,
\[ \operatorname{var}\!\bigl[\hat{\boldsymbol{\zeta}}\bigr]=\mathcal{I}(\boldsymbol{\zeta}_{0})^{-1},\qquad \operatorname{se}(\hat{\zeta}_{j})=\sqrt{\bigl[\mathcal{I}(\boldsymbol{\zeta}_{0})^{-1}\bigr]_{jj}}. \]
Una forma equivalente y útil para contrastes globales es
\[ \bigl(\hat{\boldsymbol{\zeta}}-\boldsymbol{\zeta}_{0}\bigr)^{\!\top}\! \mathcal{I}(\boldsymbol{\zeta}_{0})\, \bigl(\hat{\boldsymbol{\zeta}}-\boldsymbol{\zeta}_{0}\bigr)\;\sim\;\chi^{2}_{q}, \]
lo que sirve de base a pruebas de razón de verosimilitud y a la construcción de regiones de confianza elipsoidales.
Estas propiedades explican por qué la máxima verosimilitud domina la inferencia paramétrica moderna: brinda estimadores que, con datos suficientes, son insesgados, eficientes y de distribución conocida, permitiendo construir intervalos o realizar tests sin depender de aproximaciones ad hoc. Sin embargo, la “virtud asintótica” no garantiza un desempeño óptimo con muestras pequeñas; en tales casos puede ser prudente evaluar sesgos, utilizar ajustes (p. ej. bias‐reduction) o recurrir a métodos de remuestreo. Finalmente, la normalidad asintótica depende de que el modelo esté bien especificado; diagnósticos y contrastes de adecuación siguen siendo imprescindibles para que estas propiedades se traduzcan en inferencias fiables.
Una vez estimado el modelo mediante máxima verosimilitud surge, de manera natural, la necesidad de contrastar afirmaciones sobre sus parámetros. Esta sección describe cómo, bajo supuestos asintóticos, se construyen tres familias de contrastes equivalentes: Wald, score y razón de verosimilitud. La exposición comienza con un único parámetro y evoluciona hacia modelos multivariados, incluidos los contrastes sobre subconjuntos de parámetros.
Sea \(\zeta\) el único parámetro de interés y considérese la hipótesis nula \(H_{0}\!: \zeta=\zeta^{0}\) frente a la alternativa bilateral.
La prueba de Wald mide la distancia entre el estimador y el valor hipotético, escalada por la varianza asintótica de \(\hat\zeta\),
\[ W=\frac{(\hat\zeta-\zeta^{0})^{2}}{\widehat{\operatorname{var}}[\hat\zeta]},\qquad \widehat{\operatorname{var}}[\hat\zeta]=\frac{1}{\mathcal{I}(\hat\zeta)}. \]
Bajo \(H_{0}\), \(W\xrightarrow{d}\chi^{2}_{1}\). Emplear \(Z=\sqrt{W}\) resulta útil cuando se desea una formulación unilateral, porque \(Z\xrightarrow{d}N(0,1)\).
La prueba de score evita estimar \(\hat\zeta\); se concentra en la pendiente de la log‑verosimilitud evaluada en \(\zeta^{0}\):
\[ S=\frac{U(\zeta^{0})^{2}}{\mathcal{I}(\zeta^{0})},\qquad U(\zeta)=\frac{\mathrm d\ell}{\mathrm d\zeta}. \]
También aquí \(S\xrightarrow{d}\chi^{2}_{1}\) y \(\sqrt{S}\) converge a \(N(0,1)\), lo que permite contrastes unilaterales con la misma facilidad.
Finalmente, la prueba de razón de verosimilitud cuantifica la pérdida de ajuste que implica imponer la restricción:
\[ L=2\{\ell(\hat\zeta;y)-\ell(\zeta^{0};y)\}\xrightarrow{d}\chi^{2}_{1}. \]
Los tres estadísticos comparten, pues, la misma distribución límite y, con muestras grandes, conducen a conclusiones concordantes. Ello se ilustra con los datos de lluvia en Quilpie: al contrastar \(H_{0}\!:\mu=0.5\) los tres valores (0.0588, 0.0588, 0.0588) arrojan \(p\)‑valores próximos a 0.81, indicando que la proporción analizada no difiere del 50 %:
## --- datos -------------------------------------------------------------
## Se supone que ya existe el data‑frame `quilpie`
## con la variable binaria y (1 = lluvia > 10 mm; 0 = lo contrario)
y <- quilpie$y # vector de 0/1
n <- length(y) # tamaño muestral
## --- estimación de máxima verosimilitud -------------------------------
muhat <- mean(y) # MLE de µ en el modelo Bernoulli
var_hat <- muhat*(1 - muhat)/n # varianza asintótica de mû
## --- estadístico de Wald ----------------------------------------------
W <- (muhat - 0.5)^2 / var_hat
pW <- pchisq(W, df = 1, lower.tail = FALSE)
## --- estadístico de score ---------------------------------------------
U_mu0 <- (sum(y)/0.5)-((n-sum(y))/0.5) # U(mu0)
I_mu0 <- n / (0.5 * (1 - 0.5)) # I(mu0)
S <- U_mu0^2 / I_mu0 # S = U^2/I
pS <- pchisq(S, df = 1, lower.tail = FALSE)
## --- estadístico de razón de verosimilitud ----------------------------
logLik_mu0 <- sum(dbinom(y, size = 1, prob = 0.5, log = TRUE))
logLik_muhat <- sum(dbinom(y, size = 1, prob = muhat, log = TRUE))
L <- 2 * (logLik_muhat - logLik_mu0)
pL <- pchisq(L, df = 1, lower.tail = FALSE)
## --- resumen ----------------------------------------------------------
round(
cbind(
Estadístico = c(W, S, L),
"p‑valor" = c(pW, pS, pL)
),
4
) Estadístico p‑valor
[1,] 0.0589 0.8083
[2,] 0.0588 0.8084
[3,] 0.0588 0.8084Si \(\boldsymbol{\zeta}\in\mathbb{R}^{q}\), la nula global \(H_{0}\!: \boldsymbol{\zeta}=\boldsymbol{\zeta}^{0}\) se examina mediante extensiones matriciales:
\[ \begin{aligned} \text{Wald}&:\;W=(\hat{\boldsymbol{\zeta}}-\boldsymbol{\zeta}^{0})^{\!\top}\! \mathcal{I}(\hat{\boldsymbol{\zeta}})\, (\hat{\boldsymbol{\zeta}}-\boldsymbol{\zeta}^{0}),\\ \text{Score}&:\;S=U(\boldsymbol{\zeta}^{0})^{\!\top}\! \mathcal{I}(\boldsymbol{\zeta}^{0})^{-1} U(\boldsymbol{\zeta}^{0}),\\ \text{Razón de verosimilitud}&:\; L=2\{\ell(\hat{\boldsymbol{\zeta}};y)-\ell(\boldsymbol{\zeta}^{0};y)\}, \end{aligned} \]
y cada uno converge a \(\chi^{2}_{q}\).
Con los mismos datos de Quilpie, probar \(H_{0}\!:\beta_{0}=\beta_{1}=0\) arroja distintos valores numéricos (11.79, 15.92, 18.37), pero todos conducen a \(p<0.01\), descartando la nula:
y <- quilpie$y # 0 / 1 (lluvia > 10 mm)
x <- quilpie$SOI # covariable – Southern Oscillation
n <- length(y)
X <- cbind(1, x) # matriz de diseño completa (β0, β1)
## ----------------------------------------------------------------------
## 2. Modelo completo (β0, β1 libres) y modelo nulo (β0 = β1 = 0)
## ----------------------------------------------------------------------
fit.full <- glm(y ~ x, family = binomial(link = "logit"))
beta_hat <- coef(fit.full) # (β̂0, β̂1)
## log‑likelihood del modelo completo (la función ‘logLik’ lo entrega)
ll_full <- as.numeric(logLik(fit.full))
## log‑likelihood bajo H0 (μ = 0.5 constante ⇒ logit(μ) = 0)
ll_null <- sum(dbinom(y, size = 1, prob = 0.5, log = TRUE))
## ----------------------------------------------------------------------
## 3. Estadístico Wald
## W = (β̂ - β0)ᵀ (Var β̂)⁻¹ (β̂ - β0) con β0 = (0,0)
## ----------------------------------------------------------------------
V <- vcov(fit.full) # Var‑Cov de β̂ ( = I(β̂)⁻¹ )
Wald <- t(beta_hat) %*% solve(V) %*% beta_hat
p.W <- pchisq(Wald, df = 2, lower.tail = FALSE)
## ----------------------------------------------------------------------
## 4. Estadístico Score (Rao)
## S = U(β0)ᵀ I(β0)⁻¹ U(β0)
## donde U(β0) = Xᵀ (y - μ0) con μ0 = 0.5
## I(β0) = Xᵀ W X, W = μ0 (1-μ0) I
## ----------------------------------------------------------------------
mu0 <- 0.5
U0 <- t(X) %*% (y - mu0) # vector score bajo la nula
Wmat <- diag(mu0 * (1 - mu0), n) # 0.25 * I
I0 <- t(X) %*% Wmat %*% X
Score <- t(U0) %*% solve(I0) %*% U0
p.S <- pchisq(Score, df = 2, lower.tail = FALSE)
## ----------------------------------------------------------------------
## 5. Cociente de verosimilitudes
## L = 2 [ ℓ(β̂) – ℓ(β0) ]
## ----------------------------------------------------------------------
L.Ratio <- 2 * (ll_full - ll_null)
p.L <- pchisq(L.Ratio, df = 2, lower.tail = FALSE)
## ----------------------------------------------------------------------
## 6. Resumen
## ----------------------------------------------------------------------
out <- rbind(
Wald = c(stat = Wald, p.value = p.W),
Score = c(stat = Score, p.value = p.S),
L.Ratio = c(stat = L.Ratio, p.value = p.L)
)
print( round(out, 4) ) stat p.value
Wald 11.7952 0.0027
Score 15.9248 0.0003
L.Ratio 18.3674 0.0001Con frecuencia se desea verificar si un subconjunto \(\boldsymbol{\zeta}_{2}\) (dimensión \(q_{2}\)) aporta información adicional una vez fijado el resto \(\boldsymbol{\zeta}_{1}\). Particionando la matriz de información,
\[ \mathcal{I}(\hat{\boldsymbol{\zeta}})= \begin{pmatrix} \mathcal{I}_{11} & \mathcal{I}_{12}\\ \mathcal{I}_{21} & \mathcal{I}_{22} \end{pmatrix}, \]
el contraste \(H_{0}\!:\boldsymbol{\zeta}_{2}=\boldsymbol{\zeta}^{0}_{2}\) se basa en
\[ W = (\hat{\boldsymbol{\zeta}}_{2}-\boldsymbol{\zeta}^{0}_{2})^{\!\top} (\mathcal{I}^{22})^{-1} (\hat{\boldsymbol{\zeta}}_{2}-\boldsymbol{\zeta}^{0}_{2}), \quad S = U(\boldsymbol{\zeta}^{*})^{\!\top}\! \mathcal{I}(\boldsymbol{\zeta}^{*})^{-1} U(\boldsymbol{\zeta}^{*}), \quad L = 2\{\ell(\hat{\boldsymbol{\zeta}})-\ell(\boldsymbol{\zeta}^{*})\}, \]
donde \(\boldsymbol{\zeta}^{*}\) combina los MLE de \(\boldsymbol{\zeta}_{1}\) con el valor fijado \(\boldsymbol{\zeta}^{0}_{2}\). Todos los estadísticos tienden a \(\chi^{2}_{q_{2}}\). En el ejemplo, contrastar \(\beta_{1}=0\) (con \(q_{2}=1\)) produce \(p\)‑valores de orden \(10^{-4}\) con las tres técnicas, confirmando la relevancia de la covariable SOI.
###############################################################
## Quilpie: prueba H0 : β1 = 0 (β0 libre) — q2 = 1
###############################################################
## --- Modelo completo (ya calculado antes) -------------------
beta_hat <- coef(fit.full) # (β̂0, β̂1)
V <- vcov(fit.full) # Var‑Cov(β̂)
ll_full <- as.numeric(logLik(fit.full))
## --- Modelo reducido: solo intercepto -----------------------
fit.red <- glm(y ~ 1, family = binomial, data = quilpie)
ll_red <- as.numeric(logLik(fit.red))
## ============================================================
## 1) Wald (Z²)
## ============================================================
var_b1 <- V[2, 2]
Wald_1 <- beta_hat[2]^2 / var_b1 # χ²1
pW_1 <- pchisq(Wald_1, df = 1, lower.tail = FALSE)
## ============================================================
## 2) Score (Rao) — se evalúa en β1 = 0
## ============================================================
mu0 <- fitted(fit.red) # μ con β1 = 0
U1 <- sum((y - mu0) * x) # componente score para β1
I22 <- sum(mu0 * (1 - mu0) * x^2) # elemento (2,2) de I(β0,β1)
Score_1 <- U1^2 / I22
pS_1 <- pchisq(Score_1, df = 1, lower.tail = FALSE)
## ============================================================
## 3) Cociente de verosimilitudes
## ============================================================
LR_1 <- 2 * (ll_full - ll_red) # χ²1
pL_1 <- pchisq(LR_1, df = 1, lower.tail = FALSE)
## --- Resumen ------------------------------------------------
round(
rbind(Wald = c(stat = Wald_1, p = pW_1),
Score= c(stat = Score_1, p = pS_1),
L.Ratio=c(stat = LR_1, p = pL_1)), 5) stat.x p.x
Wald 11.66136 0.00064
Score 15.87943 0.00007
L.Ratio 18.30858 0.00002Cuando \(q_{2}=1\) la prueba de Wald adopta la forma escalar familiar
\[ W=\frac{(\hat{\beta}_{j}-\beta^{0}_{j})^{2}} {\operatorname{var}[\hat{\beta}_{j}]}, \qquad Z=\frac{\hat{\beta}_{j}-\beta^{0}_{j}} {\sqrt{\operatorname{var}[\hat{\beta}_{j}]}}. \]
En términos prácticos, \(Z\) es el estadístico que suelen devolver los paquetes de software junto con su error estándar. El contraste de razón de verosimilitud requiere volver a ajustar el modelo con la restricción impuesta, mientras que la prueba score evita ese paso calculando la pendiente en el punto nulo.
Los tres métodos son asintóticamente equivalentes, pero difieren en comodidad computacional y estabilidad para muestras pequeñas. El contraste de Wald sólo necesita el vector de estimadores y su matriz de varianzas, por lo que resulta muy simple; sin embargo, su comportamiento puede deteriorarse cuando la estimación está cerca del borde del espacio paramétrico (fenómeno de Hauck–Donner). La prueba de score es atractiva si el ajuste bajo \(H_{0}\) es trivial y, además, suele mantener mejor nivel con n moderado. La razón de verosimilitud, a costa de dos ajustes, es la opción preferida cuando se desea un criterio globalmente robusto y compatible con comparaciones entre modelos anidados. En la práctica conviene contrastar los resultados y, ante discrepancias con muestras pequeñas, favorecer score o razón de verosimilitud.
La estimación puntual mediante máxima verosimilitud resulta muy potente; sin embargo, la práctica exige acompañar cada punto con una medida explícita de precisión. En este apartado describimos cómo se construyen regiones o bandas de confianza a partir de los tres grandes estadísticos asintóticos —Wald, puntaje y cociente de verosimilitudes— haciendo hincapié en sus vínculos con la información de Fisher. En todos los casos se parte de la aproximación cuadrática de la log‑verosimilitud cuando el tamaño muestral \(n\) es grande, suposición bajo la cual los estadísticos convergen a distribuciones \(\chi^{2}\).
Supóngase que \(\boldsymbol{\zeta}\) contiene \(q\) parámetros desconocidos. Una región de confianza conjunta al nivel \(100(1-\alpha)\%\) se obtiene hallando todos los valores de \(\boldsymbol{\zeta}\) que satisfacen, respectivamente,
\[ (\hat{\boldsymbol{\zeta}}-\boldsymbol{\zeta})^{\mathsf T}\, \mathcal I(\hat{\boldsymbol{\zeta}})\, (\hat{\boldsymbol{\zeta}}-\boldsymbol{\zeta}) \;\le\; \chi^{2}_{q,\,1-\alpha}, \]
\[ U(\boldsymbol{\zeta})^{\mathsf T}\, \mathcal I(\boldsymbol{\zeta})^{-1}\, U(\boldsymbol{\zeta}) \;\le\; \chi^{2}_{q,\,1-\alpha}, \]
\[ 2\{\ell(\hat{\boldsymbol{\zeta}};\mathbf y)-\ell(\boldsymbol{\zeta};\mathbf y)\} \;\le\; \chi^{2}_{q,\,1-\alpha}. \]
La primera expresión corresponde al criterio de Wald, la segunda al de puntaje y la tercera al cociente de verosimilitudes. Hallar la envolvente exacta de estos conjuntos rara vez es trivial, pues implica resolver ecuaciones no lineales en dimensión \(q\); por ello, en la práctica se suelen recurrir a cortes bidimensionales o aproximaciones elípticas.
Cuando el interés recae en un componente \(\zeta_{j}\), las expresiones anteriores se reducen a funciones unidimensionales y es posible localizar numéricamente los dos valores frontera que satisfacen el umbral \(\chi^{2}_{1,\,1-\alpha}\). La figura 4.9 del Dunn ilustra la geometría: los límites Wald se obtienen desplazándose una distancia fija a ambos lados de \(\hat{\zeta}_{j}\), mientras que los límites puntaje y verosimilitud se definen allí donde la pendiente del log‑likelihood o la caída respecto al máximo alcanzan la cota crítica.
Para el estadístico de Wald, la forma cerrada es particularmente cómoda:
\[ \hat{\zeta}_{j}\; \pm\; z_{1-\alpha/2}\; \sqrt{\operatorname{var}[\hat{\zeta}_{j}]} , \]
donde \(z_{1-\alpha/2}\) es el cuantil estándar normal y \(\operatorname{var}[\hat{\zeta}_{j}]\) se obtiene como la entrada \((j,j)\) de \(\mathcal I(\hat{\boldsymbol{\zeta}})^{-1}\). Nótese la simetría inherente, una característica que no comparten las bandas basadas en puntaje o en el cociente de verosimilitudes, ya que la log‑verosimilitud rara vez es perfectamente cuadrática en muestras reales.
En el ejemplo de Quilpie (sección 4.11), el parámetro de interés es la pendiente \(\beta_{1}\). Los tres métodos arrojan límites muy similares: el intervalo Wald va de \(0.0624\) a \(0.2305\), el de puntaje de \(0.0655\) a \(0.2289\) y el de verosimilitud de \(0.0719\) a \(0.2425\). Las pequeñas discrepancias obedecen a la ligera asimetría de la superficie de verosimilitud, pero todos los intervalos excluyen el cero, reforzando la significación detectada previamente por los contrastes.
Los contrastes de razón de verosimilitudes funcionan cuando un modelo está contenido en otro, pero fallan si ambos modelos son no anidados. En ese contexto la práctica estándar consiste en comparar AIC y BIC, dos criterios que parten de la log‑verosimilitud máxima pero añaden un castigo por complejidad.
Sea \(\ell(\hat{\zeta}_{1},\ldots,\hat{\zeta}_{p};y)\) la log‑verosimilitud evaluada en los MLE de un modelo cualquiera. Entonces
\[ \mathrm{AIC}\;=\;-2\,\ell(\hat{\boldsymbol{\zeta}};y)\;+\;2p, \]
mientras que
\[ \mathrm{BIC}\;=\;-2\,\ell(\hat{\boldsymbol{\zeta}};y)\;+\;(\log n)\,p. \]
La primera expresión aplica un penalizador constante \(k=2\); la segunda emplea \(k=\log n\), de modo que el castigo crece con el tamaño muestral. En ambos casos cuanto más negativo sea el valor resultante, mejor se considera el equilibrio entre ajuste y parsimonia.
En la base de lluvias de Quilpie se han ajustado dos modelos binomiales para la probabilidad de superar 10 mm en julio:
Modelo 1: logit \((\mu)=\beta_{0}+\beta_{1}\,\text{SOI}\).
La log‑verosimilitud máxima es \(-39.9503\) y hay dos parámetros libres (\(p=2\)). Se obtiene
AIC = 79.9 y BIC = 84.3.
Modelo 2: logit \((\mu)=\beta_{0}+\sum_{k=2}^{5}\beta_{k}\,\text{Fase}_{k}\),
donde la SOI se reemplaza por sus cinco fases categóricas (cuatro dummies más intercepto, \(p=5\)). La log‑verosimilitud asciende a \(-31.8995\), produciendo
AIC = 75.8 y BIC = 86.9.
Las conclusiones difieren: