Mam do analizy zestaw danych z niezrównoważonymi powtarzanymi pomiarami i przeczytałem, że sposób, w jaki większość pakietów statystycznych obsługuje to z ANOVA (tj. Suma kwadratów typu III) jest błędny. Dlatego chciałbym użyć modelu mieszanych efektów do analizy tych danych. Dużo czytałem o modelach mieszanych R, ale wciąż jestem Rnowicjuszem w modelach z efektami mieszanymi i nie jestem zbyt pewny, że robię wszystko dobrze. Zauważ, że nie mogę jeszcze całkowicie rozwieść się z „tradycyjnymi” metodami i nadal potrzebuję wartości i testów post hoc.$p$

Chciałbym wiedzieć, czy poniższe podejście ma sens, czy też robię coś strasznie złego. Oto mój kod:

# load packages
library(lme4)
library(languageR)
library(LMERConvenienceFunctions)
library(coda)
library(pbkrtest)

# import data
my.data <- read.csv("data.csv")

# create separate data frames for each DV & remove NAs
region.data <- na.omit(data.frame(time=my.data$time, subject=my.data$subject, dv=my.data$dv1))

# output summary of data
data.summary <- summary(region.data)

# fit model
# "time" is a factor with three levels ("t1", "t2", "t3")
region.lmer <- lmer(dv ~ time + (1|subject), data=region.data)

# check model assumptions
mcp.fnc(region.lmer)

# remove outliers (over 2.5 standard deviations)
rm.outliers <- romr.fnc(region.lmer, region.data, trim=2.5)
region.data <- rm.outliers$data
region.lmer <- update(region.lmer)

# re-check model assumptions
mcp.fnc(region.lmer)

# compare model to null model
region.lmer.null <- lmer(dv ~ 1 + (1|subject), data=region.data)
region.krtest <- KRmodcomp(region.lmer, region.lmer.null)

# output lmer summary
region.lmer.summary <- summary(region.lmer)

# run post hoc tests
t1.pvals <- pvals.fnc(region.lmer, ndigits=10, withMCMC=TRUE)

region.lmer <- lmer(dv ~ relevel(time,ref="t2") + (1|subject), data=region.data)
t2.pvals <- pvals.fnc(region.lmer, ndigits=10, withMCMC=TRUE)

region.lmer <- lmer(dv ~ relevel(time,ref="t3") + (1|subject), data=region.data)
t3.pvals <- pvals.fnc(region.lmer, ndigits=10, withMCMC=TRUE)

# Get mcmc mean and 50/95% HPD confidence intervals for graphs
# repeated three times and stored in a matrix (not shown here for brevity)
as.numeric(t1.pvals$fixed$MCMCmean)
as.numeric(t1.pvals$fixed$HPD95lower)
as.numeric(t1.pvals$fixed$HPD95upper)
HPDinterval(as.mcmc(t1.pvals$mcmc),prob=0.5)
    HPDinterval(as.mcmc(t1.pvals$mcmc),prob=0.5)

Mam kilka konkretnych pytań:

1. Czy to prawidłowy sposób analizowania modeli efektów mieszanych? Jeśli nie, co powinienem zamiast tego zrobić.
2. Czy wykresy krytyki wyprowadzane przez mcp.fnc są wystarczająco dobre do weryfikacji założeń modelu, czy też powinienem podjąć dodatkowe kroki.
3. Rozumiem, że aby modele mieszane były ważne, dane muszą szanować założenia normalności i homoscedastyczności. Jak ocenić, co jest „w przybliżeniu normalne”, a co nie, patrząc na wykresy krytyki wygenerowane przez mcp.fnc? Czy muszę to po prostu wyczuć, czy jest to przepisany sposób robienia rzeczy? Jak solidne są modele mieszane w odniesieniu do tych założeń?
4. Muszę ocenić różnice między trzema punktami czasowymi dla ~ 20 cech (biomarkerów) badanych w mojej próbie. Dopasowanie i testowanie osobnych modeli dla każdego możliwego do zaakceptowania, o ile zgłaszam wszystkie podjęte testy (znaczące lub nie), czy też potrzebuję jakiejkolwiek formy korekty do wielu porównań.

