Załóżmy, że mam uczestników, z których każdy daje odpowiedź 20 razy, 10 w jednym stanie i 10 w innym. Dopasowuję liniowy model efektów mieszanych porównujący w każdych warunkach. Oto powtarzalny przykład symulujący tę sytuację za pomocą pakietu w :Y Tlme4
R
library(lme4)
fml <- "~ condition + (condition | participant_id)"
d <- expand.grid(participant_id=1:40, trial_num=1:10)
d <- rbind(cbind(d, condition="control"), cbind(d, condition="experimental"))
set.seed(23432)
d <- cbind(d, simulate(formula(fml),
newparams=list(beta=c(0, .5),
theta=c(.5, 0, 0),
sigma=1),
family=gaussian,
newdata=d))
m <- lmer(paste("sim_1 ", fml), data=d)
summary(m)
Model m
daje dwa ustalone efekty (punkt przecięcia i nachylenie dla warunku) oraz trzy efekty losowe (losowy punkt przechwytywania przez uczestnika, losowe nachylenie dla uczestnika dla warunku oraz korelacja między nachyleniem a nachyleniem).
Chciałbym statystycznie porównać rozmiar wariancji losowego przechwytywania według uczestników we wszystkich grupach określonych przez condition
(tj. Obliczyć komponent wariancji podświetlony na czerwono osobno w warunkach kontrolnych i eksperymentalnych, a następnie przetestować, czy różnica w wielkości komponentów jest niezerowa). Jak mam to zrobić (najlepiej w R)?
PREMIA
Powiedzmy, że model jest nieco bardziej skomplikowany: każdy uczestnik doświadcza 10 bodźców 20 razy, 10 w jednym stanie i 10 w innym. Istnieją zatem dwa zestawy skrzyżowanych efektów losowych: efekty losowe dla uczestnika i efekty losowe dla bodźca. Oto powtarzalny przykład:
library(lme4)
fml <- "~ condition + (condition | participant_id) + (condition | stimulus_id)"
d <- expand.grid(participant_id=1:40, stimulus_id=1:10, trial_num=1:10)
d <- rbind(cbind(d, condition="control"), cbind(d, condition="experimental"))
set.seed(23432)
d <- cbind(d, simulate(formula(fml),
newparams=list(beta=c(0, .5),
theta=c(.5, 0, 0, .5, 0, 0),
sigma=1),
family=gaussian,
newdata=d))
m <- lmer(paste("sim_1 ", fml), data=d)
summary(m)
Chciałbym porównać statystycznie wielkość losowej wariancji przechwytywania przez uczestnika we wszystkich grupach określonych przez condition
. Jak mam to zrobić i czy proces ten różni się od opisanego powyżej?
EDYTOWAĆ
Aby być bardziej szczegółowym na temat tego, czego szukam, chcę wiedzieć:
- Czy pytanie „czy warunkowe średnie odpowiedzi w ramach każdego warunku (tj. Losowe wartości przechwytywania w każdym warunku) są zasadniczo różne od siebie, poza tym, czego można by oczekiwać z powodu błędu próbkowania”, pytanie dobrze zdefiniowane (tj. Czy to pytanie nawet teoretycznie odpowiedzialny)? Jeśli nie, dlaczego nie?
- Jeśli odpowiedź na pytanie (1) brzmi „tak”, jak mam na nie odpowiedzieć? Wolałbym
R
implementację, ale nie jestem przywiązany dolme4
pakietu - na przykład wydaje się, żeOpenMx
pakiet ma możliwość obsługi analiz wielopoziomowych i wielopoziomowych ( https: //openmx.ssri.psu. edu / openmx-features ) i wydaje się, że jest to pytanie, na które należy odpowiedzieć w środowisku SEM.
źródło
Odpowiedzi:
Jest więcej niż jeden sposób na przetestowanie tej hipotezy. Na przykład procedura opisana przez @amoeba powinna działać. Wydaje mi się jednak, że najprostszym i najdogodniejszym sposobem przetestowania jest użycie starego dobrego testu współczynnika wiarygodności porównującego dwa modele zagnieżdżone. Jedyną potencjalnie trudną częścią tego podejścia jest umiejętność skonfigurowania pary modeli, tak aby porzucenie jednego parametru czysto przetestowało pożądaną hipotezę nierównych wariancji. Poniżej wyjaśniam, jak to zrobić.
Krótka odpowiedź
Przełącz na kodowanie kontrastu (suma na zero) dla zmiennej niezależnej, a następnie wykonaj test współczynnika prawdopodobieństwa, porównując pełny model z modelem, który wymusza korelację między losowymi nachyleniami i losowymi punktami przechwytywania na 0:
Wyjaśnienie wizualne / intuicja
Aby ta odpowiedź miała sens, musisz intuicyjnie rozumieć, co oznaczają różne wartości parametru korelacji dla obserwowanych danych. Rozważ (losowo zmieniające się) linie regresji specyficzne dla pacjenta. Zasadniczo parametr korelacji kontroluje, czy linie regresji uczestnika „rozkładają się w prawo” (korelacja dodatnia), czy „rozkładają się w lewo” (korelacja ujemna) względem punktuX= 0 , gdzie X jest niezależnym kodowanym kontrastem zmienna. Każda z tych implikuje nierówną wariancję średnich warunkowych odpowiedzi uczestników. Zilustrowano to poniżej:
Na tym wykresie pomijamy wiele obserwacji, które mamy dla każdego podmiotu w każdych warunkach i zamiast tego po prostu wykreślamy dwa losowe środki każdego podmiotu, z linią łączącą je, reprezentującą losowe nachylenie tego obiektu. (Są to dane 10 hipotetycznych podmiotów, a nie dane zamieszczone w PO).
W kolumnie po lewej stronie, gdzie występuje silna ujemna korelacja przechwytywania nachylenia, linie regresji rozkładają się w lewo względem punktuX= 0 . Jak widać wyraźnie na rysunku, to prowadzi do większej zmienności w badanych przypadkowych środków stan X= - 1 niż w stanie X= 1 .
Kolumna po prawej pokazuje odwrócone, lustrzane odbicie tego wzoru. W tym przypadku nie jest większe odchylenia w losowo badanych środków w stanX= 1 , niż w stanie X= - 1 .
Kolumna pośrodku pokazuje, co się dzieje, gdy losowe zbocza i przypadkowe przechwyty są nieskorelowane. Oznacza to, że linie regresji wychodzą w lewo dokładnie tak samo, jak w prawo, w stosunku do punktuX= 0 . Oznacza to, że wariancje badanych oznaczają, że w tych dwóch warunkach są równe.
Ważne jest tutaj, abyśmy zastosowali schemat kodowania kontrastu suma do zera, a nie fałszywe kody (to znaczy, nie ustawiając grup naX= 0 vs. X= 1 ). Jest tylko pod kontrastowego systemu kodowania, że mamy ten związek, w którym wariancje są równe wtedy i tylko wtedy, gdy korelacja osią nachylenie wynosi 0. Na rysunku poniżej próbuje zbudować że intuicji:
Ten rysunek pokazuje ten sam dokładny zestaw danych w obu kolumnach, ale z niezależną zmienną kodowaną na dwa różne sposoby. W kolumnie po lewej używamy kodów kontrastowych - dokładnie tak jest z pierwszej cyfry. W kolumnie po prawej używamy fałszywych kodów. Zmienia to znaczenie przechwyceń - teraz przechwyty reprezentują przewidywane odpowiedzi badanych w grupie kontrolnej. Dolny panel pokazuje konsekwencje tej zmiany, a mianowicie to, że korelacja przechwytywania nachylenia nie jest już nigdzie bliska zeru, mimo że dane są w głębokim znaczeniu takie same, a wariancje warunkowe są równe w obu przypadkach. Jeśli nadal nie wydaje się to mieć sensu, przestudiowanie mojej poprzedniej odpowiedzi, w której mówię więcej o tym zjawisku, może pomóc.
Dowód
Niechyja j k będzie jot tą odpowiedzią ja tego podmiotu pod warunkiem k . (Mamy tutaj tylko dwa warunki, więc k to tylko 1 lub 2.) Następnie można zapisać model mieszany
yja j k= αja+ βjaxk+ eja j k,
gdzie αja są podmiotami przechwytuje losowo i ma wariancję σ2)α , βja są losowym nachyleniem badanych i mają wariancję σ2)β , mija j k jest poziomem błędu poziomu obserwacji, a cov ( αja, βja) = σα β .
Chcemy pokazać, żevar ( αja+ βjax1) = var ( αja+ βjax2)) ⇔ σα β= 0
Począwszy od lewej strony tego wniosku, mamyvar ( αja+ βjax1)σ2)α+ x2)1σ2)β+ 2 x1σα βσ2)β( x2)1- x2)2)) + 2 σα β( x1- x2))= var ( αja+ βjax2))= σ2)α+ x2)2)σ2)β+ 2 x2)σα β= 0
Kody kontrastu suma do zera oznaczają, żex1+ x2)= 0 i x2)1= x2)2)= x2) . Następnie możemy dodatkowo zmniejszyć ostatnią linię powyższego do
σ2)β( x2)- x2)) + 2 σα β( x1+ x1)σα β= 0= 0 ,
co chcieliśmy udowodnić. (Aby ustalić inny kierunek implikacji, możemy po prostu wykonać te same kroki w odwrotnej kolejności).
Powtarzając, pokazuje to, że jeśli zmienna niezależna jest zakodowana kontrastem (suma do zera) , wówczas wariancje średnich losowych badanych w każdych warunkach są równe wtedy i tylko wtedy, gdy korelacja między losowymi nachyleniami i przypadkowymi przechwytywaniami wynosi 0. Klucz od tego wszystkiego wynika, że testowanie hipotezy zerowej, żeσα β= 0 , przetestuje hipotezę zerową o równych wariancjach opisanych przez OP.
źródło
(1 | subject)
dummy
Możesz przetestować istotność parametrów modelu za pomocą oszacowanych przedziałów ufności, dla których pakiet lme4 ma
confint.merMod
funkcję.ładowanie (patrz na przykład przedział ufności z ładowania )
profil prawdopodobieństwa (patrz na przykład jaki jest związek między prawdopodobieństwem profilu a przedziałami ufności? )
Istnieje również metoda,
'Wald'
ale jest ona stosowana tylko do stałych efektów.Istnieją również jakąś ANOVA (współczynnik prawdopodobieństwa) typu wypowiedzi w opakowaniu
lmerTest
, które jest imieniemranova
. Ale nie wydaje mi się, żeby to miało sens. Rozkład różnic w logLikelihood, kiedy hipoteza zerowa (zerowa wariancja dla efektu losowego) jest prawdziwa, nie jest rozkładem chi-kwadrat (być może, gdy liczba uczestników i prób jest wysoka, test współczynnika prawdopodobieństwa może mieć sens).Odchylenie w określonych grupach
Aby uzyskać wyniki dla wariancji w określonych grupach, możesz ponownie sparametryzować
Tam, gdzie dodaliśmy dwie kolumny do ramki danych (jest to potrzebne tylko, jeśli chcesz ocenić nieskorelowaną „kontrolę” i „eksperymentalną”, funkcja
(0 + condition || participant_id)
nie prowadziłaby do oceny różnych czynników w stanie jako nieskorelowanych)Teraz
lmer
poda wariancję dla różnych grupMożesz zastosować do nich metody profilowania. Na przykład teraz confint daje przedziały ufności dla wariancji kontrolnej i ćwiczeniowej.
Prostota
Możesz użyć funkcji wiarygodności, aby uzyskać bardziej zaawansowane porównania, ale istnieje wiele sposobów przybliżenia drogi (np. Możesz wykonać konserwatywny test anova / lrt, ale czy tego właśnie chcesz?).
W tym momencie zastanawiam się, jaki jest sens tego (nie tak powszechnego) porównania między wariancjami. Zastanawiam się, czy zaczyna być zbyt wyrafinowane. Dlaczego różnica między wariancjami zamiast stosunku między wariancjami (która odnosi się do klasycznego rozkładu F)? Dlaczego nie zgłosić przedziałów ufności? Musimy cofnąć się o krok i wyjaśnić dane oraz historię, którą powinien opowiedzieć, zanim przejdziemy do zaawansowanych ścieżek, które mogą być zbędne i luźne w kontakcie z materią statystyczną i względami statystycznymi, które w rzeczywistości są głównym tematem.
Zastanawiam się, czy należy zrobić coś więcej niż zwykłe określanie przedziałów ufności (które mogą w rzeczywistości powiedzieć znacznie więcej niż test hipotez. Test hipotez daje odpowiedź tak, nie, ale brak informacji o rzeczywistym rozprzestrzenianiu się populacji. Biorąc pod uwagę wystarczającą ilość danych, które możesz nieznaczną różnicę należy zgłosić jako znaczącą różnicę). Zagłębianie się bardziej w materię (w jakimkolwiek celu) wymaga, moim zdaniem, bardziej szczegółowego (wąsko zdefiniowanego) pytania badawczego, aby poprowadzić maszynę matematyczną do wprowadzenia odpowiednich uproszczeń (nawet jeśli dokładne obliczenia mogą być wykonalne lub gdy można to przybliżyć za pomocą symulacji / ładowania początkowego, nawet w niektórych ustawieniach nadal wymaga to odpowiedniej interpretacji). Porównaj z dokładnym testem Fishera, aby dokładnie rozwiązać (szczególne) pytanie (dotyczące tabel awaryjnych),
Prosty przykład
Aby podać przykład możliwej prostoty, pokazuję poniżej porównanie (symulacje) z prostą oceną różnicy między dwiema wariancjami grupowymi w oparciu o test F przeprowadzony przez porównanie wariancji poszczególnych średnich odpowiedzi i dokonany przez porównanie wariancje pochodne modelu mieszanego.
Możesz to zobaczyć w symulacji na poniższym wykresie, gdzie pomijając F-score na podstawie próbki oznacza, że F-score jest obliczany na podstawie przewidywanych wariancji (lub sum błędu kwadratu) z modelu.
Widać, że jest jakaś różnica. Różnica ta może wynikać z faktu, że model liniowy efektów mieszanych uzyskuje sumy błędu kwadratu (dla efektu losowego) w inny sposób. Te kwadratowe wyrażenia błędu nie są (już) dobrze wyrażone jako prosty rozkład chi-kwadrat, ale nadal są ściśle powiązane i można je aproksymować.
Model oparty na średnich jest więc bardzo dokładny. Ale jest mniej potężny. To pokazuje, że właściwa strategia zależy od tego, czego chcesz / potrzebujesz.
W powyższym przykładzie, gdy ustawisz prawą granicę ogona na 2.1 i 3.1, otrzymasz około 1% populacji w przypadku równej wariancji (odpowiednio 103 i 104 z 10 000 przypadków), ale w przypadku nierównej wariancji granice te różnią się dużo (dając 5334 i 6716 spraw)
kod:
źródło
sim_1 ~ condition + (0 + condition | participant_id)"
w którym to przypadku uzyskuje się parametryzację na dwa parametry (jeden dla każdej grupy) zamiast dwóch parametrów jeden dla przechwytywania i jeden dla efektu (który należy połączyć dla grup).car::linearHypothesisTest
( math.furman.edu/~dcs/courses/math47/R/library/car/html/... ), co pozwala użytkownikowi przetestować dowolne hipotezy za pomocą dopasowanego modelu. Musiałbym jednak użyć metody @ amoeba, aby uzyskać oba losowe przechwyty w tym samym dopasowanym modelu, aby można je było porównać z tą funkcją. Jestem również trochę niepewny co do ważności metody.Jednym stosunkowo prostym sposobem może być zastosowanie testów współczynnika wiarygodności,
anova
zgodnie z opisem wlme4
FAQ .Zaczynamy od pełnego modelu, w którym wariancje są nieograniczone (tzn. Dozwolone są dwie różne warianty), a następnie dopasowujemy jeden model ograniczony, w którym zakłada się, że dwie wariancje są równe. Po prostu je porównać z
anova()
(Zauważ, że mogę ustawićREML = FALSE
mimoREML = TRUE
zeanova(..., refit = FALSE)
jest całkowicie wykonalne ).Ten test jest jednak prawdopodobnie konserwatywny . Na przykład FAQ mówi:
Istnieje kilka alternatyw:
Symuluj prawidłową dystrybucję za pomocą
RLRsim
(jak opisano również w FAQ).Przedstawię drugą opcję w następujący sposób:
Jak widać, dane wyjściowe sugerują,
REML = TRUE
że uzyskalibyśmy dokładne wyniki. Ale pozostaje to ćwiczeniu dla czytelnika.Jeśli chodzi o premię, nie jestem pewien, czy
RLRsim
pozwala na jednoczesne testowanie wielu składników, ale jeśli tak, można to zrobić w ten sam sposób.Odpowiedź na komentarz:
Nie jestem pewien, czy na to pytanie można uzyskać rozsądną odpowiedź.
Czy losowe nachylenia wpływają na przechwytywanie losowe? W pewnym sensie może to mieć sens, ponieważ pozwalają każdemu poziomowi czynnika grupowania całkowicie idiosynkratyczny efekt dla każdego warunku. Na koniec szacujemy dwa idiosynkratyczne parametry dla dwóch warunków. Myślę jednak, że rozróżnienie między ogólnym poziomem przechwyconym przez przechwytywanie a efektem zależnym od warunków uchwyconym przez losowe nachylenie jest ważne, a wtedy losowe nachylenie nie może tak naprawdę wpływać na losowe przechwytywanie. Nadal jednak pozwala każdemu poziomowi czynnika grupowania osobno dla każdego poziomu warunku osobno.
Niemniej jednak mój test nadal spełnia oczekiwania pierwotnego pytania. Sprawdza, czy różnica wariancji między tymi dwoma warunkami wynosi zero. Jeśli wynosi zero, wówczas wariancje w obu warunkach są równe. Innymi słowy, tylko jeśli nie ma potrzeby losowego nachylenia, wariancja w obu warunkach jest identyczna. Mam nadzieję, że to ma sens.
źródło
contr.treatment
), dla których warunek kontrolny jest punktem odniesienia (tj. Dla którego obliczany jest losowy punkt przecięcia). Proponowaną przeze mnie parametryzację używam kontrastów sumarycznych (tj.contr.sum
), A przecięcie jest wielką średnią. Wydaje mi się, że bardziej sensowne jest sprawdzenie, czy różnica jest zerowa, gdy punkt przecięcia jest główną wartością zamiast warunku kontrolnego (ale napisanie go po zakończeniu sugeruje, że może być względnie nieistotne). Możesz przeczytać str. 24–26 z: singmann.org/download/publications/…condition
: umożliwia losowe przechwytywanie zmienianie się na różnych poziomachcondition
. Czy to prawda?m_full
kontram_full2b
. To znaczy: wariancje średnich warunkowych odpowiedzi uczestników w A vs. B są nierówne, jeśli losowa korelacja przechwytywania nachylenia jest niezerowa - co ważne, przy parametryzacji kodowania kontrastu suma na zero . Testowanie losowej wariancji nachylenia nie jest konieczne. Próbuję wymyślić, jak to zwięźle wyjaśnić ...Twój model
pozwala już na to, że wariancja między podmiotami w warunkach kontrolnych różni się od wariancji między podmiotami w warunkach eksperymentalnych. Można to wyjaśnić poprzez równoważną ponowną parametryzację:
Macierz kowariancji losowej ma teraz prostszą interpretację:
Tutaj dwie wariancje są dokładnie tymi dwiema wariantami, którymi jesteś zainteresowany: wariancją [między podmiotami] średnich warunkowych odpowiedzi w warunkach kontrolnych i takich samych w warunkach eksperymentalnych. W symulowanym zbiorze danych są to 0,25 i 0,21. Różnica jest podana przez
i jest równy 0,039. Chcesz sprawdzić, czy różni się znacznie od zera.
EDYCJA: Zdałem sobie sprawę, że test permutacji, który opisuję poniżej, jest niepoprawny; nie będzie działać zgodnie z przeznaczeniem, jeśli środki w warunkach kontrolnych / eksperymentalnych nie są takie same (ponieważ wtedy obserwacji nie można wymieniać pod zerą). Lepszym pomysłem może być ładowanie przedmiotów (lub przedmiotów / przedmiotów w przypadku premii) i uzyskanie przedziału ufności dla
delta
.Spróbuję naprawić poniższy kod, aby to zrobić.
Oryginalna sugestia oparta na permutacji (niepoprawna)
Często stwierdzam, że można zaoszczędzić sobie wielu kłopotów, wykonując test permutacji. Rzeczywiście, w tym przypadku konfiguracja jest bardzo łatwa. Permutujmy warunki kontrolne / eksperymentalne dla każdego osobnika osobno; wtedy wszelkie różnice w wariancjach powinny zostać wyeliminowane. Powtarzanie tego wiele razy da zerowy rozkład różnic.
(Nie programuję w języku R; wszyscy zachęcamy do ponownego napisania następujących elementów w lepszym stylu R.)
nrep
Dokładnie taką samą logikę można zastosować w przypadku bonusu.
źródło
sim_1 ~ 0 + condition + (0 + dummy(condition, "control") + dummy(condition, "experimental") | participant_id)
sformułowania i uzyskać taki sam wynik jak w mojej odpowiedzi.Odmiana przecięcia, która odpowiada wielkiemu środkowi, jest
Wariant kontrastu jest
A kowariancja między przechwyceniem a kontrastem jest
który, jak @Jake Westfall wyprowadził nieco inaczej, testuje hipotezę równych wariancji, gdy porównujemy model bez tego parametru kowariancji z modelem, w którym parametr kowariancji jest nadal zawarty / nie ustawiony na zero.
W szczególności wprowadzenie innego skrzyżowanego losowego czynnika grupującego (takiego jak bodźce) nie zmienia porównania modelu, które należy wykonać, tj.
anova(mod1, mod2)
(Opcjonalnie z argumentem,refit = FALSE
gdy używasz oszacowania REML), gdziemod1
imod2
są zdefiniowane tak jak @Jake Westfall.który testuje hipotezę, że wariancje w dwóch warunkach są równe i że są one równe (dodatniej) kowariancji między tymi dwoma warunkami.
Gdy mamy dwa warunki, model pasujący do macierzy kowariancji z dwoma parametrami w (dodatniej) złożonej strukturze symetrycznej można również zapisać jako
lub (przy użyciu zmiennej / współczynnika kategorialnego
condition
)z
Poniżej widzimy, że
mod1
,mod3
imod4
otrzymując równoważne pasuje:mod4
źródło