El problema:
He leído en otros puestos que predict no está disponible para efectos mixtos lmer {lme4} en [R].
He intentado explorar este tema con un conjunto de datos de juguete...
Antecedentes:
El conjunto de datos está adaptado de esta fuente y disponible como...
require(gsheet)
data <- read.csv(text =
gsheet2text('https://docs.google.com/spreadsheets/d/1QgtDcGJebyfW7TJsB8n6rAmsyAnlz1xkT3RuPFICTdk/edit?usp=sharing',
format ='csv'))Estas son las primeras filas y cabeceras:
> head(data)
Subject Auditorium Education Time Emotion Caffeine Recall
1 Jim A HS 0 Negative 95 125.80
2 Jim A HS 0 Neutral 86 123.60
3 Jim A HS 0 Positive 180 204.00
4 Jim A HS 1 Negative 200 95.72
5 Jim A HS 1 Neutral 40 75.80
6 Jim A HS 1 Positive 30 84.56Tenemos algunas observaciones repetidas ( Time ) de una medida continua, es decir, el Recall de algunas palabras, y varias variables explicativas, incluyendo efectos aleatorios ( Auditorium donde se realizó la prueba; Subject nombre); y efectos fijos como por ejemplo Education , Emotion (la connotación emocional de la palabra para recordar), o $\small \text{mgs.}$ de Caffeine ingerido antes de la prueba.
La idea es que es fácil de recordar para los sujetos cableados e hipercafeinados, pero la capacidad disminuye con el tiempo, quizás debido al cansancio. Las palabras con connotación negativa son más difíciles de recordar. La educación tiene un efecto predecible, e incluso el auditorio influye (quizás uno era más ruidoso, o menos cómodo). He aquí un par de gráficos exploratorios:
Diferencias en la tasa de recuerdo en función de Emotional Tone , Auditorium y Education :
Al ajustar las líneas en la nube de datos para la llamada:
fit1 <- lmer(Recall ~ (1|Subject) + Caffeine, data = data)
Me sale esta trama:
library(ggplot2)
p <- ggplot(data, aes(x = Caffeine, y = Recall, colour = Subject)) +
geom_point(size=3) +
geom_line(aes(y = predict(fit1)),size=1)
print(p)mientras que el siguiente modelo:
fit2 <- lmer(Recall ~ (1|Subject/Time) + Caffeine, data = data)
incorporando Time y un código paralelo obtiene una trama sorprendente:
p <- ggplot(data, aes(x = Caffeine, y = Recall, colour = Subject)) +
geom_point(size=3) +
geom_line(aes(y = predict(fit2)),size=1)
print(p)
La pregunta:
¿Cómo es que el predict operan en este lmer ¿Modelo? Evidentemente, está teniendo en cuenta la Time variable, lo que da lugar a un ajuste mucho más estrecho, y el zig-zag que trata de mostrar esta tercera dimensión de Time retratado en la primera trama.
Si llamo a predict(fit2) Me sale 132.45609 para la primera entrada, que corresponde al primer punto. Aquí está el head del conjunto de datos con la salida de predict(fit2) adjunta como última columna:
> data$predict = predict(fit2)
> head(data)
Subject Auditorium Education Time Emotion Caffeine Recall predict
1 Jim A HS 0 Negative 95 125.80 132.45609
2 Jim A HS 0 Neutral 86 123.60 130.55145
3 Jim A HS 0 Positive 180 204.00 150.44439
4 Jim A HS 1 Negative 200 95.72 112.37045
5 Jim A HS 1 Neutral 40 75.80 78.51012
6 Jim A HS 1 Positive 30 84.56 76.39385Los coeficientes de fit2 son:
$`Time:Subject`
(Intercept) Caffeine
0:Jason 75.03040 0.2116271
0:Jim 94.96442 0.2116271
0:Ron 58.72037 0.2116271
0:Tina 70.81225 0.2116271
0:Victor 86.31101 0.2116271
1:Jason 59.85016 0.2116271
1:Jim 52.65793 0.2116271
1:Ron 57.48987 0.2116271
1:Tina 68.43393 0.2116271
1:Victor 79.18386 0.2116271
2:Jason 43.71483 0.2116271
2:Jim 42.08250 0.2116271
2:Ron 58.44521 0.2116271
2:Tina 44.73748 0.2116271
2:Victor 36.33979 0.2116271
$Subject
(Intercept) Caffeine
Jason 30.40435 0.2116271
Jim 79.30537 0.2116271
Ron 13.06175 0.2116271
Tina 54.12216 0.2116271
Victor 132.69770 0.2116271Mi mejor apuesta fue...
> coef(fit2)[[1]][2,1]
[1] 94.96442
> coef(fit2)[[2]][2,1]
[1] 79.30537
> coef(fit2)[[1]][2,2]
[1] 0.2116271
> data$Caffeine[1]
[1] 95
> coef(fit2)[[1]][2,1] + coef(fit2)[[2]][2,1] + coef(fit2)[[1]][2,2] * data$Caffeine[1]
[1] 194.3744¿Cuál es la fórmula para llegar en su lugar a 132.45609 ?
EDITAR para acceder rápidamente... La fórmula para calcular el valor previsto (según la respuesta aceptada se basaría en el ranef(fit2) de salida:
> ranef(fit2)
$`Time:Subject`
(Intercept)
0:Jason 13.112130
0:Jim 33.046151
0:Ron -3.197895
0:Tina 8.893985
0:Victor 24.392738
1:Jason -2.068105
1:Jim -9.260334
1:Ron -4.428399
1:Tina 6.515667
1:Victor 17.265589
2:Jason -18.203436
2:Jim -19.835771
2:Ron -3.473053
2:Tina -17.180791
2:Victor -25.578477
$Subject
(Intercept)
Jason -31.513915
Jim 17.387103
Ron -48.856516
Tina -7.796104
Victor 70.779432... para el primer punto de entrada:
> summary(fit2)$coef[1]
[1] 61.91827 # Overall intercept for Fixed Effects
> ranef(fit2)[[1]][2,]
[1] 33.04615 # Time:Subject random intercept for Jim
> ranef(fit2)[[2]][2,]
[1] 17.3871 # Subject random intercept for Jim
> summary(fit2)$coef[2]
[1] 0.2116271 # Fixed effect slope
> data$Caffeine[1]
[1] 95 # Value of caffeine
summary(fit2)$coef[1] + ranef(fit2)[[1]][2,] + ranef(fit2)[[2]][2,] +
summary(fit2)$coef[2] * data$Caffeine[1]
[1] 132.4561El código de este puesto es aquí .
