TP7 - Modelos de efectos mixtos

En este práctico vamos a trabajar con modelos lineales de efectos mixtos. Vamos a partir de los modelos lineales tradicionales y vamos a ver cómo surge la necesidad de incluir efectos aleatorios para analizar ciertos datasets.

Carga de bibliotecas

library(magrittr)
library(purrr)
library(tidyverse)

Antes de comenzar, vamos a cargar una de las bibliotecas que nos permite trabajar con modelos de efectos mixtos:

library(nlme)

Esta biblioteca contiene una serie de funciones que permiten resolver estos modelos, y además trae ciertos datasets interesantes para entender cómo funcionan los modelos y las funcionalidades de la biblioteca.

Ejercicio 1

Vamos a comenzar usando un dataset incorporado en la biblioteca nlme:

print(Rail)
#> Grouped Data: travel ~ 1 | Rail
#>    Rail travel
#> 1     1     55
#> 2     1     53
#> 3     1     54
#> 4     2     26
#> 5     2     37
#> 6     2     32
#> 7     3     78
#> 8     3     91
#> 9     3     85
#> 10    4     92
#> 11    4    100
#> 12    4     96
#> 13    5     49
#> 14    5     51
#> 15    5     50
#> 16    6     80
#> 17    6     85
#> 18    6     83

Este dataset contiene el tiempo que tarda una señal de ultrasonido en recorrer la longitud de un riel. Es un ejemplo muy sencillo tomado de la ingeniería (lo cuál muestra que este tipo de modelos puede ser útil en cualquier área de la ciencia y la tecnología). Muy importante!! No es necesario saber cómo ni para qué fin se obtuvo el dataset. Podemos pensar que tenemos una cierta medición (el tiempo de viaje) que fue obtenida muchas veces para distintas unidades de análisis (los rieles) tomadas aleatoriamente. El dataset nos muestra los valores del tiempo de viaje para cada riel.

Nota: Lo que vamos a ver a continuación sigue la presentación de Pinheiro y Bates en su libro ``Mixed-effect Models in S and S-PLUS´´ (capítulo 1).

  1. Grafique tiempo de viaje para cada riel. Interprete el gráfico. ¿Qué información podemos obtener?

  2. Ajuste un modelo lineal de efectos fijos:

\[ y_{ij} = \beta + \epsilon_{ij} \]

donde

  • \(y_{ij}\) es la observacion \(i\) para el riel \(j\) (en este caso, \(i = 1...M\) con \(M = 6\); y \(j = 1...n_i\) con \(n_i = 3\)). Tener en cuenta que cada riel podria tener una cantidad diferente de mediciones.
  • \(\beta\) es el tiempo medio para toda la poblacion de rieles (lo que queremos estimar).
  • \(\epsilon_{ij}\) es el término de error; representa realizaciones independiente a partir de una distribución Gaussiana con media 0 y desvío estándar \(\sigma\), es decir, \(N(0, \sigma)\). La variabilidad de este término es desconocida, por lo que también queremos estimar \(\sigma\).

Discuta los resultados obtenidos (\(\beta\) y \(\sigma\)). Recuerde que para calcular el modelo puede usar la función lm, y para ver los resultados, summary(). Recuerde también guardar el resultado del modelo en una variable del entorno R; esto permitirá procesar posteriormente las estimaciones del modelo.

  1. Obtenga los resiudos y grafíquelos. Para esto puede ser útil definir un nuevo tibble que contenga una copia de los datos originales y nos permita agregar una columna con los residuos (usando mutate).

    Discuta la información que le dan los residuos. ¿Están distribuídos en torno al 0? ¿Son normales? ¿Existe correlación entre los residuos provenientes de un mismo riel?

  2. Ajuste un modelo lineal de efectos fijos (mediante lm) que tenga en cuenta las diferencias entre cada riel:

\[y_{ij} = \beta_{i} + \epsilon_{ij}\]

(con \(i\) y \(j\) igual que antes), donde:

  • \(\beta_i\) es (la estimación de) el tiempo medio para el riel \(i\)-ésimo.
  • \(y_{ij}\) y \(\epsilon_{ij}\) están definidos igual que en el modelo anterior.

Discuta las estimaciones obtenidas (cada \(\beta_i\) y \(\sigma\)): ¿La estimación de la variabilidad es mayor o menor que en el caso anterior? ¿A qué se debe? ¿Qué representa cada \(\beta_i\)? ¿Cuántos grados de libertad tiene el modelo?

¿Es posible estimar parámetros de la población de rieles a partir de este modelo? ¿De qué manera?

  1. Obtenga los residuos, grafíquelos, y compare con el modelo previo.

  2. Analice los datos mediante un modelo de efectos mixtos, es decir, un modelo que considere explíticamente que los rieles fueron obtenidos aleatoriamente de una población de rieles, y que le permita estimar la variabilidad de dicha población. Para motivar dicho modelo, considere el modelo anterior reescrito de la siguiente manera:

\[ y_{ij} = \bar{\beta} + (\bar{\beta} - \beta_i) + \epsilon_{ij} \]

donde \(\hat{\beta}\) es el promedio de los \(\beta_i\). El término \(\bar{\beta} - \beta_i\) representa el desvío del riel \(i\)-ésimo respecto a (la estimación de la media de) la población de rieles \(\bar{\beta}\). Como los rieles fueron obtenidos aleatoriamente, esperamos que \(\bar{\beta} - \beta_i\) se comporte como una variable aleatoria Gaussiana centrada en 0 y con cierto desvío estándar que caracteriza la variabilidad de la población de rieles.

