Jaka jest minimalna zalecana liczba grup dla współczynnika efektów losowych?

Używam modelu mieszanego w R( lme4) do analizy niektórych danych z powtarzanymi pomiarami. Mam zmienną odpowiedzi (zawartość włókna w kale) i 3 stałe efekty (masa ciała itp.). Moje badanie ma tylko 6 uczestników, z 16 powtarzanymi pomiarami dla każdego z nich (chociaż dwóch ma tylko 12 powtórzeń). Podmiotami są jaszczurki, którym podawano różne kombinacje jedzenia w różnych „zabiegach”.

Moje pytanie brzmi: czy mogę użyć ID przedmiotu jako losowego efektu?

Wiem, że jest to zwykły sposób działania w modelach mieszanych efektów podłużnych, w celu uwzględnienia losowo dobranej natury badanych oraz faktu, że obserwacje w obrębie badanych będą bardziej skorelowane niż między badanymi. Lecz traktowanie identyfikatora podmiotu jako efektu losowego obejmuje oszacowanie średniej i wariancji dla tej zmiennej.

• Ponieważ mam tylko 6 osób (6 poziomów tego czynnika), czy to wystarczy, aby uzyskać dokładną charakterystykę średniej i wariancji?

• Czy fakt, że mam kilka powtórzonych pomiarów dla każdego przedmiotu, pomaga w tym względzie (nie rozumiem, jak to się liczy)?

• Wreszcie, jeśli nie mogę użyć identyfikatora podmiotu jako efektu losowego, czy włączenie go jako efektu stałego pozwoli mi kontrolować fakt, że powtarzałem pomiary?

Edycja: Chciałbym tylko wyjaśnić, że kiedy mówię „czy mogę” użyć identyfikatora podmiotu jako losowego efektu, mam na myśli „czy to dobry pomysł”. Wiem, że mogę dopasować model tylko z 2 poziomami, ale na pewno byłoby to nie do obrony? Pytam, w którym momencie sensownie jest myśleć o traktowaniu pacjentów jako efektów losowych? Wydaje się, że literatura mówi, że 5-6 poziomów jest dolną granicą. Wydaje mi się, że szacunki średniej i wariancji efektu losowego nie byłyby bardzo precyzyjne, dopóki nie będzie 15+ poziomów czynników.

Krótka odpowiedź: Tak, możesz użyć ID jako efektu losowego z 6 poziomami.

Nieco dłuższa odpowiedź: FAQ GLMM @ BenBolkera mówi (między innymi), co następuje pod nagłówkiem „ Czy powinienem traktować czynnik xxx jako stały czy losowy? ”:

Jednym ze szczególnie istotnych punktów dla „nowoczesnego” oszacowania modelu mieszanego (zamiast „klasycznej” metody chwilowej) jest to, że dla celów praktycznych musi istnieć rozsądna liczba poziomów efektów losowych (np. Bloków) - więcej niż Co najmniej 5 lub 6.

Jesteś więc w dolnej granicy, ale po jej prawej stronie.

Próbując obliczyć minimalną liczbę grup dla modelu wielopoziomowego, zajrzałem do książki Analiza danych za pomocą regresji i modeli Mulitilevel / Hierarchical autorstwa Gelmana i Hilla (2007).

Wydaje się, że zajmują się tym tematem w rozdziale 11, sekcja 5 (strona 247), gdzie piszą, że gdy jest <5 grup, wówczas modele wielopoziomowe zwykle niewiele dodają w stosunku do klasycznych modeli. Wydaje się jednak, że piszą, że stosowanie modelu wielopoziomowego jest niewielkie.

Ci sami autorzy wydają się powracać do tego tematu w rozdziale 12, sekcja 9 (strony 275–276). Tam piszą, że porady dotyczące minimalnej liczby grup dla modelu wielopoziomowego są błędne. Tam ponownie mówią, że modele wielopoziomowe często dodają niewiele więcej niż modele klasyczne, gdy liczba grup jest niewielka. Piszą jednak również, że modele wielopoziomowe nie powinny robić nic gorszego niż regresja bez puli (gdzie wydaje się, że brak puli oznacza, że ​​wskaźniki klasowe są stosowane w regresji klasycznej).

Na stronach 275–276 autorzy mają konkretny podrozdział dotyczący jednej lub dwóch grup (np. Mężczyzna kontra kobieta). Piszą tutaj, że zazwyczaj wyrażają model w klasycznej formie. Twierdzą jednak, że modelowanie wielopoziomowe może być przydatne nawet w przypadku jednej lub dwóch grup. Piszą, że w przypadku jednej lub dwóch grup modelowanie wielopoziomowe sprowadza się do regresji klasycznej.

Odnoszę z tego wrażenie, że regresja klasyczna stanowi jeden koniec ciągłości modeli, tj. Szczególny przypadek modelu wielopoziomowego.

