W modelu wielopoziomowym, jakie są praktyczne implikacje oszacowania w porównaniu z niedoszacowaniem parametrów korelacji efektu losowego?

W modelu wielopoziomowym, jakie są praktyczne i związane z interpretacją implikacje oszacowania w porównaniu z niedoszacowaniem parametrów korelacji efektu losowego? Praktycznym powodem pytania jest to, że w ramce Lmer w R nie ma zaimplementowanej metody szacowania wartości p za pomocą technik MCMC, gdy dokonuje się szacunków w modelu korelacji między parametrami.

Na przykład, patrząc na ten przykład (części cytowane poniżej), jakie są praktyczne implikacje M2 w porównaniu z M3. Oczywiście w jednym przypadku P5 nie zostanie oszacowane, aw drugim tak się stanie.

pytania

1. Ze względów praktycznych (chęć uzyskania wartości p za pomocą technik MCMC) można dopasować model bez korelacji między efektami losowymi, nawet jeśli P5 jest zasadniczo niezerowe. Jeśli to zrobi, a następnie oszacuje wartości p za pomocą techniki MCMC, czy wyniki są interpretowalne? (Wiem, @Ben Bolker wspominał wcześniej, że „łączenie testów istotności z MCMC jest trochę niespójne statystycznie, chociaż rozumiem potrzebę tego (uzyskanie przedziałów ufności jest bardziej możliwe)” , więc jeśli to pozwoli ci lepiej spać w nocy udawaj, że mówiłem o przedziałach ufności).
2. Jeśli nie uda się oszacować P5, czy to to samo, co stwierdzenie, że wynosi 0?
3. Jeśli P5 naprawdę jest niezerowy, to w jaki sposób wpływa to na oszacowane wartości P1-P4?
4. Jeśli P5 naprawdę jest niezerowy, to w jaki sposób wpływa to na oszacowanie błędu dla P1-P4?
5. Jeśli P5 naprawdę jest niezerowy, to w jaki sposób błędne są interpretacje modelu, który nie zawiera P5?

Pożyczanie od odpowiedzi @ Mike'a Lawrence'a (osoby bardziej kompetentne niż ja mogę to zastąpić pełną notacją modelu, nie jestem do końca pewien, czy mogę to zrobić z rozsądną wiernością):

M2: V1 ~ (1|V2) + V3 + (0+V3|V2)(Szacunki P1 - P4)

M3: V1 ~ (1+V3|V2) + V3(Oszacowania P1-P5)

Parametry, które można oszacować:

P1 : Globalny przechwytywanie

P2 : Przechwytuje efekt losowy dla V2 (tj. Dla każdego poziomu V2, odchylenie przechwytywania tego poziomu od globalnego przechwytywania)

P3 : Pojedyncza globalna ocena wpływu (nachylenia) V3

P4 : Wpływ V3 na każdym poziomie V2 (a dokładniej stopień, w jakim efekt V3 na danym poziomie odbiega od globalnego efektu V3), jednocześnie wymuszając zerową korelację między odchyleniami przechwytującymi a odchyleniami efektu V3 między poziomami V2.

P5 : Korelacja między odchyleniami przechwytującymi a odchyleniami V3 na różnych poziomach V2

Odpowiedzi uzyskane z wystarczająco dużej i szerokiej symulacji wraz z towarzyszącym kodem w języku R za pomocą lmera byłyby dopuszczalne.

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.

i

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.

1. Stałe parametry mają dokładnie takie same wartości w obu modelach.
2. 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.
3. 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).
4. Kowariancja i wariancja ustalonych efektów są różne.
5. 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.
6. 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)
    }

Placidia udzieliła już dokładnej odpowiedzi przy użyciu danych symulowanych na podstawie sleepstudyzestawu danych. Oto kolejna (mniej rygorystyczna) odpowiedź, która również wykorzystuje sleepstudydane.

Widzimy, że można wpływać na oszacowaną korelację między przypadkowym przechwyceniem a losowym nachyleniem poprzez „przesunięcie” losowej zmiennej predykcyjnej. Spójrz na wyniki z modeli fm1i fm2poniżej:

