Rozważ dane dotyczące uśpienia zawarte w lme4. Bates omawia to w swojej książce online o lme4. W rozdziale 3 rozważa dwa modele danych.
M0:Reaction∼1+Days+(1|Subject)+(0+Days|Subject)
i
MA:Reaction∼1+Days+(Days|Subject)
W badaniu wzięło udział 18 osób, badanych przez okres 10 dni pozbawionych snu. Czasy reakcji obliczono na początku i na kolejne dni. Istnieje wyraźny wpływ między czasem reakcji a czasem trwania deprywacji snu. Istnieją również znaczące różnice między podmiotami. Model A dopuszcza możliwość interakcji między przypadkowym efektem przechwytywania a efektem nachylenia: wyobraźmy sobie, powiedzmy, że osoby o słabym czasie reakcji cierpią bardziej dotkliwie z powodu braku snu. Oznaczałoby to dodatnią korelację efektów losowych.
W przykładzie Batesa nie było widocznej korelacji z wykresu Kraty i nie było znaczącej różnicy między modelami. Jednak, aby zbadać postawione wyżej pytanie, postanowiłem wziąć dopasowane wartości z badania snu, zwiększyć korelację i spojrzeć na wydajność obu modeli.
Jak widać z obrazu, długi czas reakcji wiąże się z większą utratą wydajności. Korelacja zastosowana w symulacji wyniosła 0,58
Symulowałem 1000 próbek, stosując metodę symulacji w lme4, w oparciu o dopasowane wartości moich sztucznych danych. Dopasowałem M0 i Ma do każdego z nich i spojrzałem na wyniki. Oryginalny zestaw danych zawierał 180 obserwacji (10 dla każdego z 18 pacjentów), a symulowane dane mają tę samą strukturę.
Najważniejsze jest to, że różnica jest bardzo mała.
- Stałe parametry mają dokładnie takie same wartości w obu modelach.
- Losowe efekty są nieco inne. Istnieje 18 losowych efektów przechwytywania i 18 nachyleń dla każdej symulowanej próbki. Dla każdej próbki efekty te są zmuszane do dodania do 0, co oznacza, że średnia różnica między dwoma modelami wynosi (sztucznie) 0. Jednak wariancje i kowariancje różnią się. Mediana kowariancji pod MA wynosiła 104, wobec 84 pod M0 (wartość rzeczywista, 112). Wariancje nachyleń i przecięć były większe pod M0 niż MA, przypuszczalnie w celu uzyskania dodatkowego miejsca na drgania potrzebnego przy braku parametru wolnego kowariancji.
- Metoda ANOVA dla lmera daje statystykę F do porównywania modelu Slope z modelem tylko z przypadkowym przechwytywaniem (brak efektu z powodu braku snu). Oczywiście wartość ta była bardzo duża w obu modelach, ale zazwyczaj była (ale nie zawsze) większa w przypadku MA (średnia 62 vs. średnia 55).
- Kowariancja i wariancja ustalonych efektów są różne.
- Mniej więcej w połowie przypadków wie, że MA jest poprawne. Mediana wartości p dla porównania M0 z MA wynosi 0,0442. Pomimo obecności znaczącej korelacji i 180 zrównoważonych obserwacji właściwy model zostałby wybrany tylko w około połowie przypadków.
- Prognozowane wartości różnią się w obu modelach, ale bardzo nieznacznie. Średnia różnica między przewidywaniami wynosi 0, przy sd 2,7. Sd samych przewidywanych wartości wynosi 60,9
Dlaczego tak się dzieje? @gung domyślił się, że nieuwzględnienie możliwości korelacji zmusza losowe efekty do nieskorelowania. Być może powinno; ale w tej implementacji dozwolone są korelacje efektów losowych, co oznacza, że dane są w stanie wyciągnąć parametry we właściwym kierunku, niezależnie od modelu. Błędność niewłaściwego modelu pojawia się w prawdopodobieństwie, dlatego możesz (czasami) rozróżnić dwa modele na tym poziomie. Model efektów mieszanych zasadniczo dopasowuje regresje liniowe do każdego podmiotu, pod wpływem tego, co według modelu powinno być. Zły model wymusza dopasowanie mniej prawdopodobnych wartości niż w przypadku właściwego modelu. Ale parametry na koniec dnia zależą od dopasowania do rzeczywistych danych.
Oto mój nieco niezręczny kod. Chodziło o dopasowanie danych z badania snu, a następnie zbudowanie symulowanego zestawu danych o tych samych parametrach, ale z większą korelacją dla efektów losowych. Ten zestaw danych został przekazany do simulation.lmer () w celu symulacji 1000 próbek, z których każda pasowała w obie strony. Po sparowaniu dopasowanych obiektów mogłem wyciągnąć różne cechy dopasowania i porównać je za pomocą testów t lub cokolwiek innego.
# Fit a model to the sleep study data, allowing non-zero correlation
fm01 <- lmer(Reaction ~ 1 + Days +(1+Days|Subject), data=sleepstudy, REML=FALSE)
# Now use this to build a similar data set with a correlation = 0.9
# Here is the covariance function for the random effects
# The variances come from the sleep study. The covariance is chosen to give a larger correlation
sigma.Subjects <- matrix(c(565.5,122,122,32.68),2,2)
# Simulate 18 pairs of random effects
ranef.sim <- mvrnorm(18,mu=c(0,0),Sigma=sigma.Subjects)
# Pull out the pattern of days and subjects.
XXM <- model.frame(fm01)
n <- nrow(XXM) # Sample size
# Add an intercept to the model matrix.
XX.f <- cbind(rep(1,n),XXM[,2])
# Calculate the fixed effects, using the parameters from the sleep study.
yhat <- XX.f %*% fixef(fm01 )
# Simulate a random intercept for each subject
intercept.r <- rep(ranef.sim[,1], each=10)
# Now build the random slopes
slope.r <- XXM[,2]*rep(ranef.sim[,2],each=10)
# Add the slopes to the random intercepts and fixed effects
yhat2 <- yhat+intercept.r+slope.r
# And finally, add some noise, using the variance from the sleep study
y <- yhat2 + rnorm(n,mean=0,sd=sigma(fm01))
# Here is new "sleep study" data, with a stronger correlation.
new.data <- data.frame(Reaction=y,Days=XXM$Days,Subject=XXM$Subject)
# Fit the new data with its correct model
fm.sim <- lmer(Reaction ~ 1 + Days +(1+Days|Subject), data=new.data, REML=FALSE)
# Have a look at it
xyplot(Reaction ~ Days | Subject, data=new.data, layout=c(6,3), type=c("p","r"))
# Now simulate 1000 new data sets like new.data and fit each one
# using the right model and zero correlation model.
# For each simulation, output a list containing the fit from each and
# the ANOVA comparing them.
n.sim <- 1000
sim.data <- vector(mode="list",)
tempReaction <- simulate(fm.sim, nsim=n.sim)
tempdata <- model.frame(fm.sim)
for (i in 1:n.sim){
tempdata$Reaction <- tempReaction[,i]
output0 <- lmer(Reaction ~ 1 + Days +(1|Subject)+(0+Days|Subject), data = tempdata, REML=FALSE)
output1 <- lmer(Reaction ~ 1 + Days +(Days|Subject), data=tempdata, REML=FALSE)
temp <- anova(output0,output1)
pval <- temp$`Pr(>Chisq)`[2]
sim.data[[i]] <- list(model0=output0,modelA=output1, pvalue=pval)
}