10 Métodos de estimación

En este capítulo se muestran los métodos ML (maximum likelihood) y REML (restricted maximum likelihood) para estimar los parámetros de un modelo mixto.

10.1 ML y REML

El método de máxima verosimilitud restringida (o residual o reducida) (REML) es una alternativa de estimación de máxima verosimilitud que no basa las estimaciones en un ajuste de máxima verosimilitud de toda la información, sino que utiliza una función de verosimilitud calculada a partir de una conjunto de datos transformado, de modo que los parámetros molestos no tengan ningún efecto.

En el caso de la estimación del componente de varianza, el conjunto de datos original se reemplaza por un conjunto de contrastes calculados a partir de los datos, y la función de verosimilitud se calcula a partir de la distribución de probabilidad de estos contrastes, de acuerdo con el modelo para el conjunto de datos completo. En particular, REML se utiliza como método para ajustar modelos lineales mixtos. En contraste con la estimación de máxima verosimilitud anterior, REML puede producir estimaciones insesgadas de parámetros de varianza y covarianza.

Ejemplo: comparando REML y ML usando lme4

En este ejemplo vamos a simular datos de un modelo lineal mixto con intercepto y pendiente aleatoria, luego vamos a comparar las estimaciones de los parámetros del modelo usando REML y ML.

Solución.

Los datos los vamos a simular del siguiente modelo.

\[\begin{align*} y_{ij} | b_0, b_1 &\sim N(\mu_{ij}, \sigma^2_y) \\ \mu_{ij} &= 4 - 5 x_{ij} + b_{0i} + b_{1i} x_{ij} \\ \sigma^2_y &= 4 \\ \left ( \begin{matrix} b_{0} \\ b_{1} \end{matrix} \right ) &\sim N\left ( \left ( \begin{matrix} 0 \\ 0 \end{matrix} \right ), \left ( \begin{matrix} \sigma^2_{b0}=40 & \sigma_{b01}=3 \\ \sigma_{b01}=3 & \sigma^2_{b1}=50 \end{matrix} \right ) \right ) \\ x_{ij} &\sim U(0, 5) \end{align*}\]

Vamos a simular datos balanceados del anterior modelo considerando cinco grupos y diez observaciones por grupo. Para hacer esto vamos a usar la siguiente función.

gen_data <- function() {
  ni <- 10
  G  <- 5
  nobs <- ni * G # Numero total de observaciones
  grupo <- factor(rep(x=1:G, each=ni)) # Para crear la variable grupal
  obs <- rep(x=1:ni, each=G) # Para identificar las obs por grupo
  x <- runif(n=nobs, min=0, max=5) # La covariable
  library(MASS) # Libreria para simular obs de Normal bivariada
  Sigma <- matrix(c(40, 3, # Matriz de var-cov
                    3, 50), ncol=2, nrow=2)
  b <- mvrnorm(n=G, mu=rep(0, 2), Sigma=Sigma) # Simulando b0 y b1
  b <- apply(b, MARGIN=2, function(c) rep(c, each=ni)) # Replicando
  b0 <- as.vector(b[, 1]) # Separando los b0
  b1 <- as.vector(b[, 2]) # Separando los b1
  media <- 4 - 5 * x + b0 + b1 * x # La media
  y <- rnorm(n=nobs, mean=media, sd=2) # La variable respuesta
  datos <- data.frame(grupo, obs, x, y) # El dataframe
}

set.seed(123456)     # Fijando la semilla para replicar el ejemplo
datos <- gen_data()  # Creando los datos

Ahora vamos a ajustar dos modelos, uno con estimación REML y otro con estimación ML.

library(lme4)
fit_reml <- lmer(y ~ 1 + x + (1 + x| grupo), data = datos, REML = TRUE)
fit_ml   <- lmer(y ~ 1 + x + (1 + x| grupo), data = datos, REML = FALSE)

Para extraer sólo los efectos fijos de los dos modelos ajustados hacemos lo siguiente:

fixef(fit_reml)

## (Intercept)           x 
##   11.449474   -6.236896

fixef(fit_ml)