En base a esto, vamos a definir un nuevo modelo, reemplazando \(\hat{\beta}\) por el tiempo medio de viaje sobre la poblacion de rieles; y las desviaciones \(\beta_i - \bar{\beta}\) por una variable aleatoria que describa a la distribución de rieles (la cuál vamos a tratar de estimar).

El nuevo modelo será

\[y_{ij} = \beta + b_i + \epsilon_{ij}\]

donde:

  • \(\beta\) es (una estimación de) el tiempo medio sobre la poblacion de rieles.
  • \(b_i\) es (una realizacion de) una variable aleatoria que representa la variabilidad de cada riel presente en nuestra muestra respecto a la poblacion de rieles.
  • y \(\epsilon_{ij}\) es lo mismo de siempre.

Para terminar de definir el modelo, describimos a nuestras variables aleatorias:

  • La variabilidad en la población de rieles: \(b_i = N(0, \sigma_b)\).
  • La variabilidad en las mediciones para cada riel, dentro de la población: \(\epsilon_ij = N(0, \sigma)\).

Este es un modelo jerárquico, ya que su estructura tiene incluída la jerarquía que estaba presente en los datos.

Para esto, vamos a usar la función lme provista por la biblioteca nlme:

fm1Rail.lme <- lme(travel ~ 1, data = Rail, random = ~ 1 | Rail)

summary(fm1Rail.lme)
#> Linear mixed-effects model fit by REML
#>   Data: Rail 
#>       AIC      BIC   logLik
#>   128.177 130.6766 -61.0885
#> 
#> Random effects:
#>  Formula: ~1 | Rail
#>         (Intercept) Residual
#> StdDev:    24.80547 4.020779
#> 
#> Fixed effects:  travel ~ 1 
#>             Value Std.Error DF  t-value p-value
#> (Intercept)  66.5  10.17104 12 6.538173       0
#> 
#> Standardized Within-Group Residuals:
#>         Min          Q1         Med          Q3         Max 
#> -1.61882658 -0.28217671  0.03569328  0.21955784  1.61437744 
#> 
#> Number of Observations: 18
#> Number of Groups: 6

  1. Obtenga los parámetros \(\beta\), \(\sigma_b\) y \(\sigma\), y compare con los resultados anteriores.

  2. Obtenga los residuos del modelo considerando: i) únicamente los efectos fijos y ii) los efectos fijos más los aleatorios. Grafique y compare con los casos anteriores.

Para esto puede ser útil la función resid(), usando el argumento level (0 para los efectos fijos y 1 para fijos y aleatorios)

  1. Obtenga las predicciones del modelo considerando los efectos fijos, y los fijos más los aleatorios. Para esto, use la función predict() (con el argumento level igual que en el caso anterior)

Ejercicio 2

Vamos a simular un experimento en el cual se obtienen mediciones de una variable continua \(y\) en función de otra variable continua \(x\). Para obtener las mediciones, se utilizan ciertas unidades de análisis, que son tomadas aleatoriamente de una población. Para cada unidad de análisis, se obtienen mediciones de \(y\) en función de los valores de \(x\) seleccionados (medidas repetidas). La relación entre \(x\) e \(y\) es lineal, con pendientes y ordenadas propias de cada unidad de análisis muestreada.

Los resultados de esta simulación se guardan en un tibble:

n_uni = 5  # unidades de analisis

# parametros que caracterizan a la pendiente y la ordenada de la oblación de undiades de analisis
a=3
b=1
sigma_a= 1
sigma_b= 3

# variabilidad de y dentro de cada unidad de análisis
sigma = 0.5

# valores de la variable independiente x
x = 0:6

# muestreo a y b
pendientes = rnorm(n=n_uni, mean=a, sd=sigma_a)
ordenadas  = rnorm(n=n_uni, mean=b, sd=sigma_b)

# tibble con la info de las unidades de analisis (inaccesible in real life)
unidades = tibble(id_uni = 1:n_uni, a = pendientes, b=ordenadas, s=sigma)

# tibble con la simulación: y en función de x para cada u.d.a.
datos = unidades %>%
  uncount(length(x)) %>%
  mutate(x = rep(x, times=n_uni)) %>%
  mutate(y = x*a + b + rnorm(n=length(x), mean=0, sd=s)) %>%
  select(-a,-b,-s)

datos$id_uni %<>% as_factor()

  1. Grafique los datos. Discuta la jerarquía de los datos a partir del gráfico. ¿Qué parámetros de la generación de los datos afectan sus conclusiones?

  2. Aplique a los datos un modelo lineal de efectos fijos (usando lm). Obtenga los parámetros, residuos y predicciones. Grafique los residuos y las predicciones para cada unidad de análisis.

  3. Aplique a los datos un modelo lineal de efectos mixtos (usando lme). Discuta si corresponde aplicar solo ordenadas o pendientes aleatorias, o ambas. Para cada caso:

  • Obtenga los parámetros, residuos y predicciones.
  • Grafique los residuos y las predicciones para cada unidad de análisis.
m1.lme = lme(y ~ x, random = ~ x|id_uni, data=datos)
m2.lme = lme(y ~ x, random = ~ 1|id_uni, data=datos)
m3.lme = lme(y ~ x, random = ~ x-1|id_uni, data=datos)