library(lmer)

#Fit Models
fm1 <- lmer(Reaction ~ Days + (Days | Subject), sleepstudy)
k <- 3 # Shift "Days" by an arbitrary amount
fm2 <- lmer(Reaction ~ I(Days + k) + (I(Days + k)| Subject), sleepstudy)

fm1 # Model Output
# Linear mixed model fit by REML ['lmerMod']
# Formula: Reaction ~ Days + (Days | Subject)
# Data: sleepstudy
# REML criterion at convergence: 1743.628
# Random effects:
#   Groups   Name        Std.Dev. Corr
# Subject  (Intercept) 24.740       
# Days         5.922   0.07
# Residual             25.592       
# Number of obs: 180, groups:  Subject, 18
# Fixed Effects:
#   (Intercept)         Days  
# 251.41        10.47

fm2 # Model Output
# Linear mixed model fit by REML ['lmerMod']
# Formula: Reaction ~ I(Days + k) + (I(Days + k) | Subject)
# Data: sleepstudy
# REML criterion at convergence: 1743.628
# Random effects:
#   Groups   Name        Std.Dev. Corr 
# Subject  (Intercept) 29.498        
# I(Days + k)  5.922   -0.55
# Residual             25.592        
# Number of obs: 180, groups:  Subject, 18
# Fixed Effects:
#   (Intercept)  I(Days + k)  
# 220.00        10.47

# Random effects from both models
cbind(ranef(fm1)$Subject,ranef(fm2)$Subject)
# (Intercept)        Days (Intercept) I(Days + k)
# 308   2.2585654   9.1989719 -25.3383538   9.1989727
# 309 -40.3985769  -8.6197032 -14.5394628  -8.6197043
# 310 -38.9602458  -5.4488799 -22.6136027  -5.4488807
# 330  23.6904985  -4.8143313  38.1334933  -4.8143315
# 331  22.2602027  -3.0698946  31.4698868  -3.0698946
# 332   9.0395259  -0.2721707   9.8560377  -0.2721706
# 333  16.8404311  -0.2236244  17.5113040  -0.2236243
# 334  -7.2325792   1.0745761 -10.4563076   1.0745761
# 335  -0.3336958 -10.7521591  31.9227854 -10.7521600
# 337  34.8903508   8.6282840   9.0054946   8.6282850
# 349 -25.2101104   1.1734142 -28.7303527   1.1734141
# 350 -13.0699567   6.6142050 -32.9125736   6.6142054
# 351   4.5778352  -3.0152572  13.6236077  -3.0152574
# 352  20.8635924   3.5360133  10.2555505   3.5360138
# 369   3.2754530   0.8722166   0.6588028   0.8722167
# 370 -25.6128694   4.8224646 -40.0802641   4.8224648
# 371   0.8070397  -0.9881551   3.7715053  -0.9881552
# 372  12.3145393   1.2840297   8.4624492   1.2840300

Na podstawie wyników modelu widzimy, że korelacja wariancji losowej uległa zmianie. Jednak nachylenia (stałe i losowe) pozostały takie same, jak oszacowano wariancję resztkową. Wartości przecięcia (stałe i losowe) zmieniły się w odpowiedzi na zmienną przesuniętą.

Dekorelacja losowej kowariancji przechwytywania-nachylenia dla LMM została omówiona w notatkach z wykładów dr Jacka Weissa tutaj . Weiss zauważa, że ​​zmniejszenie korelacji wariancji w ten sposób może czasem pomóc między innymi w konwergencji modeli.

Powyższy przykład zmienia losową korelację (parametr „P5”). Częściowo odnosząc się do Q3 PO, z powyższego wynika, że:

#   Parameter           Status
=================================
P1  Fixed Intercept     Affected
P2  Random Intercepts   Affected
P3  Fixed Slope         Not Affected
P4  Random Slopes       Not Affected
P5  Random Correlation  Affected