Estimación via máxima verosimilitud

Especialización en Estadística Aplicada

Roberto Trespalacios

• El estimador \( \hat{\boldsymbol{\beta}}_{GLM} \) satisface las condiciones de optimalidad del Teorema de Gauss-Markov cuando \( \boldsymbol{V} \) es

• En general, la matriz \( \boldsymbol{V} \) es desconocida y se debe estimar con los datos disponibles.

• En este caso el estimador \( \hat{\boldsymbol{\beta}}_{GLM} \) no es insesgado y su varianza no tiene una fórmula exacta.

• La estimación de los parámetros en el modelo demanda métodos iterativos que oscilan entre la estimación de \( \boldsymbol{\beta} \) y \( \boldsymbol{V} \) hasta lograr convergencia.

• Es importante notar que el estimador de mínimos cuadrados ordinarios \( \hat{\boldsymbol{\beta}}_{LS} \) se puede usar para estimar \( \boldsymbol{\beta} \) en el modelo \( \boldsymbol{Y} = \boldsymbol{X\beta} + \boldsymbol{\varepsilon} \); ya que es insesgado y consistente.

• El estimador es menos eficiente \( \hat{\boldsymbol{\beta}}_{LS} \) que el estimador \( \hat{\boldsymbol{\beta}}_{GLM} \).

• Las implicaciones de este resultado es que en la práctica los errores estándar del estimador LS ordinario podrían ser mayores que los del estimador GLS.

• Note que los estimadores coinciden cuando \( \boldsymbol{V} = \boldsymbol{I} \).

Recordemos que la función de densidad(probabilidad) para la distribución normal viene dada por:

\[ P(x) = \frac{1}{\sqrt{2\pi}\sigma}e^{-(x-\mu)^2/(2\sigma^2)} \]

• Tenemos los datos a los que nos gustaría ajustar el PDF, por lo que conocemos el vector \( x \).

• Necesitamos encontrar \( \mu \) y \( \sigma \).

• Como se conoce el vector \( x \), variaremos los parámetros \( \mu \) y \( \sigma \) para maximizar el valor de la función de verosimilitud \( L(\mu,sigma) \).

\[ L(x|\mu,\sigma)= \frac{1}{\sqrt{2\pi}\sigma}e^{-(x-\mu)^2/(2\sigma^2)} \]

Supongamos que las variables \( x_1,x_2, \dots,x_n \) tienen distribución normal con parametros \( \mu \) y \( \sigma \). entonces la función de verosimilitud toma la forma: \( L(x_1,x_2, \dots,x_n\mu,\sigma) \). Por simplicidad solo escribiremos \( L(\mu,\sigma) \).

\[ \begin{align*} L(\mu, \sigma) &= \frac{1}{\sqrt{2\pi}\sigma}e^{-(x_1-\mu)^2/(2\sigma^2)} \times \frac{1}{\sqrt{2\pi}\sigma}e^{-(x_2-\mu)^2/(2\sigma^2)} \times \cdots \times \frac{1}{\sqrt{2\pi}\sigma}e^{-(x_n-\mu)^2/(2\sigma^2)}\\ &= \prod_{i=1}^n \frac{1}{\sqrt{2\pi}\sigma}e^{-(x_i-\mu)^2/(2\sigma^2)}\\ &= \left[\frac{1}{\sqrt{2\pi}\sigma}\right]^{n} \prod_{i=1}^n e^{-(x_i-\mu)^2/(2\sigma^2)}\\ &= \frac{(2\pi)^{-2/n}}{\sigma^n} \exp\left[ -\frac{ \sum_{i=1}^n(x_i-\mu)^2}{2\sigma^2}\right] \end{align*} \]

Continuamos tomando el logaritmo de la función de verosimilitud porque es más fácil trabajar con sumas en lugar de productos.

\[ \ln[L(\mu,\sigma)] = -\frac{n}{2}\ln(2\pi) - \frac{n}{2}\ln(\sigma^2) - \frac{1}{2\sigma^2}\sum_{i=1}^{n} (x_i-\mu)^2 \]

El primer término es una constante podemos descomponer el logaritmo. Del mismo modo hacemos con \( n \) del numerador en el segundo término. Nos quedamos con

