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

27

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


entrez la description de l'image ici

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

entrez la description de l'image ici


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:

avec le code suivant (notez l'appel pour <code> predire (fit1) </code> dedans):

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)

entrez la description de l'image ici


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 .

Antoni Parellada
la source
3
Notez qu'il y a une predictfonction dans ce package depuis la version 1.0-0, publiée le 2013-08-01. Voir la page des nouvelles du package dans CRAN . S'il n'y en avait pas eu, vous n'auriez pas pu obtenir de résultats predict. N'oubliez pas que vous pouvez voir le code R avec lme4 ::: predire.merMod à l'invite de commande R et inspecter la source pour toutes les fonctions compilées sous-jacentes dans le package source pour lme4.
EdM
1
Merci, je vois qu'il y a la fonctionnalité avec les effets aléatoires ignorés ou inclus. Savez-vous où je peux trouver comment il est calculé? Si je tape ?predictsur la console [r], j'obtiens les prévisions de base pour {stats} ...
Antoni Parellada
@EdM ... oui, c'est nouveau pour moi ... Merci. Je n'avais pas appelé predict.merMod, cependant ... Comme vous pouvez le voir sur l'OP, j'ai simplement appelé predict...
Antoni Parellada
1
Chargez le lme4package, puis tapez lme4 ::: Predict.merMod pour voir la version spécifique au package. La sortie de lmerest stockée dans un objet de classe merMod.
EdM
4
L'une des beautés de R est qu'une fonction comme predictsait quoi faire en fonction de la classe de l'objet sur lequel elle est appelée à agir. Vous appeliez predict.merMod, vous ne le saviez tout simplement pas.
EdM

Réponses:

25

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:

82.06+17,39+33.04=132,49

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.

EdM
la source
Merci! C'était exceptionnel! Je ne pense pas qu'il soit inutile de le laisser ouvert ... C'est la réponse, n'est-ce pas?
Antoni Parellada
Je suis l'indication sur ranefde l' examen du code de R dans lme4. J'ai clarifié la présentation à quelques endroits.
EdM
(+1) Remarque: les effets aléatoires temporels imbriqués dans la personne semblent étranges.
Michael M
@MichaelM: Oui, les données présentées semblent être un plan croisé (temps x sujet) plutôt qu'un plan imbriqué, mais c'est ainsi que le PO a soulevé la question de savoir comment interpréter la lme4sortie. De plus, les données présentées semblaient être à titre d'illustration plutôt que d'être une véritable étude à analyser.
EdM