Na podstawie powyższego mam wrażenie, że regresja klasyczna i modelowanie wielopoziomowe zwrócą prawie identyczne oszacowania, gdy istnieją tylko dwie grupy, i że stosowanie modeli wielopoziomowych z tylko jedną, dwiema, trzema, pięcioma lub pięcioma grupami jest w porządku.

Spróbuję zmodyfikować tę odpowiedź w przyszłości za pomocą Rkodu i małego zestawu danych, porównując szacunki uzyskane przy obu podejściach przy użyciu dwóch grup.

Dla tego, co jest warte, przeprowadziłem trochę badań symulacyjnych, aby sprawdzić stabilność oszacowania wariancji dla stosunkowo prostego LMM (używając sleepstudyzestawu danych dostępnego przez lme4). Pierwszy sposób generuje wszystkie możliwe kombinacje tematów dla ngroupsliczby tematów i dostosowuje model dla każdej możliwej kombinacji. Drugi obejmuje kilka losowych podzbiorów podmiotów.

library(lme4)
library(ggplot2)
library(tidyr)

m0 <- lmer(Reaction ~ Days + (1|Subject), data = sleepstudy,
           control = lmerControl(optimizer = "nloptwrap"))
# set the number of factor levels
ngroups <- 3:18 
# generate all possible combinations
combos <- lapply(X = ngroups, 
                 FUN = function(x) combn(unique(sleepstudy$Subject), x)) 

# allocate output (sorry, this code is entirely un-optimized)
out <- list(matrix(NA, ncol(combos[[1]]), 1), matrix(NA, ncol(combos[[2]]), 1),
            matrix(NA, ncol(combos[[3]]), 1), matrix(NA, ncol(combos[[4]]), 1),
            matrix(NA, ncol(combos[[5]]), 1), matrix(NA, ncol(combos[[6]]), 1),
            matrix(NA, ncol(combos[[7]]), 1), matrix(NA, ncol(combos[[8]]), 1),
            matrix(NA, ncol(combos[[9]]), 1), matrix(NA, ncol(combos[[10]]), 1),
            matrix(NA, ncol(combos[[11]]), 1), matrix(NA, ncol(combos[[12]]), 1),
            matrix(NA, ncol(combos[[13]]), 1), matrix(NA, ncol(combos[[14]]), 1),
            matrix(NA, ncol(combos[[15]]), 1), matrix(NA, ncol(combos[[16]]), 1))
# took ~ 2.5 hrs on my laptop, commented out for safety
#system.time(for(ii in 1:length(combos)) {
#    for(jj in 1:ncol(combos[[ii]])) {
#    sls <- sleepstudy[sleepstudy$Subject %in% combos[[ii]][,jj],]
#    out[[ii]][jj] <- attr(VarCorr(update(m0, data = sls))$Subject, 'stddev')
#        }
#    })

# pad with zeros, not all were equal
# from http://stackoverflow.com/questions/11148429/r-convert-asymmetric-list-to-matrix-number-of-elements-in-each-sub-list-diffe
max.len <- max(sapply(out, length))
corrected.list <- lapply(out, function(x) {c(x, rep(NA, max.len - length(x)))})
mat <- do.call(rbind, corrected.list)
mat <- data.frame(t(mat))
names(mat) <- paste0('s',3:18)
mat <- gather(mat, run, value)

ggplot(mat, aes(x = value, fill = run)) + 
    geom_histogram(bins = 60) +
    geom_vline(xintercept = 37.12, linetype =  'longdash', 
               aes(colour = 'original')) +
    facet_wrap(~run, scales = 'free_y') +
    scale_x_continuous(breaks = seq(0, 100, by = 20)) + 
    theme_bw() + 
    guides(fill = FALSE)

Kropkowana czarna linia jest oryginalnym oszacowaniem wariancji punktowej, a aspekty reprezentują różne liczby podmiotów ( s3będące grupami trzech podmiotów, s4będących czterema itp.).

I alternatywny sposób:

ngroups <- 3:18
reps <- 500
out2<- matrix(NA, length(ngroups), reps)

for (ii in 1:length(ngroups)) {
    for(j in 1:reps) {
        sls <- sleepstudy[sleepstudy$Subject %in% sample(unique(sleepstudy$Subject), ngroups[i], replace = FALSE),]
        out2[i,j] <- attr(VarCorr(update(m0, data = sls))$Subject, 'stddev')
    }
}
out2 <- data.frame(t(out2))
names(out2) <- paste0('s',3:18)
out2 <- gather(out2, run, value)

ggplot(out2, aes(x = value, fill = run)) + 
    geom_histogram(bins = 60) +
    geom_vline(xintercept = 37.12, linetype =  'longdash', 
               aes(colour = 'original')) +
    facet_wrap(~run, scales = 'free_y') +
    scale_x_continuous(breaks = seq(0, 100, by = 20)) + 
    theme_bw() + 
    guides(fill = FALSE)

Wydaje się (w każdym razie dla tego przykładu), że wariancja tak naprawdę się nie stabilizuje, dopóki nie będzie co najmniej 14 osób, jeśli nie później.

