przewidywanie () Funkcja dla mniejszych modeli efektów mieszanych

Problem:

Przeczytałem w innych postach, które predictnie są dostępne dla lmermodeli z efektami mieszanymi {lme4} w [R].

Próbowałem zgłębić ten temat za pomocą zestawu danych o zabawkach ...

Tło:

Zestaw danych jest dostosowany z tego źródła i dostępny jako ...

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

Oto pierwsze wiersze i nagłówki:

> 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

Mamy kilka powtarzających się obserwacji ( Time) ciągłego pomiaru, a mianowicie Recallczęstość niektórych słów, i kilka zmiennych objaśniających, w tym efekty losowe ( Auditoriumgdzie miał miejsce test; Subjectnazwa); i trwałe działanie , takie jak Education, Emotion(emocjonalne konotację słowa do zapamiętania) lub z spożywana przed badaniem.$\small \text{mgs.}$Caffeine

Chodzi o to, że łatwo jest zapamiętać hiperkofeinowe przewodowe obiekty, ale zdolność maleje z czasem, być może z powodu zmęczenia. Słowa o negatywnym znaczeniu są trudniejsze do zapamiętania. Edukacja ma przewidywalny efekt, a nawet audytorium odgrywa pewną rolę (być może jedna była bardziej głośna lub mniej wygodna). Oto kilka fabuł eksploracyjnych:

Różnice w częstości przywołania jako funkcję Emotional Tone, Auditoriumi Education:

Podczas dopasowywania linii w chmurze danych dla połączenia:

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

Dostaję tę fabułę:

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)

podczas gdy następujący model:

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

włączenie Timei równoległy kod daje zaskakujący wątek:

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

Pytanie:

Jak działa predictfunkcja w tym lmermodelu? Najwyraźniej bierze to pod uwagę Timezmienną, co skutkuje znacznie ściślejszym dopasowaniem oraz zygzakiem, który próbuje wyświetlić ten trzeci wymiar Timezobrazowany na pierwszym wykresie.

Jeśli zadzwonię predict(fit2), dostanę 132.45609pierwszy wpis, który odpowiada pierwszemu punktowi. Oto headzbiór danych z wyjściem predict(fit2)dołączonym jako ostatnia kolumna:

> 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

Współczynniki dla fit2:

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

Mój najlepszy zakład to ...

> 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

Jaką formułę wybrać zamiast tego 132.45609?

EDYCJA dla szybkiego dostępu ... Wzór do obliczenia przewidywanej wartości (zgodnie z przyjętą odpowiedzią będzie oparty na danych ranef(fit2)wyjściowych:

> 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

... dla pierwszego punktu wejścia:

> 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

Kod tego posta jest tutaj .

Łatwo się pomylić prezentacją współczynników podczas rozmowy coef(fit2). Spójrz na podsumowanie 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

Istnieje ogólny punkt przechwytywania wynoszący 61,92 dla modelu, ze współczynnikiem kofeiny wynoszącym 0,212. Tak więc dla kofeiny = 95 przewidujesz średnie wycofanie 82,06.

Zamiast używać coef, użyj, ranefaby uzyskać różnicę dla każdego przechwytywania z efektem losowym od średniego przechwytywania na wyższym poziomie zagnieżdżenia:

> 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

Wartości dla Jima w czasie = 0 będą różnić się od tej średniej wartości 82,06 sumą zarówno jego, jak Subject i jego Time:Subjectwspółczynników:

który moim zdaniem mieści się w błędzie zaokrąglenia wynoszącym 132,46.

Wartości zwracane przez przechwytują coefzdają się reprezentować ogólną osią oraz dodatkowo Subjectlub Time:Subjectspecyficzne różnice, więc jest to trudniejsze do pracy z osobami; jeśli spróbujesz wykonać powyższe obliczenia z coefwartościami, podwójnie policzymy ogólny punkt przecięcia.