27 Coeficiente \(R^2\)
En este capítulo vamos mostrar como se puede calcular la medida \(R^2\) propuesta por Nakagawa, Schielzeth, and O’Hara (2013) y Nakagawa, Johnson, and Schielzeth (2017) para modelos lineales generalizados mixtos.
La función r.squaredGLMM() del paquete MuMIn de (R-MuMIn?) permite calcular \(R^2\) para varios tipos de modelos con efectos mixtos. La función puede calcular \(R^2\) para modelos de la clase merMod, lme, glmmTMB, glmmADMB, glmmPQL, cpglm(m) y
(g)lm.
Nakagawa, Johnson, and Schielzeth (2017) define dos \(R^2\) para LMM (linear mixed models) de la siguiente manera:
- \(R^2_{m}\): representa la varianza explicada por los efectos fijos.
- \(R^2_{c}\): representa la varianza explicada por todo el modelo, incluyendo ambos efectos fijos y aleatorios.
Las expresiones para obtener ambos \(R^2\) son:
\[ R^2_{m} = \frac{\sigma^2_f}{\sigma^2_f + \sigma^2_\alpha + \sigma^2_\epsilon}, \]
\[ R^2_{c} = \frac{\sigma^2_f + \sigma^2_\alpha}{\sigma^2_f + \sigma^2_\alpha + \sigma^2_\epsilon}, \]
donde \(\sigma^2_f\) es la varianza asociada a los efectos fijos, \(\sigma^2_\alpha\) es la varianza de los efectos aleatorios y \(\sigma^2_\epsilon\) es la varianza asociada a los valores de la variable respuesta.
Nakagawa, Johnson, and Schielzeth (2017) en su artículo proponen el cálculo de \(R^2\) para otros GLMM (generalized linear mixed models). Los autores afirman que el reto al proponer un \(R^2\) para los GLMM radica en la definición de \(\sigma^2_\epsilon\). En la implementación de \(R^2\) para GLMM hecha por (R-MuMIn?) se pueden obtener tres formas para cuantificar \(\sigma^2_\epsilon\), esas formas son: “delta”, “lognormal” y “trigamma”.
En los siguientes ejemplos vamos a ilustrar el uso de la función r.squaredGLMM().
Ejemplo: modelo normal con intercepto aleatorio
En este ejemplo vamos a simular observaciones \(n_i=50\) observaciones para \(G=10\) grupos (en total 500 obs) que tengan la estructura mostrada abajo y en la cual la varianza de los efectos aleatorio \(b_0\) (\(\sigma^2_{b0}=625\)) es más grande que la varianza de las observaciones (\(\sigma^2_y=16\)).
El objetivo del ejemplo es calcular el \(R^2\) para el modelo ajustado correctamente.
El modelo para simular los datos es el siguiente.
\[\begin{align*} y_{ij} | b_0 &\sim N(\mu_{ij}, \sigma^2_y) \\ \mu_{ij} &= 4 - 6 x_{ij} + b_{0i} \\ \sigma^2_y &= 16 \\ b_{0} &\sim N(0, \sigma^2_{b0}=625) \\ x_{ij} &\sim U(-5, 6) \end{align*}\]
El vector de parámetros de este modelo es \(\boldsymbol{\Theta}=(\beta_0=4, \beta_1=-6, \sigma_y=4, \sigma_{b0}=25)^\top\).
El código para simular las 500 observaciones se muestra a continuación.
ni <- 50
G <- 10
nobs <- ni * G
grupo <- factor(rep(x=1:G, each=ni))
obs <- rep(x=1:ni, times=G)
x <- runif(n=nobs, min=-5, max=6)
b0 <- rnorm(n=G, mean=0, sd=sqrt(625)) # Intercepto aleatorio
b0 <- rep(x=b0, each=ni) # El mismo intercepto aleatorio pero repetido
media <- 4 - 6 * x + b0
y <- rnorm(n=nobs, mean=media, sd=sqrt(16))
datos <- data.frame(obs, grupo, b0, x, media, y)Ahora vamos a ajustar el modelo correcto así:
A continuación vamos a calcular los \(R^2\).
##
## Attaching package: 'MuMIn'## The following object is masked from 'package:randomForest':
##
## importance## Warning: 'r.squaredGLMM' now calculates a revised statistic. See the help page.## R2m R2c
## [1,] 0.4287911 0.9826349En la salida anterior se reporta la medida marginal y la medida condicional del \(R^2\). De los resultados se observa que el \(R^2_{m}\) es mayor que el \(R^2_{c}\), eso significa que el modelo completo (con efectos fijos y aleatorio) logra explicar un 98% de la varianza. En otras palabras, significa que si valió la pena ajustar el modelo con esos efectos aleatorios. Era de esperarse que esto sucediera porque los datos se simularon con \(\sigma^2_{b0}=625\).
Ejemplo: modelo normal sin intercepto aleatorio
En este ejemplo vamos a repetir a simular datos como el ejemplo anterior pero con \(\sigma^2_{b0}=0\). La idea es ver si el \(R^2\) logra detectar que no se necesita incluir un intercepto aleatorio al ajustar el modelo.
El código para simular las 500 observaciones se muestra a continuación.
ni <- 50
G <- 10
nobs <- ni * G
grupo <- factor(rep(x=1:G, each=ni))
obs <- rep(x=1:ni, times=G)
x <- runif(n=nobs, min=-5, max=6)
b0 <- rnorm(n=G, mean=0, sd=0) # Intercepto aleatorio
b0 <- rep(x=b0, each=ni) # El mismo intercepto aleatorio pero repetido
media <- 4 - 6 * x + b0
y <- rnorm(n=nobs, mean=media, sd=sqrt(16))
datos <- data.frame(obs, grupo, b0, x, media, y)Ahora vamos a ajustar el modelo así:
A continuación vamos a calcular los \(R^2\).
## R2m R2c
## [1,] 0.9571738 0.9578291De la salida anterior vemos que ambos \(R^2\) son muy parecidos. Esto significa que el porcentaje de la varianza explicada por un modelo con ambos tipos de efectos (fijos y aleatorios) es igual a la varianza explicada por un modelo sólo con efectos fijos, en otras palabras, esto indica que NO vale la pena incluir efectos aleatorios en el modelo. Este resultado era el esperado porque los datos fueron simulados con \(\sigma^2_{b0}=0\).
Ejemplo: modelo clase lme
En este ejemplo vamos a ajustar un modelo para explicar la desde la hipófisis hasta la fisura pterigomaxilar (mm), en función del sexo y la edad de los pacientes (Subject). La base de datos a usar es Orthodont. Luego se desea calcular el \(R^2\) para el modelo ajustado.
El código que vamos a usar es el siguiente.
library(nlme)
library(MuMIn)
data(Orthodont, package = "nlme")
fm1 <- lme(distance ~ Sex * age, ~ 1 | Subject, data = Orthodont)
r.squaredGLMM(fm1)## R2m R2c
## [1,] 0.4098416 0.7827266De la salida anterior se observa que 78.27% de la variabilidad observada se logra explicar con el modelo ajustado con efectos fijos y un intercepto aleatorio.
Ejemplo: modelo clase lme4
En este ejemplo vamos a calcular \(R^2\) para el modelo gamma presentado en el ejemplo 1 sobre conductores que está en el capítulo 24.
A continuación se presenta el mismo código usado en el ejemplo 1 del capítulo 24. En la última línea de este código se calcula el \(R^2\).
library(hglm)
data(semiconductor)
library(lme4)
mod1 <- glmer(y ~ x1 + x2 + x3 + x4 + x5 + x6 + (1 | Device),
data = semiconductor,
family = Gamma(link = log))
library(MuMIn)
r.squaredGLMM(mod1)## R2m R2c
## delta 0.7112290 0.7749512
## lognormal 0.7185672 0.7829469
## trigamma 0.7033419 0.7663575De la salida anterior se observa que el \(R^2_c\) es ligeramente mayor al \(R^2_m\). Los valores de \(R^2_c\) son alrededor de 77%, esto indica que el modelo mod1 con efectos fijos y un intercepto aletorio, logra explicar un 77% de la variabilidad observada. Este resultado es una evidencia en favor del modelo mixto frente al modelo sólo con efectos fijos.