\[ \ln[L(\mu,\sigma)] = - \frac{1}{2}\ln(\sigma^2) - \frac{1}{2\sigma^2}\sum_{i=1}^{n} (x_i-\mu)^2 \]

Derivando \( \ln[L(\mu,\sigma)] \) respecto a \( mu \) y \( \sigma \) y luego despejando para, \( \mu \) y \( \sigma \) tenemos que:

\[ \frac{\partial (\ln L)}{\partial \mu} = \frac{ \sum_{i=1}^n(x_i-\mu)}{\sigma^2} =0 \]

luego \[ \hat\mu= \frac{\sum_{i=1}^n x_i}{n} \].

Similarmente, para \( \sigma \) tenemos que:

\[ \frac{\partial (\ln L)}{\partial \sigma} = -\frac{n}{\sigma} - \frac{ \sum_{i=1}^n(x_i-\mu)^2}{\sigma^3} =0 \]

luego,

\[ \hat\sigma=\sqrt{\frac{\sum_{i=1}^n(x_i-\mu)^2}{n}} \].

De acuerdo con el modelo clásico (\( \boldsymbol{Y} = \boldsymbol{X\beta} + \boldsymbol{\varepsilon} \)) las observaciones \( Y_i \) son independientes y normalmente distribuidas. Por lo tanto, se puede definir la función de verosimilitud:

