¿REML o ML para comparar dos modelos de efectos mixtos con efectos fijos diferentes, pero con el mismo efecto aleatorio?

Question

Antecedentes: Nota: Mi conjunto de datos y el código R se incluyen a continuación del texto

Deseo utilizar el AIC para comparar dos modelos de efectos mixtos generados con el lme4 en R. Cada modelo tiene un efecto fijo y un efecto aleatorio. El efecto fijo difiere entre los modelos, pero el efecto aleatorio sigue siendo el mismo entre los modelos. He descubierto que si utilizo REML=TRUE , model2 tiene la puntuación AIC más baja, pero si utilizo REML=FALSE , model1 tiene la puntuación AIC más baja.

Soporte para el uso de ML:

Zuur et al. (2009; p. 122) sugieren que "Para comparar modelos con efectos fijos anidados (pero con la misma estructura aleatoria), debe utilizarse la estimación ML y no REML". Esto me indica que debería utilizar ML ya que mis efectos aleatorios son los mismos en ambos modelos, pero mis efectos fijos difieren. [Zuur et al. 2009. Mixed Effect Models and Extensions in Ecology with R. Springer].

Soporte para el uso de REML:

Sin embargo, observo que cuando utilizo ML, la varianza residual asociada a los efectos aleatorios difiere entre los dos modelos ( model1 = 136.3; model2 = 112,9), pero cuando utilizo REML, es el mismo entre los modelos (modelo1=modelo2=151,5). Esto me da a entender que debería utilizar REML para que la varianza residual aleatoria sea la misma entre los modelos con la misma variable aleatoria.

Pregunta:

¿No tiene más sentido utilizar REML que ML para las comparaciones de modelos en los que los efectos fijos cambian y los efectos aleatorios permanecen iguales? Si no es así, ¿puede explicar por qué o indicarme otra literatura que lo explique mejor?

# Model2 "wins" if REML=TRUE:
REMLmodel1 = lmer(Response ~ Fixed1 + (1|Random1),data,REML = TRUE)
REMLmodel2 = lmer(Response ~ Fixed2 + (1|Random1),data,REML = TRUE)
AIC(REMLmodel1,REMLmodel2)
summary(REMLmodel1)
summary(REMLmodel2)

# Model1 "wins" if REML=FALSE:
MLmodel1 = lmer(Response ~ Fixed1 + (1|Random1),data,REML = FALSE)
MLmodel2 = lmer(Response ~ Fixed2 + (1|Random1),data,REML = FALSE)
AIC(MLmodel1,MLmodel2)
summary(MLmodel1)
summary(MLmodel2)

Conjunto de datos:

Response    Fixed1  Fixed2  Random1
5.20    A   A   1
32.50   A   A   1
6.57    A   A   2
24.77   A   B   3
41.69   A   B   3
34.29   A   B   4
1.80    A   B   4
10.00   A   B   5
15.56   A   B   5
4.44    A   C   6
21.65   A   C   6
9.20    A   C   7
4.11    A   C   7
12.52   B   D   8
0.25    B   D   8
27.34   B   D   9
11.54   B   E   10
0.86    B   E   10
0.68    B   E   11
4.00    B   E   11

Answer 1

Zuur et al., y Faraway (del comentario de @janhove más arriba) tienen razón; el uso de métodos basados en la verosimilitud (incluido el AIC) para comparar dos modelos con diferentes efectos fijos que se ajustan mediante REML generalmente conducirá a un sinsentido.

Lejos (2006) Ampliación del modelo lineal con R (p. 156):

La razón es que REML estima los efectos aleatorios considerando combinaciones lineales de los datos que eliminan los efectos fijos. Si se modifican estos efectos fijos, las probabilidades de los dos modelos no serán directamente comparables

Estas dos preguntas profundizan en la cuestión: Permitió comparaciones de modelos de efectos mixtos (efectos aleatorios principalmente) ; REML vs ML stepAIC

Answer 2

Daré un ejemplo para ilustrar por qué la probabilidad REML no puede utilizarse para cosas como las comparaciones AIC. Imaginemos que tenemos un modelo normal de efectos mixtos. Dejemos que $X$ denota la matriz de diseño y supone que esta matriz tiene rango completo. Podemos encontrar una reparametrización del espacio de valores medios, dada por la matriz $\tilde{X}$ . Las dos matrices abarcan el mismo subespacio lineal de $\mathbb{R}^n$ . Así, las columnas de $\tilde{X}$ pueden escribirse como combinaciones lineales de las columnas de $X$ . Por lo tanto, podemos encontrar una matriz cuadrática, $B$ , de tal manera que

$\tilde{X} = XB$ .

Además, $B$ tiene rango completo (esto se puede demostrar suponiendo que no lo tuviera; entonces tampoco lo tendría $X$ , una contradicción). Esto significa que $B$ es invertible.

Si empezamos utilizando la segunda parametrización del espacio de valores medios y dejamos que $V$ sea una matriz de covarianza entonces consideremos el criterio REML que debemos maximizar (estoy omitiendo una constante)

$ |V|^{-1/2}|\tilde{X}'V^{-1}\tilde{X}|^{-1/2}\exp((y-\tilde{X}\tilde{\beta})'V^{-1}(y-\tilde{X}\tilde{\beta})/2) $ ,

sobre el conjunto de parámetros, donde $\beta = (\tilde{X}V^{-1}\tilde{X})^{-1}y$ . Utilizando el hecho de que $X = \tilde{X}B$ podemos darnos cuenta de que esto se puede reescribir como

$ |B||V|^{-1/2}||X'V^{-1}X|^{-1/2}|\exp((y-X\bar{\beta})'V^{-1}(y-X\bar{\beta})/2) $ ,

donde $\bar{\beta} = (XV^{-1}X)^{-1}y$ . Esta es la probabilidad REML para la otra parametrización veces el determinante de $|B|$ .

Por lo tanto, tenemos un ejemplo de dos parametrizaciones diferentes del mismo modelo, que dan diferentes valores de probabilidad, suponiendo que $|B| \neq 1$ (dicha matriz puede encontrarse fácilmente). El mismo valor del parámetro maximizará el criterio en ambos casos, pero el valor de la probabilidad será diferente. Esto muestra que hay un elemento arbitrario en el valor de la verosimilitud y, por lo tanto, ilustra por qué no se puede utilizar el valor de la verosimilitud para la comparación entre modelos con diferentes efectos fijos: se podrían cambiar los resultados simplemente cambiando la parametrización del espacio del valor medio en uno de los modelos.

Este es un ejemplo de por qué no se debe utilizar REML cuando se comparan modelos con diferentes efectos fijos. Sin embargo, REML suele estimar mejor los parámetros de los efectos aleatorios, por lo que a veces se recomienda utilizar ML para las comparaciones y REML para estimar un único modelo (quizás el definitivo).