Problem:
Przeczytałem w innych postach, które predict
nie są dostępne dla lmer
modeli 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 Recall
częstość niektórych słów, i kilka zmiennych objaśniających, w tym efekty losowe ( Auditorium
gdzie miał miejsce test; Subject
nazwa); i trwałe działanie , takie jak Education
, Emotion
(emocjonalne konotację słowa do zapamiętania) lub z spożywana przed badaniem.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
, Auditorium
i 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 Time
i 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 predict
funkcja w tym lmer
modelu? Najwyraźniej bierze to pod uwagę Time
zmienną, co skutkuje znacznie ściślejszym dopasowaniem oraz zygzakiem, który próbuje wyświetlić ten trzeci wymiar Time
zobrazowany na pierwszym wykresie.
Jeśli zadzwonię predict(fit2)
, dostanę 132.45609
pierwszy wpis, który odpowiada pierwszemu punktowi. Oto head
zbió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 .
źródło
predict
w 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ówpredict
. 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łowymlme4
.?predict
na konsoli [r], otrzymuję podstawową prognozę {statystyki} ...predict.merMod
jednak ... Jak widać na OP, zadzwoniłem po prostupredict
...lme4
pakiet, a następnie wpisz lme4 ::: Prognozuj.merMod, aby zobaczyć wersję specyficzną dla pakietu. Dane wyjściowelmer
są przechowywane w obiekcie klasymerMod
.predict
wie, 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ś.Odpowiedzi:
Łatwo się pomylić prezentacją współczynników podczas rozmowy
coef(fit2)
. Spójrz na podsumowanie fit2: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,ranef
aby uzyskać różnicę dla każdego przechwytywania z efektem losowym od średniego przechwytywania na wyższym poziomie zagnieżdżenia: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 jegoTime:Subject
współczynników:który moim zdaniem mieści się w błędzie zaokrąglenia wynoszącym 132,46.
Wartości zwracane przez przechwytują
coef
zdają się reprezentować ogólną osią oraz dodatkowoSubject
lubTime:Subject
specyficzne różnice, więc jest to trudniejsze do pracy z osobami; jeśli spróbujesz wykonać powyższe obliczenia zcoef
wartościami, podwójnie policzymy ogólny punkt przecięcia.źródło
ranef
od zbadaniu kodu R wlme4
. Wyjaśniłem prezentację w kilku miejscach.lme4
wyników. Przedstawione dane wydawały się raczej ilustracyjne niż prawdziwe badanie, które należy przeanalizować.