## (Intercept)           x 
##   11.449454   -6.235161

Los efectos fijos verdaderos son \(\beta_0=4\) y \(\beta_1=-5\). Al comparar las estimaciones REML y ML vemos que son cercanas entre sí y próximas de los parámetros.

Para extraer sólo las componentes de varianza de los dos modelos hacemos lo siguiente:

(vc_reml <- VarCorr(fit_reml))

##  Groups   Name        Std.Dev. Corr 
##  grupo    (Intercept) 6.5259        
##           x           4.0121   0.466
##  Residual             2.0379

(vc_ml   <- VarCorr(fit_ml))

##  Groups   Name        Std.Dev. Corr 
##  grupo    (Intercept) 5.8025        
##           x           3.5803   0.475
##  Residual             2.0377

Las componentes de varianza son \(\sigma_{b0}=6.32\) (obtenida como \(\sqrt{40}\)), \(\sigma_{b1}=7.07\) (obtenida como \(\sqrt{50}\)), \(\sigma_{b01}=3\) y \(\sigma_y=2\).

Vamos a calcular el MSE (mean squared error) para las dos estimaciones REML y ML.

a <- vc_reml[[1]]
b <- vc_ml[[1]]

mean((a[upper.tri(a, diag = TRUE)] - c(40, 3, 50))^2)

## [1] 413.6231

mean((b[upper.tri(b, diag = TRUE)] - c(40, 3, 50))^2)

## [1] 489.9286

De los anteriores resultados vemos que MSE es menor cuando se usa REML.

Ejemplo: comparando REML y ML usando nlme

Vamos a ajustar nuevamente los modelos pero usando el paquete nlme.

library(nlme)
fit_reml <- lme(y ~ 1 + x, random = ~ 1 + x | grupo, data = datos, method = "REML")
fit_ml   <- lme(y ~ 1 + x, random = ~ 1 + x | grupo, data = datos, method = "ML")

Para extraer sólo los efectos fijos de los dos modelos ajustados hacemos lo siguiente:

fixef(fit_reml)

## (Intercept)           x 
##   11.449474   -6.236896

fixef(fit_ml)

## (Intercept)           x 
##   11.449459   -6.235162

Los efectos fijos verdaderos son \(\beta_0=4\) y \(\beta_1=-5\). Al comparar las estimaciones REML y ML vemos que son cercanas entre sí, y próximas de los parámetros.

Para extraer sólo las componentes de varianza de los dos modelos hacemos lo siguiente:

(vc_reml <- VarCorr(fit_reml))

## grupo = pdLogChol(1 + x) 
##             Variance  StdDev   Corr  
## (Intercept) 42.588627 6.525996 (Intr)
## x           16.097479 4.012166 0.466 
## Residual     4.153191 2.037938

(vc_ml   <- VarCorr(fit_ml))

## grupo = pdLogChol(1 + x) 
##             Variance  StdDev   Corr  
## (Intercept) 33.668543 5.802460 (Intr)
## x           12.819501 3.580433 0.475 
## Residual     4.152341 2.037729

Las componentes de varianza son \(\sigma_{b0}=6.32\) (obtenida como \(\sqrt{40}\)), \(\sigma_{b1}=7.07\) (obtenida como \(\sqrt{50}\)), \(\sigma_{b01}=3\) y \(\sigma_y=2\).

Cuando se tienen muchos grupos y muchas repeticiones por grupo las estimaciones con REML y ML son muy cercanas.

Ejercicios

Repita el primer ejemplo de este capítulo para diez grupos y veinte observaciones por grupo. ¿Cuál es el valor de MSE para REML y ML?
Repita el primer ejemplo de este capítulo para treinta grupos y cincuenta observaciones por grupo. ¿Cuál es el valor de MSE para REML y ML?
Repita el primer ejemplo sin fijar la semilla y haga 100 réplicas con las siguientes combinaciones. Calcule el promedio de los MSE obtenidos.

\(G\)	\(n_i\)	\(\overline{MSE}_{REML}\)	\(\overline{MSE}_{ML}\)
5	10
10	15
20	20
30	25