Aby być nieco bardziej precyzyjnym w odniesieniu do eksperymentu, oto kilka szczegółów. Obserwowaliśmy wielu uczestników, którzy poddawali się leczeniu. Zmierzyliśmy wiele biomarkerów przed rozpoczęciem leczenia i po dwóch punktach czasowych. Chciałbym zobaczyć, czy między tymi trzema punktami czasowymi istnieją różnice w tych biomarkerach.

Opieram większość tego, co robię tutaj, na tym samouczku , ale wprowadziłem pewne zmiany w zależności od moich potrzeb i rzeczy, które czytam. Wprowadzone przeze mnie zmiany to:

1. przypisać współczynnik „czasu” do uzyskania porównań t1-t2, t2-t3 i t1-t3 z pvals.fnc (z pakietu languageR)
2. porównaj mój model mieszany z modelem zerowym przy użyciu przybliżonego testu F opartego na podejściu Kenwarda-Rogera (przy użyciu pakietu pbkrtest) zamiast testu współczynnika wiarygodności (ponieważ czytam, że Kenward-Roger jest teraz lepiej postrzegany)
3. Użyj pakietu LMERConvenienceFunctions, aby sprawdzić założenia i usunąć wartości odstające (ponieważ czytam, że modele mieszane są bardzo wrażliwe na wartości odstające)

Twoje pytanie jest trochę „duże”, więc zacznę od ogólnych komentarzy i wskazówek.

Trochę czytania w tle i przydatne pakiety

Prawdopodobnie powinieneś rzucić okiem na niektóre wprowadzenie do samouczka dotyczące korzystania z modeli mieszanych, polecam zacząć od Baayen i wsp. (2008) - autor jest Baayen languageR. Barr i wsp. (2013) omawiają niektóre kwestie ze strukturą efektów losowych, a Ben Bolker jest jednym z obecnych twórców lme4. Baayen i in. Są szczególnie dobrzy do twoich pytań, ponieważ spędzają dużo czasu omawiając stosowanie testów ładowania / permutacji (rzeczy za nimi mcp.fnc).

Literatura

• Baayen, RH, D. Davidson und D. Bates (2008). Modelowanie efektów mieszanych ze skrzyżowanymi efektami losowymi dla przedmiotów i przedmiotów . Journal of Memory and Language, 59 (390–412).
• Barr, Dale J, R. Levy, C. Scheepers und HJ Tily (2013). Struktura efektów losowych do testowania hipotez potwierdzających: zachowaj maksymalną wartość . Journal of Memory and Language, 68: 255–278.
• Bolker, Benjamin M, ME Brooks, CJ Clark, SW Geange, JR Poulsen, MHH Stevens und J.-SS White (2009). Uogólnione liniowe modele mieszane: praktyczny przewodnik po ekologii i ewolucji . Trendy Ecol Evol, 24 (3): 127–35.

Florian Jaeger ma również sporo artykułów na temat praktycznych modeli mieszanych na blogu swojego laboratorium .

Istnieje również wiele przydatnych pakietów R do wizualizacji i testowania modeli mieszanych, takich jak lmerTesti effects. effectsPakiet jest szczególnie miłe, ponieważ pozwala wykreślić liniowe przedziały ufności regresji a dzieje się za kulisami.

Dobroć dopasowania i porównywanie modeli

$p$-wartości są naprawdę kiepskim sposobem porównywania modeli mieszanych (i ogólnie nie są świetną metodą do niczego). Istnieje powód, dla którego Doug Bates (oryginalny autor lme4) nie uważa, że ​​konieczne jest ich włączenie. Jeśli naprawdę chcesz pójść tą drogą (a błagam, abyś tego nie robił), wspomniane powyżej lmerTestzapewnia dodatkowe możliwości obliczania i radzenia sobie z nimi.

