przewidywanie () Funkcja dla mniejszych modeli efektów mieszanych

27

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


wprowadź opis zdjęcia tutaj

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

wprowadź opis zdjęcia tutaj


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

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

Dostaję tę fabułę:

z następującym kodem (zwróć uwagę na wezwanie do <code> prognozy (fit1) </code>):

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)

wprowadź opis zdjęcia tutaj


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 .

Antoni Parellada
źródło
3
Zauważ, że predictw tym pakiecie jest funkcja od wersji 1.0-0, wydanej 01.01.2013. Zobacz stronę aktualności pakietu w CRAN . Gdyby tego nie było, nie uzyskałbyś żadnych rezultatów predict. Nie zapominaj, że możesz zobaczyć kod R za pomocą lme4 ::: Prognozuj.merMod w wierszu polecenia R i sprawdź źródło dla wszystkich podstawowych kompilowanych funkcji w pakiecie źródłowym lme4.
EdM
1
Dzięki, widzę, że funkcjonalność z losowymi efektami została pominięta lub uwzględniona. Czy wiesz, gdzie mogę znaleźć sposób obliczania? Jeśli piszę ?predictna konsoli [r], otrzymuję podstawową prognozę {statystyki} ...
Antoni Parellada,
@EdM ... tak, to dla mnie nowość ... Dziękuję. Nie zadzwoniłem predict.merModjednak ... Jak widać na OP, zadzwoniłem po prostu predict...
Antoni Parellada,
1
Załaduj lme4pakiet, a następnie wpisz lme4 ::: Prognozuj.merMod, aby zobaczyć wersję specyficzną dla pakietu. Dane wyjściowe lmersą przechowywane w obiekcie klasy merMod.
EdM
4
Jedną z zalet R jest to, że funkcja taka predictwie, co robić w zależności od klasy obiektu, na który ma działać. Dzwoniłeś predict.merMod, po prostu o tym nie wiedziałeś.
EdM

Odpowiedzi:

25

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

82,06+17,39+33.04=132,49

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.

EdM
źródło
Dziękuję Ci! To było znakomite! Nie sądzę, że istnieje jakikolwiek punkt w pozostawiając ją otworzyć ... To odpowiedź, prawda?
Antoni Parellada,
Dostałem wskazówkę o ranefod zbadaniu kodu R w lme4. Wyjaśniłem prezentację w kilku miejscach.
EdM
(+1) Uwaga: efekty losowe w czasie zagnieżdżone w osobie wyglądają dziwnie.
Michael M
@MichaelM: Tak, przedstawione dane wydają się być przekreślonym (Time x Subject), a nie zagnieżdżonym projektem, ale w ten sposób OP podniósł kwestię interpretacji lme4wyników. Przedstawione dane wydawały się raczej ilustracyjne niż prawdziwe badanie, które należy przeanalizować.
EdM