Modelos de Efectos Aleatorios para un sólo factor

Ejemplo 2

Consideremos ahora el Ejercicio 5.1 de Kuehl (2000) sobre un estudio de herencia con animales de carne. Cinco sementales (animales machos) fueron apareados cada uno con un grupo separado de madres (animales hembras). Se registraron los pesos al nacer de ocho terneros machos (de diferentes madres) en cada uno de los cinco grupos de padres. Queremos averiguar qué parte de la variación total se debe a la variación entre diferentes toros.

pesos <- c(61, 100,  56, 113,  99, 103,  75,  62,  ## padre 1
            75, 102,  95, 103,  98, 115,  98,  94,  ## padre 2
            58,  60,  60,  57,  57,  59,  54, 100,  ## padre 3
            57,  56,  67,  59,  58, 121, 101, 101,  ## padre 4
            59,  46, 120, 115, 115,  93, 105,  75)  ## padre 5

padre <- factor(rep(1:5, each = 8))

animales <- data.frame(padre, pesos)
str(animales) # ver la estructura del dataframe
## 'data.frame':    40 obs. of  2 variables:
##  $ padre: Factor w/ 5 levels "1","2","3","4",..: 1 1 1 1 1 1 1 1 2 2 ...
##  $ pesos: num  61 100 56 113 99 103 75 62 75 102 ...

Realizamos una inspección visual de los datos

library(ggplot2)
p <- ggplot(animales, aes(padre, pesos))
p + geom_point() + theme_classic()

p2 <- ggplot(animales, aes(padre, pesos))+
  geom_boxplot()
p2

A primera vista, parece que la variación entre diferentes toros no es muy grande.

Ahora ajustamos el modelo de efectos aleatorios con la funcion lmer función en el paquete lme4 ( el cual debes instalar) .Queremos tener un efecto aleatorio por toro. Esto se puede especificar con la notación (1|padre)en la fórmula del modelo (El “ 1” significa que solo queremos tener una llamada intercepción aleatoria (αi por padre).

#install.packages(lme4)
library(lme4)
## Loading required package: Matrix
fit1 <- lmer(pesos ~ (1 | padre), data = animales)
summary(fit1)
## Linear mixed model fit by REML ['lmerMod']
## Formula: pesos ~ (1 | padre)
##    Data: animales
## 
## REML criterion at convergence: 358.2
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.9593 -0.7459 -0.1581  0.8143  1.9421 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  padre    (Intercept) 116.7    10.81   
##  Residual             463.8    21.54   
## Number of obs: 40, groups:  padre, 5
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)   82.550      5.911   13.96

Del resumen podemos leer en la tabla etiquetada Random effects las varianzas padre (116.7) y residuos (463.8), la varianza de yij por lo tanto se estima como 116.7 + 463.8 = 580.5 Por lo tanto, sólo cerca del 20 % (116.7/580.5) de la varianza total del peso al nacer se debe al padre (esta es la correlación intraclase). Esto confirma nuestra inspección visual inicial . Es una buena práctica verificar si la estructura de agrupación se interpretó correctamente. En la fila Number of obs podemos leer que tenemos 40 observaciones en total y un factor de agrupación dado por padre el cual tiene 5 niveles.

Debajo Fixed effects encontramos la estimación 82.55. Es una estimación del peso esperado al nacer de un ternero macho de un toro seleccionado al azar (seleccionado al azar de la población total de todos los toros).

Se pueden obtener intervalos de confianza aproximados del 95% para los parámetros con la función confint.

confint(fit1, oldNames = FALSE)
## Computing profile confidence intervals ...
##                         2.5 %   97.5 %
## sd_(Intercept)|padre  0.00000 24.61580
## sigma                17.32943 27.76544
## (Intercept)          69.83802 95.26197

El correspondiente intervalo de confianza para la varianza (bastante impreciso). La fila etiquetada con (Intercept)contiene el intervalo de confianza correspondiente para la media de peso a nacer.

¿Qué pasaría si usáramos la aovfunción “ordinaria” aquí?

Primero obtendremos el ANOVA para un diseño completamente al azar

options(contrasts = c("contr.sum", "contr.poly"))
fit2 <- aov(pesos  ~ padre, data= animales)
confint(fit2)
##                 2.5 %   97.5 %
## (Intercept)  75.63725 89.46275
## padre1      -12.75051 14.90051
## padre2        1.12449 28.77551
## padre3      -33.25051 -5.59949
## padre4      -18.87551  8.77551

Vemos que este intervalo de confianza es más corto que el del modelo de efectos aleatorios. La razón es la interpretación diferente. El modelo de efectos aleatorios permite hacer inferencias sobre la población de todos los toros (de los cuales hemos visto cinco hasta ahora) mientras que el modelo de efectos fijos permite hacer inferencias sobre estos cinco toros específicos . Por lo tanto, enfrentamos un problema más difícil con el modelo de efectos aleatorios, es por eso que tenemos menos confianza en nuestra estimación, lo que da como resultado intervalos de confianza más amplios en comparación con el modelo de efectos fijos.

