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

26

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.

Chris
źródło

Odpowiedzi:

21

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.

Henrik
źródło
12

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.

Mark Miller
źródło
10

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.). wprowadź opis zdjęcia tutaj

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)

wprowadź opis zdjęcia tutaj

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.

alexforrence
źródło
1
+1. Oczywiście im mniejsza liczba osób, tym większa wariancja estymatora wariancji. Ale nie sądzę, żeby to miało znaczenie. Pytanie brzmi: jaka liczba osób pozwala uzyskać sensowne wyniki? Jeśli zdefiniujemy „niesensowny” wynik jako uzyskanie zerowej wariancji, wówczas w twojej symulacji zdarza się to dość często przy n = 5 lub mniej. Począwszy od n = 6 lub n = 7, prawie nigdy nie otrzymujesz dokładnej 0 wariancji, tzn. Model jest zbieżny z rozwiązaniem nie zdegenerowanym. Mój wniosek byłby taki, że n = 6 jest dopuszczalne na granicy.
ameba mówi Przywróć Monikę
1
BTW jest to zgodne z rpubs.com/bbolker/4187 .
ameba mówi Przywróć Monikę
8

„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.

StasK
źródło
1
Uważam, że odpowiedź na pytanie brzmi w pewnym stopniu „jak długi jest kawałek sznurka”. Ale nie pokładałbym dużej pewności w oszacowaniu średniej lub wariancji na podstawie próby mniejszej niż 15-20, więc ta sama ogólna zasada nie miałaby zastosowania do poziomów losowego efektu. Nigdy nie widziałem, aby ktoś uwzględniał identyfikator podmiotu jako stały i losowy efekt w badaniach podłużnych - czy to powszechna praktyka?
Chris
Oprócz niewielkiej liczby badanych w modelu mieszanym ich losowe efekty nie są obserwowane, więc musisz drażnić ich z danych, i prawdopodobnie potrzebujesz stosunkowo więcej danych, aby zrobić to niezawodnie niż w samym szacowaniu średniej i średniej wariancja, gdy wszystko jest obserwowane. Zatem 42 vs. 15-20 :). Myślę, że miałem na myśli losowe zbocza, ponieważ masz poprawne dane identyfikacyjne badanych, podając jedynie efekty losowe, w przeciwnym razie nie zostaną zidentyfikowane. Nawiasem mówiąc, ekonomiści nie wierzą w efekty losowe i publikują prawie wyłącznie to, co nazywają „efektami stałymi”, tj. Szacunki wewnątrz podmiotu.
StasK,
2
+1 @StasK za bardzo dobrą odpowiedź na pytanie, z którym bardzo trudno sobie poradzić. Wydaje mi się, że jest w tym odrobina niepotrzebnego sarkazmu i możesz zastanowić się nad edytowaniem swojej odpowiedzi, aby trochę bardziej szanować PO.
Michael R. Chernick,
@Michael, prawdopodobnie masz rację, że jest to markotna odpowiedź i prawdopodobnie niepotrzebnie. OP zaakceptował odpowiedź, którą chcieli usłyszeć, więc otrzymał rezolucję w tej sprawie. Poważniejsza odpowiedź wskazywałaby albo na dobry dowód symulacji, albo na analizę asymptotyczną wyższego rzędu, ale niestety nie znam takich odniesień.
StasK
3
Jeśli chodzi o wartość, myślę, że magiczna liczba „42” nie dotyczy tego, kiedy losowe efekty są uzasadnione, ale kiedy można uciec bez martwienia się o korekty o skończonych rozmiarach (np. Myśląc o skutecznych mianownikach stopniach swobody / korektach Kenwarda-Rogera / inne podobne podejścia).
Ben Bolker,
7

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ę.

Tom Wenseleers
źródło
4

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.

Charles
źródło
Czy zostało to przegłosowane z jakiegoś powodu związanego z jego treścią techniczną?
N Brouwer
Z jakiego rodzaju danymi pracujesz? Nie jestem pewien, dlaczego jesteś zaskoczony, słysząc, że model będzie zbieżny z 12-16 miarami na osobę. Nie mogę komentować błędu systematycznego w wynikowych modelach, ale nigdy nie miałem problemów z konwergencją w lme4modelach mieszanych i często uruchamiam je na próbkach o podobnych rozmiarach jak OP (pracuję również z zestawami danych biologicznych).
RTbecard
1

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.

m_e_s
źródło