\[ L(\boldsymbol{\beta}, \sigma^2; \boldsymbol{y})=\frac{1}{(2\pi\sigma^2)^{n/2}}\prod_{i=1}^n\exp\left[ -\frac{(y_i-\boldsymbol{x}'_i\boldsymbol{\beta})^2}{2\sigma^2}\right] \]

Los estimadores de máxima verosimilitud de los parámetros del modelo están dados por

\[ \hat{\boldsymbol{\beta}}_{ML} = (\boldsymbol{X}'\boldsymbol{X})^{-1}\boldsymbol{X}'\boldsymbol{Y} \]

y

\[ \hat{\sigma}_{ML}=\sqrt{\frac{\sum_i^n(y_i-\hat{y}_i)^2}{n}} \]

Podemos ver, que para

\[ \hat{\boldsymbol{\beta}}_{ML} = (\boldsymbol{X}'\boldsymbol{X})^{-1}\boldsymbol{X}'\boldsymbol{Y} \]

y

\[ \hat{\sigma}_{ML}=\sqrt{\frac{\sum_i^n(y_i-\hat{y}_i)^2}{n}} \]

tenemos que:

• El estimador ML de \( \boldsymbol{\beta} \) es igual que el estimador LS pero el estimador ML de \( \sigma^2 \) tiene un sesgo negativo (anteriormente vimos que \( s^2 \) es un estimador insesgado para \( \sigma^2 \)).

• Este sesgo se debe a que la estimación de \( \boldsymbol{\beta} \) no tiene en cuenta al calcular \( \sigma^2_{ML} \) (no ajusta el denominador o los grados de libertad).

El problema con el método de ML se puede corregir si se usan los datos transformados que no dependan de entonces se puede definir el logaritmo de la los coeficientes del modelo. Esto es, si definimos \( r_i = y_i - \boldsymbol{x}'_i\hat{\boldsymbol\beta} \) función de verosimilitud basada en \( r_i \) la cual se conoce como función de verosimilitud restringida (REML).

\[ l_{REML}(\sigma^2;\boldsymbol y) = \frac{n-p}{2}log(\sigma^2)-\frac{1}{2\sigma^2}\sum_{i=1}^nr_i^2 \]

y

\[ \hat{\sigma}^2_{REML}=\frac{\sum_i^n(y_i-\hat{y}_i)^2}{n-p} \]

el cual coincide con el estimador de mínimos cuadrados. El estimador REML para \( \boldsymbol \beta \) es igual al estimador ML. Estas coincidencias son ciertas para el caso del modelo clásico (observaciones independientes y varianza constante) pero no surgen en modelos más complejos.

La matriz de covarianza del estimador de los coeficientes en el modelo clásico esta dada por

\[ Var(\hat {\boldsymbol \beta}) = \sigma^2 (\boldsymbol{X}'\boldsymbol{X})^{-1} \]

Un estimador para la matriz de covarianza esta dado por

\[ \hat{Var}(\hat {\boldsymbol \beta}) = \hat{\sigma}^2 (\boldsymbol{X}'\boldsymbol{X})^{-1} \]

donde \( \hat{\sigma}^2 \) puede ser el estimador LS, ML o REML de \( \sigma^2 \), dependiendo del método de estimación usado. También note que a pesar de que \( \hat{\boldsymbol\beta} \) basado en ML o REML son iguales sus matrices de covarianza estimadas no son iguales debido a la estimación de \( \sigma^2 \).

Usando la normalidad de \( \hat{\boldsymbol\beta} \) es posible llevar a cabo inferencia acerca de los parámetros en el modelo (\( \boldsymbol\beta \) y \( \sigma^2 \)). Esto incluye tanto llevar a cabo pruebas de hipótesis como intervalos de confianza usando las distribuciones clásicas derivadas de la normalidad tales como \( t \), \( \chi^2 \) y \( F \).

Métodos Tradicionales y Modernos de Selección de Variables

La mayoría de libros tradicionales de regresión hacen referencia a las técnicas de selección de variables tales como:

• Forward

• Backward

• Stepwise

entre otras. A pesar de sus críticas, estas técnicas están implementadas en la mayoría de programas estadísticos tanto comerciales como gratuitos.

• Algunas técnicas modernas de selección de modelos usan métodos de penalización para estimar los parámetros del modelo para aquellas variables que no son significativas las estimaciones de los coeficientes para estas variables son iguales a cero.

• Por ejemplo, el famoso método o estimador LASSO (Least Absolute Selection and Shrinkage Operator) estima los coeficientes y al mismo tiempo selecciona variables al arrojar coeficientes iguales a 0 para las variables no relevantes en el modelo.

• Para el modelo clásico el LASSO básicamente resuelve el siguiente problema que involucra la norma \( L_1 \) o valor absoluto en la penalidad de los coeficientes

\[ Q =\min_{\beta}\left(\frac{1}{2}SSE + \lambda \sum_k |\beta_k| \right) \]

donde \( \lambda \) es una constante de regularización también llamado un parámetro de ajuste (tuning parameter). Cuando \( \lambda = 0 \) entonces el problema de la función de minimización \( Q \), se reduce al método de mínimos cuadrados tradicional.

• Otros ejemplos de estos métodos de regularización usan otras funciones de penalización tales como SCAD, LASSO adaptivo, entre otros.

• Estos métodos están especialmente diseñados para problemas con dimensiones grandes en los cuales \( n << p \). Este tipo de problemas se ven frecuentemente en genética.

• El método LASSO y algunos de sus semejantes se encuentran implementados en algunas librerías de R.

Veamos si los datos Hitters estan completos (no tienen NA)

1 - (sum(complete.cases(Hitters)) / nrow(Hitters))

[1] 0.18

El dataset contiene aproximadamente un 18\% de observaciones a las que les falta al menos una variable (missing value). Existen técnicas para manejar missing values pero en este caso simplemente se van a excluir los registros que no estén completos.

datos.Hitters <- na.omit(Hitters)
dim(datos.Hitters)

[1] 263  20

los datos Boston de la librería MASS, contiene las siguientes columnas:

• crim: tasa de criminalidad per cápita por ciudad.
• zn: proporción de tierra residencial dividida en zonas para lotes de más de 25,000 pies cuadrados.
• indus: proporción de acres de negocios no minoristas por pueblo.
• chas: Variable dummy de Charles River (= 1 si el tramo limita el río, 0 de lo contrario).
• nox: concentración de óxidos de nitrógeno (partes por 10 millones).
• rm: número promedio de habitaciones por vivienda.
• age: proporción de unidades ocupadas por sus propietarios construidas antes de 1940.
• dis: media ponderada de las distancias a cinco centros de empleo de Boston.

• rad: índice de accesibilidad a carreteras radiales.
• tax: tasa de impuesto a la propiedad a valor completo por $ 10,000.
• ptratio: relación alumnos por profesor por ciudad.
• black: 1000 (Bk - 0.63) ^ 2 donde Bk es la proporción de negros por ciudad.+
• lstat: estado más bajo de la población (porcentaje).
• medv: valor medio de las viviendas ocupadas por sus propietarios en $ 1000.

realizar el mismo procedimiento para la selección de variables LASSO y comparar con cualquiera de los métodos anteriormente vistos en clase. (usar medv como variable respuesta y las demás variables como explicativas)

Cuando se quieren comparar dos modelos anidados se puede usar la teoría de verosimilitud.

• Supongamos que queremos comparar dos modelos \( M_0 \) y \( M_1 \) donde \( M_0 \subset M_1 \) , es decir, los parámetros en el modelo \( M_0 \) se obtienen imponiendo restricciones a los parámetros de \( M_1 \) .

• En otras palabras, el modelo \( M_0 \) es un caso especial del modelo \( M_1 \) . Por ejemplo, supongamos que tenemos los modelos lineales

\[ Y_i = \beta_0 + \beta_1X_1 + \beta_2X_2 +\varepsilon_i \qquad (M_1) \]

y

\[ Y_i = \beta_0 + \beta_1X_1 +\varepsilon_i \qquad (M_0) \]

• Note que podemos obtener el modelo \( M_0 \) del modelo \( M_1 \) al agregar la restricción \( \beta_2 = 0 \).

• En otras palabras, podemos probar cuál de los dos modelos es preferido usando los datos simplemente probando la hipótesis nula que

\[ H_0 : \beta_2 = 0vs. H_1 : \beta_2 \neq 0 \]

Este tipo de hipótesis se pueden probar usando la razón de verosimilitud. Suponga que \( l_0 \) y \( l_1 \) es el logaritmo de la verosimilitud bajo \( H_0 \) y \( H_1 \) , respectivamente.

• Bajo ciertas condiciones de regularidad(existencia de derivadas parciales) tenemos que:

\[ \Lambda_{LRT} = -2[l_0(\hat{\boldsymbol\theta}_0;datos) - l_1(\hat{\boldsymbol\theta}_1;datos)] \stackrel{asint.}{\sim} \chi_{p_1-p_2}^2, \]

donde \( p_0 \) y \( p_1 \) es el número de parámetros en \( \theta_0 \) y \( \theta_0 \), respectivamente.

• Note que la distribución de \( \Lambda \) LRT es aproximada o asintótica (es decir, cuando \( n \rightarrow \infty \)).

• En el caso de los modelos lineales vistos anteriormente, \( \theta \) incluiría tanto \( \beta \) como \( \sigma^2 \) .

• Lo anterior significa que si rechazamos \( H_0 \) (modelo reducido) entonces la evidencia sugiere que el modelo completo en \( H_1 \) es preferido. De lo contrario, la evidencia sugiere que el modelo reducido es preferido.

El uso del LRT debe ser cuidadoso especialmente cuando se ajustan modelos más complejos tales como modelos mixtos. (nose verán en clase)

Más adelante veremos algunos casos en los cuales la distribución asintótica del estadístico de prueba no se cumple (boundary space problem).

Finalmente, cuando queremos comparar dos modelos no anidados entonces recurrimos a los criterios de información. En el caso general estos criterios se definen como

\[ IC = -2L(\hat{\boldsymbol\theta};datos) + penalidad \]

Dependiendo de la penalidad que se use el criterio toma diferentes nombres. Asumamos que \( n \) es el número de individuos en la muestra y \( p \) es el número de parámetros en el modelo. Estos son los criterios comúnmente usados:

• \( AIC = -2L(\hat{\boldsymbol\theta};datos) + p \quad \) (Akaike’s Information Criterion)
• \( BIC = -2L(\hat{\boldsymbol\theta};datos) + p\ln(n) \quad \) (Bayesian Information Criterion)

La idea es seleccionar el modelo con el valor más pequeño de BIC o AIC.

• Las definiciones anteriores pueden variar dependiendo de la referencia pero en esencia van a diferir por una constante. Por eso es importante que tenga claro la fórmula que esta usando en cada problema o programa estadístico.