Por supuesto, también deberíamos comprobar los supuestos del modelo. Aquí, esto incluye además la normalidad de los efectos aleatorios (aunque esto es difícil de verificar con solo cinco observaciones). Obtenemos el diagrama de Tukey-Anscombe con la plot función.

plot(fit1)

Para obtener diagramas QQ de los efectos aleatorios y los residuos, primero debemos extraerlos y luego usarlos qqnormcomo de costumbre.

par(mfrow = c(1, 2)) # esto configura la ventana que muestra el grafico
qqnorm(ranef(fit1)$padre[,"(Intercept)"], 
       main = "Random effects")
qqnorm(resid(fit1), main = "Residuals")

Modelos de Efectos Aleatorios para dos factores

Una gran empresa seleccionó al azar a 5 empleados y 6 lotes de material de origen del proceso de producción. El material de cada lote se dividió en 15 piezas que se asignaron al azar a los diferentes empleados de modo que cada empleado construyera 3 especímenes de prueba de cada lote. La respuesta fue el puntaje de calidad correspondiente. El objetivo era cuantificar las diferentes fuentes de variación.

library(readxl)
Datos1 <- read_excel("efectos_aleatorios.xlsx", 
    sheet = "Ej3")
View(Datos1)

Podemos, por ejemplo, visualizar este conjunto de datos con un gráfico de interacción.

with(Datos1, interaction.plot(x.factor = lote, 
                               trace.factor = empleado, 
                               response = score))

Parece haber variación entre los diferentes empleados y entre los diferentes lotes. La interacción no parece ser muy pronunciada. Establezcamos un modelo para estos datos: Denotamos el puntaje de calidad observado de la ésima muestra de empleados usando material del lote por . Un modelo natural es entonces donde es el efecto aleatorio (principal) del empleado , es el efecto aleatorio (principal) del lote.

El efecto aleatorio (principal) de empleado es la variabilidad entre diferentes empleados y el efecto aleatorio (principal) de lote es la variabilidad entre diferentes lotes. El término de interacción aleatoria puede, por ejemplo, interpretarse como inconsistencias de calidad de los empleados en diferentes lotes.

Ajustemos tal modelo en R. Queremos tener un efecto aleatorio por empleado (= (1 | empleado)), un efecto aleatorio por lote (= (1 | lote)) y un efecto aleatorio por combinación de empleado y lote (= (1 | empleado:lote)).

fit4 <- lmer(score ~ (1 | empleado) + (1 | lote) + 
                      (1 | empleado:lote), data = Datos1)
summary(fit4)
## Linear mixed model fit by REML ['lmerMod']
## Formula: score ~ (1 | empleado) + (1 | lote) + (1 | empleado:lote)
##    Data: Datos1
## 
## REML criterion at convergence: 167.5
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.8144 -0.5933 -0.1190  0.7073  2.8565 
## 
## Random effects:
##  Groups        Name        Variance Std.Dev.
##  empleado:lote (Intercept) 0.02349  0.1533  
##  lote          (Intercept) 0.51764  0.7195  
##  empleado      (Intercept) 1.54473  1.2429  
##  Residual                  0.22655  0.4760  
## Number of obs: 90, groups:  empleado:lote, 30; lote, 6; empleado, 5
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)  27.2478     0.6313   43.16

De la salida obtenemos las varianzas 1.544 (empleado), 0.517 (lote), 0.023 (interacción de empleado y lote) y 0.226 (término de error). Por lo tanto, la varianza total es 2,31

Vemos que la mayor contribución a la varianza es la variabilidad entre diferentes empleados que contribuye a aproximadamente 67% de la varianza total. Todas estas son estimaciones puntuales. Los intervalos de confianza en la escala de las desviaciones estándar se pueden obtener llamando a confint.

confint(fit4, oldNames = FALSE)
## Computing profile confidence intervals ...
##                                   2.5 %     97.5 %
## sd_(Intercept)|empleado:lote  0.0000000  0.3527999
## sd_(Intercept)|lote           0.4100058  1.4754427
## sd_(Intercept)|empleado       0.6694340  2.4993877
## sigma                         0.4020266  0.5741431
## (Intercept)                  25.9189910 28.5765641

La incertidumbre es bastante grande. Además, no hay evidencia estadística de que el término de interacción aleatoria sea realmente necesario, ya que el intervalo de confianza correspondiente contiene cero. Para obtener un intervalo de confianza más estrecho para la desviación estándar entre diferentes empleados, necesitaríamos muestrear a más empleados. Básicamente, cada empleado contribuye con una observación para estimar la desviación estándar (o varianza) entre empleados. El mismo razonamiento se aplica al lote.