Lepszymi sposobami porównywania modeli są logarytmiczne prawdopodobieństwo lub różne kryteria teoretyczne, takie jak AIC i BIC. W przypadku AIC i BIC ogólna zasada jest taka, że ​​„im mniejszy, tym lepszy”, ale trzeba tam uważać, ponieważ obaj karzą za więcej stopni swobody modelu i nie ma „absolutnej” dobrej lub złej wartości. Aby uzyskać ładne podsumowanie modeli AIC i modeli wiarygodności dziennika, możesz użyć anova()funkcji, która została przeciążona do przyjmowania merobiektów. Uwaga:$\chi^2$podane wartości są ważne tylko dla porównań między modelami zagnieżdżonymi . Niemniej jednak wynik jest całkiem przyjemny, ponieważ jest tak tabelaryczny, nawet w przypadku innych porównań. W przypadku modeli zagnieżdżonych masz fajny$\chi^2$ przetestuj i nie potrzebujesz $p$-wartości, aby bezpośrednio porównać dwa modele. Wadą tego jest to, że nie jest od razu jasne, jak dobre jest twoje dopasowanie.

Aby spojrzeć na poszczególne efekty (tj. Rzeczy, które zobaczyłbyś w tradycyjnej ANOVA), powinieneś spojrzeć na $t$-wartości dla ustalonych efektów w modelach (które są częścią summary()jeśli się nie mylę). Ogólnie wszystko$|t|>2$jest uważany za dobry (więcej szczegółów w Baayen i in.). Możesz również uzyskać dostęp do ustalonych efektów bezpośrednio za pomocą funkcji pomocnika fixef().

Powinieneś także upewnić się, że żaden z twoich ustalonych efektów nie jest zbyt silnie skorelowany - wielokoliniowość narusza założenia modelu. Florian Jaeger napisał o tym trochę , a także kilka możliwych rozwiązań.

Wiele porównań

Odpowiem na twoje pytanie 4 bardziej bezpośrednio. Odpowiedź jest w zasadzie taka sama jak dobra praktyka z tradycyjną ANOVA, niestety wydaje się, że jest to miejsce, w którym większość badaczy jest bardzo niepewna. (Jest to to samo, co tradycyjna ANOVA, ponieważ zarówno liniowe modele mieszanych efektów, jak i ANOVA są oparte na ogólnym modelu liniowym, po prostu modele mieszane mają dodatkowy termin dla efektów losowych.) Jeśli zakładasz, że okna czasowe tworzą różnicę i chcesz porównać efekty czasu, należy uwzględnić czas w swoim modelu. Nawiasem mówiąc, zapewni to również wygodny sposób oceny, czy czas zrobił różnicę: czy istnieje główny (stały) efekt czasu? Wadą tej trasy jest to, że Twój model stanie się dużo bardziej skomplikowany, a singiel „super” model z czasem jako parametrem prawdopodobnie potrwa dłużej niż trzy mniejsze modele bez czasu jako parametru. Rzeczywiście, klasyczny przykład samouczka dla modeli mieszanych tosleepstudy który wykorzystuje czas jako parametr.

Nie mogłem do końca powiedzieć, czy różne cechy mają być predyktorami, czy zmiennymi zależnymi. Jeśli mają być predyktorami, możesz rzucić je wszystkie w jeden model i zobaczyć, który ma największy$t$-wartość, ale ten model byłby niezwykle złożony, nawet jeśli nie zezwalasz na interakcje i obliczenie zajmuje dużo czasu. Szybszym i potencjalnie łatwiejszym sposobem byłoby obliczenie różnych modeli dla każdego predyktora. Polecam użycie foreachpętli do równoległego umieszczenia jej na tylu rdzeniach, ile masz - lme4nie wykonuje ona równolegle swoich operacji macierzowych, więc możesz obliczyć inny model dla każdego rdzenia równolegle. (Zobacz krótkie wprowadzenie tutaj ). Następnie możesz porównać AIC i BIC każdego modelu, aby ustalić, który predyktor jest najlepszy. (W tym przypadku nie są zagnieżdżone, więc „$\chi^2$ Statystyczny.)

Jeśli twoje cechy są zmienną zależną, to i tak będziesz musiał obliczyć różne modele, a następnie możesz użyć AIC i BIC do porównania wyników.