Fonction Predict () pour les modèles d'effets mixtes lmer

Le problème:

J'ai lu dans d'autres articles qui predictne sont pas disponibles pour les lmermodèles d' effets mixtes {lme4} dans [R].

J'ai essayé d'explorer ce sujet avec un jeu de données jouet ...

Contexte:

L'ensemble de données est adapté de cette source et disponible en ...

require(gsheet)
data <- read.csv(text = 
     gsheet2text('https://docs.google.com/spreadsheets/d/1QgtDcGJebyfW7TJsB8n6rAmsyAnlz1xkT3RuPFICTdk/edit?usp=sharing',
        format ='csv'))

Ce sont les premières lignes et en-têtes:

> 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.56

Nous avons quelques observations répétées ( Time) d'une mesure continue, à savoir le Recalltaux de certains mots, et plusieurs variables explicatives, y compris des effets aléatoires ( Auditoriumoù le test a eu lieu; Subjectnom); et les effets fixes , tels que Education, Emotion(la connotation émotionnelle du mot à retenir), ou de ingérée avant l'essai.$\small \text{mgs.}$Caffeine

L'idée est qu'il est facile de se souvenir pour les sujets câblés hyper-caféinés, mais la capacité diminue avec le temps, peut-être en raison de la fatigue. Les mots à connotation négative sont plus difficiles à retenir. L'éducation a un effet prévisible, et même l'auditorium joue un rôle (peut-être que l'un était plus bruyant ou moins confortable). Voici quelques parcelles exploratoires:

Les différences dans les taux de rappel en fonction de Emotional Tone, Auditoriumet Education:

Lors de l'ajustement des lignes sur le nuage de données pour l'appel:

fit1 <- lmer(Recall ~ (1|Subject) + Caffeine, data = data)

Je reçois ce complot:

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)

tandis que le modèle suivant:

fit2 <- lmer(Recall ~ (1|Subject/Time) + Caffeine, data = data)

incorporer Timeet un code parallèle obtient une intrigue surprenante:

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 question:

Comment predictfonctionne la fonction dans ce lmermodèle? Évidemment, il prend en compte la Timevariable, ce qui entraîne un ajustement beaucoup plus serré et le zigzag qui essaie d'afficher cette troisième dimension de Timedépeint dans le premier tracé.

Si j'appelle, predict(fit2)j'obtiens 132.45609la première entrée, qui correspond au premier point. Voici headle jeu de données avec la sortie de predict(fit2)attaché comme dernière colonne:

> 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.39385

Les coefficients pour fit2sont:

$`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.2116271

Mon meilleur pari était ...

> 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

Quelle est la formule à adopter à la place 132.45609?

EDIT pour un accès rapide ... La formule pour calculer la valeur prédite (selon la réponse acceptée serait basée sur la ranef(fit2)sortie:

> 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

... pour le premier point d'entrée:

> 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.4561

Le code de cet article est ici .

Il est facile de se perdre dans la présentation des coefficients lorsque vous appelez coef(fit2). Regardez le résumé de fit2:

> summary(fit2)
Linear mixed model fit by REML ['lmerMod']
Formula: Recall ~ (1 | Subject/Time) + Caffeine
   Data: data

REML criterion at convergence: 444.5

Scaled residuals: 
 Min       1Q   Median       3Q      Max 
-1.88657 -0.46382 -0.06054  0.31430  2.16244 

Random effects:
 Groups       Name        Variance Std.Dev.
 Time:Subject (Intercept)  558.4   23.63   
 Subject      (Intercept) 2458.0   49.58   
 Residual                  675.0   25.98   
Number of obs: 45, groups:  Time:Subject, 15; Subject, 5

Fixed effects:
Estimate Std. Error t value
(Intercept) 61.91827   25.04930   2.472
Caffeine     0.21163    0.07439   2.845

Correlation of Fixed Effects:
 (Intr)
Caffeine -0.365

Il y a une interception globale de 61,92 pour le modèle, avec un coefficient de caféine de 0,212. Donc, pour la caféine = 95, vous prédisez un rappel moyen de 82,06.

Au lieu d'utiliser coef, utilisez ranefpour obtenir la différence de chaque interception à effet aléatoire à partir de l'interception moyenne au niveau d'imbrication immédiatement supérieur:

> 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

Les valeurs de Jim à Time = 0 différeront de cette valeur moyenne de 82,06 par la somme de son Subject et de ses Time:Subjectcoefficients:

qui, je pense, se situe dans l'erreur d'arrondi de 132,46.

Les valeurs d'interception renvoyées par coefsemblent représenter l'interception globale plus la Subjectou Time:Subjectdes différences spécifiques, il est donc plus difficile de travailler avec celles-ci; si vous essayez de faire le calcul ci-dessus avec les coefvaleurs, vous compteriez deux fois l'ordonnée à l'origine globale.