„Przeważnie nieszkodliwa ekonometria” autorstwa Angrista i Pischke'a ma sekcję zatytułowaną „Mniej niż 42 gromady”, w której częściowo żartobliwie mówią:

Dlatego, zgodnie z ... powiedzeniem, że odpowiedź na życie, wszechświat i wszystko ma 42, uważamy, że pytanie brzmi: ile klastrów wystarcza do wiarygodnego wnioskowania przy użyciu standardowej korekty klastrów [podobnie jak estymator wariancji w GEE]?

Mój instruktor ekonometrii odpowiadał na pytania takie jak twoje: „Ameryka jest wolnym krajem, możesz robić, co chcesz. Ale jeśli chcesz opublikować swój artykuł, musisz być w stanie obronić to, co zrobiłeś. „ Innymi słowy, prawdopodobnie będziesz w stanie uruchomić kod R lub Stata lub HLM lub Mplus lub SAS PROC GLIMMIX z 6 podmiotami (i przełączyć się na te alternatywne pakiety, jeśli jeden z wybranych nie uruchomi tego), ale prawdopodobnie będziesz mieć bardzo trudny czas obrony tego podejścia i uzasadnienia testów asymptotycznych.

Uważam, że domyślnie uwzględnienie zmiennej jako losowego nachylenia oznacza również uwzględnienie tego jako efektu stałego i musisz przeskakiwać przez wiele pętli składniowych, jeśli chcesz mieć to jako efekt losowy za pomocą zero. To rozsądny wybór, który dokonali dla ciebie programiści.

Można również użyć mieszanego modelu bayesowskiego - w takim przypadku niepewność w oszacowaniu efektów losowych jest w pełni uwzględniona w obliczeniach wiarygodnych przedziałów prognozy 95%. Na przykład nowy pakiet brmsi funkcja R brmumożliwiają bardzo łatwe przejście z lme4częstego modelu mieszanego do modelu Bayesowskiego, ponieważ ma on prawie identyczną składnię.

Nie użyłbym modelu efektów losowych z tylko 6 poziomami. Modele wykorzystujące 6-poziomowy efekt losowy mogą być czasami uruchamiane przy użyciu wielu programów statystycznych i czasami dają obiektywne szacunki, ale:

1. Myślę, że w społeczności statystycznej panuje arbitralny pogląd, że 10-20 to liczba minimalna. Jeśli chcesz opublikować swoje badania, radzimy poszukać czasopisma bez przeglądu statystycznego (lub być w stanie uzasadnić swoją decyzję dość wyrafinowanym językiem).
2. Przy tak małej liczbie klastrów wariancja między klastrami może być źle oszacowana. Złe oszacowanie wariancji między klastrami zwykle przekłada się na słabą ocenę standardowego błędu współczynników zainteresowania. (modele efektów losowych opierają się na liczbie klastrów teoretycznie zmierzających do nieskończoności).
3. Często modele po prostu się nie zbiegają. Czy próbowałeś uruchomić swój model? Zaskoczyłbym tylko 12-16 miarami na temat, że modele są zbieżne. Kiedy udało mi się przekonać tego rodzaju model do zbiegania się, miałem setki pomiarów na klaster.

Ten problem został rozwiązany w większości standardowych podręczników w terenie, a w pewnym stopniu rozwiązałeś je w swoim pytaniu. Nie sądzę, że dam ci nowe informacje.

Minęło dużo czasu od pierwotnego pytania, ale pomyślałem, że mogę dodać kilka punktów dotyczących wyboru modelu.

1 - Dopóki model jest zidentyfikowany (tzn. Masz stopnie swobody w przestrzeni parametrów), powinieneś być w stanie WYPRÓBOWAĆ, aby dopasować model. W zależności od metody optymalizacji model może być zbieżny lub nie. W każdym razie nie próbowałbym uwzględnić więcej niż 1 lub 2 losowych efektów i zdecydowanie nie więcej niż 1 interakcja między poziomami. W konkretnym przypadku przedstawionego tutaj problemu, jeśli podejrzewamy, że interakcja między charakterystycznymi cechami jaszczurki (np. Wiekiem, rozmiarem itp.) A grupą cech charakterystycznych leczenia / miary 6 może nie wystarczyć do dokonania wystarczająco dokładnych szacunków.

2 - Jak wspomniano w kilku odpowiedziach, problemem może być konwergencja. Jednak moje doświadczenie jest takie, że chociaż dane z nauk społecznych mają ogromny problem z konwergencją z powodu problemów pomiarowych, nauki przyrodnicze, a zwłaszcza biochemiczne powtarzane pomiary, mają znacznie mniejsze błędy standardowe. Wszystko zależy od procesu generowania danych. W danych społecznych i ekonomicznych musimy pracować na różnych poziomach abstrakcji. W przypadku danych biologicznych i chemicznych, a z pewnością astronomicznych błąd pomiaru danych stanowi mniejszy problem.