predict () -Funktion für ältere Mixed-Effects-Modelle

Das Problem:

Ich habe in anderen Posts gelesen , predictdie für lmerModelle mit gemischten Effekten {lme4} in [R] nicht verfügbar sind .

Ich habe versucht, dieses Thema mit einem Spielzeugdatensatz zu untersuchen ...

Hintergrund:

Der Datensatz ist aus dieser Quelle angepasst und als ...

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

Dies sind die ersten Zeilen und Überschriften:

> 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

Wir haben einige wiederholte Beobachtungen ( Time) einer kontinuierlichen Messung, nämlich die RecallRate einiger Wörter, und mehrere erklärende Variablen, einschließlich zufälliger Effekte ( Auditoriumwo der Test stattgefunden hat; SubjectName); und festen Effekten , wie zum Beispiel Education, Emotion(die emotionale Konnotation des Wortes zu erinnern) oder $\small \text{mgs.}$die Caffeinevor dem Test aufgenommen.

Die Idee ist, dass es leicht ist, sich an hyperkoffeinhaltige verdrahtete Probanden zu erinnern, aber die Fähigkeit nimmt mit der Zeit ab, möglicherweise aufgrund von Müdigkeit. Wörter mit negativer Konnotation sind schwerer zu merken. Bildung hat einen vorhersehbaren Effekt, und sogar das Auditorium spielt eine Rolle (vielleicht war eines lauter oder unangenehmer). Hier sind ein paar Versuchspläne:

Unterschiede in der Recall - Rate in Abhängigkeit von Emotional Tone, Auditoriumund Education:

Beim Anpassen von Leitungen in der Datenwolke für den Anruf:

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

Ich bekomme diese Handlung:

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)

während das folgende Modell:

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

Wenn Sie Timeeinen Parallelcode einbauen, erhalten Sie eine überraschende Darstellung:

p <- ggplot(data, aes(x = Caffeine, y = Recall, colour = Subject)) +
  geom_point(size=3) +
  geom_line(aes(y = predict(fit2)),size=1) 
print(p)

Die Frage:

Wie funktioniert die predictFunktion in diesem lmerModell? Offensichtlich berücksichtigt es die TimeVariable, was zu einer viel engeren Anpassung führt, und das Zick-Zack-Verhältnis, das versucht, diese dritte Dimension der TimeDarstellung in der ersten Darstellung anzuzeigen .

Wenn ich anrufe predict(fit2)bekomme ich 132.45609für den ersten Eintrag, der dem ersten Punkt entspricht. Hier ist der headDatensatz mit der Ausgabe von predict(fit2)als letzte Spalte angehängt:

> 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

Die Koeffizienten für fit2sind:

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

Meine beste Wette war ...

> 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

Was ist die Formel, um stattdessen zu bekommen 132.45609?

EDIT für den schnellen Zugriff ... Die Formel zur Berechnung des vorhergesagten Wertes (entsprechend der akzeptierten Antwort würde auf der ranef(fit2)Ausgabe basieren :

> 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

... für den ersten Einstiegspunkt:

> 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

Der Code für diesen Beitrag ist hier .

Die Darstellung von Koeffizienten beim Aufrufen kann leicht zu Verwirrung führen coef(fit2). Schauen Sie sich die Zusammenfassung von fit2 an:

> 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

Es gibt einen Gesamtabschnitt von 61,92 für das Modell mit einem Koffeinkoeffizienten von 0,212. Für Koffein = 95 prognostizieren Sie also einen durchschnittlichen Rückruf von 82,06.

Anstatt zu verwenden coef, verwenden Sie ranef, um die Differenz der einzelnen Zufallseffektabschnitte von den mittleren Abschnitten auf der nächsthöheren Verschachtelungsebene zu ermitteln:

> 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

Die Werte für Jim zum Zeitpunkt = 0 unterscheiden sich von diesem Durchschnittswert von 82,06 um die Summe seiner Subject und seiner Time:SubjectKoeffizienten:

was meiner Meinung nach innerhalb des Rundungsfehlers von 132,46 liegt.

Die von zurückgegebenen Intercept-Werte coefscheinen den gesamten Intercept zuzüglich des zu repräsentierenSubject oder Time:Subjectspezifischer Unterschiede , daher ist es schwieriger, mit diesen zu arbeiten. Wenn Sie versuchen würden, die obige Berechnung mit den coefWerten durchzuführen, würden Sie den Gesamtabschnitt doppelt zählen.