14

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 : $N$ $Y$ $Y$ lme4R

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 mdaje 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 Rimplementację, ale nie jestem przywiązany do lme4pakietu - na przykład wydaje się, że OpenMxpakiet 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 mixed-model lme4-nlme

— Patrick S. Forscher
źródło

1

@MarkWhite, zaktualizowałem pytanie w odpowiedzi na twoje komentarze. Chodzi mi o to, że chcę porównać odchylenie standardowe przechwytywania uczestnika, gdy dają odpowiedzi w warunkach kontrolnych, w porównaniu do odpowiedzi w warunkach eksperymentalnych. Chcę to zrobić statystycznie, tj. Sprawdzić, czy różnica w odchyleniu standardowym przechwyceń jest różna od 0.

— Patrick S. Forscher,

2

Napisałem odpowiedź, ale będę spał na niej, ponieważ nie jestem pewien, czy jest ona bardzo przydatna. Pytanie sprowadza się do tego, że nie sądzę, że można zrobić to, o co prosisz. Losowym efektem przechwytywania jest wariancja w środkach uczestników, gdy są oni w stanie kontrolnym. Nie można więc patrzeć na wariancję obserwacji obserwowanych w warunkach eksperymentalnych. Przechwyty są definiowane na poziomie osoby, a warunek na poziomie obserwacji. Jeśli próbujesz porównać wariancje między warunkami, pomyślałbym o warunkowo heteroscedastycznych modelach.

— Mark White,

2

Pracuję nad poprawką i ponownym przesłaniem artykułu, w którym mam uczestników, którzy udzielają odpowiedzi na zestawy bodźców. Każdy uczestnik jest narażony na wiele warunków, a każdy bodziec odbiera odpowiedź w wielu warunkach - innymi słowy, moje badanie emuluje konfigurację, którą opisuję w moim opisie „BONUS”. Na jednym z moich wykresów wydaje się, że średnia odpowiedź uczestnika ma większą zmienność w jednym z warunków niż w innych. Recenzent poprosił mnie o sprawdzenie, czy to prawda.

— Patrick S. Forscher,

2

Zobacz tutaj stats.stackexchange.com/questions/322213, aby dowiedzieć się, jak skonfigurować model lme4 z różnymi parametrami wariancji dla każdego poziomu zmiennej grupującej. Nie jestem pewien, jak wykonać test hipotez, czy dwa parametry wariancji są równe; osobiście zawsze wolałbym bootstrapować nad badanymi i bodźcami, aby uzyskać przedział ufności, lub może ustawić jakiś test hipotezowy podobny do permutacji (oparty na ponownym próbkowaniu).

— ameba mówi Przywróć Monikę

3

Zgadzam się z komentarzem @ MarkWhite, że pytanie „czy losowe warianty przechwytywania zasadniczo różnią się od siebie ...” jest w najlepszym razie niejasne, a w najgorszym przypadku bezsensowne, ponieważ przecięcie musi dotyczyć wartości Y w jednej określonej grupie ( grupa przypisała wartość 0), więc porównywanie „przechwyceń” w grupach ściśle mówiąc nie ma sensu. Myślę, że lepszym sposobem na sformułowanie twojego pytania, tak jak je rozumiem, byłoby coś takiego: „czy wariancje średnich warunkowych odpowiedzi uczestników w warunku A vs. warunku B są nierówne?”

— Jake Westfall,

6

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:

# switch to numeric (not factor) contrast codes
d$contrast <- 2*(d$condition == 'experimental') - 1

# reduced model without correlation parameter
mod1 <- lmer(sim_1 ~ contrast + (contrast || participant_id), data=d)

# full model with correlation parameter
mod2 <- lmer(sim_1 ~ contrast + (contrast | participant_id), data=d)

# likelihood ratio test
anova(mod1, mod2)

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 punktu $X=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 punktu $X=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 stan $X=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 punktu $X=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 na $X=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

Niech $y_{ijk}$ będzie $j$ tą odpowiedzią $i$ 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

y_{ja jot k} = α_{ja} + β_{ja} x_{k} + {mi}_{ja jot k},

$y_{ijk} = \alpha_i + \beta_ix_k + e_{ijk},$ gdzie

α_{i}

$\alpha_i$ są podmiotami przechwytuje losowo i ma wariancję

σ_{α}^{2}

$\sigma^2_\alpha$ ,

β_{i}

$\beta_i$ są losowym nachyleniem badanych i mają wariancję

σ_{β}^{2}

$\sigma^2_\beta$ ,

e_{i j k}

$e_{ijk}$ jest poziomem błędu poziomu obserwacji, a

cov (α_{i}, β_{i}) = σ_{α β}

$\text{cov}(\alpha_i, \beta_i)=\sigma_{\alpha\beta}$ .

Chcemy pokazać, że

var (α_{ja} + β_{ja} x_{1}) = var (α_{ja} + β_{ja} x_{2)}) \Leftrightarrow σ_{α β} = 0.

$\text{var}(\alpha_i + \beta_ix_1) = \text{var}(\alpha_i + \beta_ix_2) \Leftrightarrow \sigma_{\alpha\beta}=0.$

Począwszy od lewej strony tego wniosku, mamy

\begin{aligned} var (α_{ja} + β_{ja} x_{1}) & = var (α_{ja} + β_{ja} x_{2)}) \\ σ_{α}^{2)} + x_{1}^{2)} σ_{β}^{2)} + 2) x_{1} σ_{α β} & = σ_{α}^{2)} + x_{2)}^{2)} σ_{β}^{2)} + 2) x_{2)} σ_{α β} \\ σ_{β}^{2)} (x_{1}^{2)} - x_{2)}^{2)}) + 2) σ_{α β} (x_{1} - x_{2)}) & = 0. \end{aligned}

$\begin{aligned} \text{var}(\alpha_i + \beta_ix_1) &= \text{var}(\alpha_i + \beta_ix_2) \\ \sigma^2_\alpha + x^2_1\sigma^2_\beta + 2x_1\sigma_{\alpha\beta} &= \sigma^2_\alpha + x^2_2\sigma^2_\beta + 2x_2\sigma_{\alpha\beta} \\ \sigma^2_\beta(x_1^2 - x_2^2) + 2\sigma_{\alpha\beta}(x_1 - x_2) &= 0. \end{aligned}$

Kody kontrastu suma do zera oznaczają, że $x_1 + x_2 = 0$ i $x_1^2 = x_2^2 = x^2$ . Następnie możemy dodatkowo zmniejszyć ostatnią linię powyższego do

\begin{aligned} σ_{β}^{2)} (x^{2)} - x^{2)}) + 2) σ_{α β} (x_{1} + x_{1}) & = 0 \\ σ_{α β} & = 0, \end{aligned}

$\begin{aligned} \sigma^2_\beta(x^2 - x^2) + 2\sigma_{\alpha\beta}(x_1 + x_1) &= 0 \\ \sigma_{\alpha\beta} &= 0, \end{aligned}$ 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 $\sigma_{\alpha\beta} = 0$ , przetestuje hipotezę zerową o równych wariancjach opisanych przez OP.

$x_1=0$ $x_2=1$

var (α_{ja}) = var (α_{ja} + β_{ja}) \Leftrightarrow σ_{α β} = - \frac{σ_{β}^{2)}}{2)} .

$\text{var}(\alpha_i) = \text{var}(\alpha_i + \beta_i) \Leftrightarrow \sigma_{\alpha\beta} = -\frac{\sigma^2_\beta}{2}.$

— Jake Westfall
źródło

To już wspaniała odpowiedź, dziękuję! Myślę, że jest to najbliższe odpowiedzi na moje pytanie, więc przyjmuję je i daję nagrodę (wkrótce wygaśnie), ale chciałbym zobaczyć algebraiczne uzasadnienie, jeśli masz na to czas i energię.

— Patrick S. Forscher,

1

@ PatrickS.Forscher Właśnie dodałem dowód

— Jake Westfall

1

a

$a$

- a

$-a$ (1 | subject)dummy

a

$a$

- a

$-a$

1

Widzę teraz, że masz rację @amoeba, dzięki za wyjaśnienie. Zmienię odpowiednio moją odpowiedź.

— Jake Westfall,

1

@amoeba Masz rację, że jest możliwe, że BLUP mogą wyjść skorelowane nawet bez parametru korelacji w modelu. Uważam jednak, że do celów testowych procedura nadal działa zgodnie z przeznaczeniem (np. Ma nominalny poziom błędu typu 1), ponieważ tylko model z parametrem korelacji jest w stanie uwzględnić to w funkcji prawdopodobieństwa i tym samym „otrzymać kredyt” za to . Oznacza to, że nawet jeśli BLUP wychodzą skorelowane w prostszym modelu, nadal jest tak, jakby efekty nie były skorelowane pod względem całkowitego prawdopodobieństwa, więc test LR będzie działał. Myślę :)

— Jake Westfall

6

Możesz przetestować istotność parametrów modelu za pomocą oszacowanych przedziałów ufności, dla których pakiet lme4 ma confint.merModfunkcję.

ładowanie (patrz na przykład przedział ufności z ładowania )

> confint(m, method="boot", nsim=500, oldNames= FALSE)
Computing bootstrap confidence intervals ...
                                                           2.5 %     97.5 %
sd_(Intercept)|participant_id                         0.32764600 0.64763277
cor_conditionexperimental.(Intercept)|participant_id -1.00000000 1.00000000
sd_conditionexperimental|participant_id               0.02249989 0.46871800
sigma                                                 0.97933979 1.08314696
(Intercept)                                          -0.29669088 0.06169473
conditionexperimental                                 0.26539992 0.60940435

profil prawdopodobieństwa (patrz na przykład jaki jest związek między prawdopodobieństwem profilu a przedziałami ufności? )

> confint(m, method="profile", oldNames= FALSE)
Computing profile confidence intervals ...
                                                          2.5 %     97.5 %
sd_(Intercept)|participant_id                         0.3490878 0.66714551
cor_conditionexperimental.(Intercept)|participant_id -1.0000000 1.00000000
sd_conditionexperimental|participant_id               0.0000000 0.49076950
sigma                                                 0.9759407 1.08217870
(Intercept)                                          -0.2999380 0.07194055
conditionexperimental                                 0.2707319 0.60727448

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 imieniem ranova. 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ć

# different model with alternative parameterization (and also correlation taken out) 
fml1 <- "~ condition + (0 + control + experimental || participant_id) "

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)}

#adding extra columns for control and experimental
d <- cbind(d,as.numeric(d$condition=='control'))
d <- cbind(d,1-as.numeric(d$condition=='control'))
names(d)[c(4,5)] <- c("control","experimental")

Teraz lmerpoda wariancję dla różnych grup

> m <- lmer(paste("sim_1 ", fml1), data=d)
> m
Linear mixed model fit by REML ['lmerModLmerTest']
Formula: paste("sim_1 ", fml1)
   Data: d
REML criterion at convergence: 2408.186
Random effects:
 Groups           Name         Std.Dev.
 participant_id   control      0.4963  
 participant_id.1 experimental 0.4554  
 Residual                      1.0268  
Number of obs: 800, groups:  participant_id, 40
Fixed Effects:
          (Intercept)  conditionexperimental  
               -0.114                  0.439

Możesz zastosować do nich metody profilowania. Na przykład teraz confint daje przedziały ufności dla wariancji kontrolnej i ćwiczeniowej.

> confint(m, method="profile", oldNames= FALSE)
Computing profile confidence intervals ...
                                    2.5 %     97.5 %
sd_control|participant_id       0.3490873 0.66714568
sd_experimental|participant_id  0.3106425 0.61975534
sigma                           0.9759407 1.08217872
(Intercept)                    -0.2999382 0.07194076
conditionexperimental           0.1865125 0.69149396

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.

$j$

{\hat{Y}}_{ja, jot} \sim N. (μ_{jot}, σ_{jot}^{2)} + \frac{σ_{ϵ}^{2)}}{10})

$\hat{Y}_{i,j} \sim N(\mu_j, \sigma_j^2 + \frac{\sigma_{\epsilon}^2}{10})$

$\sigma_\epsilon$ $\sigma_{j}$ $j = \lbrace 1,2 \rbrace$

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.

^{$\sigma_{j=1} = \sigma_{j=2} = 0.5$ $\sigma_\epsilon=1$}

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

$\sigma_{j=1} \neq \sigma_{j=2}$ $\hat{Y}_{i,j}$ $\sigma_j$ $\sigma_\epsilon$

^{$\sigma_{j=1} = 0.5$ $\sigma_{j=2} = 0.25$ $\sigma_\epsilon=1$}

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:

set.seed(23432)

# different model with alternative parameterization (and also correlation taken out)
fml1 <- "~ condition + (0 + control + experimental || participant_id) "
fml <- "~ condition + (condition | participant_id)"

n <- 10000

theta_m <- matrix(rep(0,n*2),n)
theta_f <- matrix(rep(0,n*2),n)

# initial data frame later changed into d by adding a sixth sim_1 column
ds <- expand.grid(participant_id=1:40, trial_num=1:10)
ds <- rbind(cbind(ds, condition="control"), cbind(ds, condition="experimental"))
  #adding extra columns for control and experimental
  ds <- cbind(ds,as.numeric(ds$condition=='control'))
  ds <- cbind(ds,1-as.numeric(ds$condition=='control'))
  names(ds)[c(4,5)] <- c("control","experimental")

# defining variances for the population of individual means
stdevs <- c(0.5,0.5) # c(control,experimental)

pb <- txtProgressBar(title = "progress bar", min = 0,
                    max = n, style=3)
for (i in 1:n) {

  indv_means <- c(rep(0,40)+rnorm(40,0,stdevs[1]),rep(0.5,40)+rnorm(40,0,stdevs[2]))
  fill <- indv_means[d[,1]+d[,5]*40]+rnorm(80*10,0,sqrt(1)) #using a different way to make the data because the simulate is not creating independent data in the two groups 
  #fill <- suppressMessages(simulate(formula(fml), 
  #                     newparams=list(beta=c(0, .5), 
  #                                    theta=c(.5, 0, 0), 
  #                                    sigma=1), 
  #                     family=gaussian, 
  #                     newdata=ds))
  d <- cbind(ds, fill)
  names(d)[6] <- c("sim_1")


  m <- lmer(paste("sim_1 ", fml1), data=d)
  m
  theta_m[i,] <- m@theta^2

  imeans <- aggregate(d[, 6], list(d[,c(1)],d[,c(3)]), mean)
  theta_f[i,1] <- var(imeans[c(1:40),3])
  theta_f[i,2] <- var(imeans[c(41:80),3])

  setTxtProgressBar(pb, i)
}
close(pb)

p1 <- hist(theta_f[,1]/theta_f[,2], breaks = seq(0,6,0.06))       
fr <- theta_m[,1]/theta_m[,2]
fr <- fr[which(fr<30)]
p2 <- hist(fr, breaks = seq(0,30,0.06))



plot(-100,-100, xlim=c(0,6), ylim=c(0,800), 
     xlab="F-score", ylab = "counts [n out of 10 000]")
plot( p1, col=rgb(0,0,1,1/4), xlim=c(0,6), ylim=c(0,800), add=T)  # means based F-score
plot( p2, col=rgb(1,0,0,1/4), xlim=c(0,6), ylim=c(0,800), add=T)  # model based F-score
fr <- seq(0, 4, 0.01)
lines(fr,df(fr,39,39)*n*0.06,col=1)
legend(2, 800, c("means based F-score","mixed regression based F-score"), 
       fill=c(rgb(0,0,1,1/4),rgb(1,0,0,1/4)),box.col =NA, bg = NA)
legend(2, 760, c("F(39,39) distribution"), 
       lty=c(1),box.col = NA,bg = NA)
title(expression(paste(sigma[1]==0.5, " , ", sigma[2]==0.5, " and ", sigma[epsilon]==1)))

— Sextus Empiricus
źródło

Jest to przydatne, ale nie wydaje się odpowiadać na pytanie, jak porównywać wariancje w dwóch warunkach.

— ameba mówi Przywróć Monikę

@amoeba Odkryłem, że ta odpowiedź stanowi sedno problemu (o testowaniu losowych składników wariancji). To, czego dokładnie chce OP, jest trudne do odczytania w całym tekście. Do czego odnosi się „wariancje losowego przechwytywania”? (liczba mnoga w odniesieniu do przechwytywania mnie myli) Jednym z możliwych przypadków może być użycie modelu, 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).

— Sextus Empiricus

Każdy osobnik ma pewną średnią odpowiedź w warunku A i pewną średnią odpowiedź w warunku B. Pytanie brzmi, czy wariancja między osobnikami w A jest różna od wariancji między osobnikami w B.

— mówi ameba Przywróć Monikę

To nie kończy zadania postawionego w tytule „Porównaj losowy komponent wariancji między poziomami zmiennej grupującej”. Zauważyłem, że w treści pytania znalazła się myląca literówka, którą naprawiłem. Starałem się również doprecyzować treść pytania.

— Patrick S. Forscher,

Odpowiedź na pytanie może być możliwa 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.

— Patrick S. Forscher,

5

Jednym stosunkowo prostym sposobem może być zastosowanie testów współczynnika wiarygodności, anovazgodnie z opisem w lme4FAQ .

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 = FALSEmimo REML = TRUEze anova(..., refit = FALSE)jest całkowicie wykonalne ).

m_full <- lmer(sim_1 ~ condition + (condition | participant_id), data=d, REML = FALSE)
summary(m_full)$varcor
 # Groups         Name                  Std.Dev. Corr  
 # participant_id (Intercept)           0.48741        
 #                conditionexperimental 0.26468  -0.419
 # Residual                             1.02677     

m_red <- lmer(sim_1 ~ condition + (1 | participant_id), data=d, REML = FALSE)
summary(m_red)$varcor
 # Groups         Name        Std.Dev.
 # participant_id (Intercept) 0.44734 
 # Residual                   1.03571 

anova(m_full, m_red)
# Data: d
# Models:
# m_red: sim_1 ~ condition + (1 | participant_id)
# m_full: sim_1 ~ condition + (condition | participant_id)
#        Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
# m_red   4 2396.6 2415.3 -1194.3   2388.6                         
# m_full  6 2398.7 2426.8 -1193.3   2386.7 1.9037      2      0.386

Ten test jest jednak prawdopodobnie konserwatywny . Na przykład FAQ mówi:

Należy pamiętać, że testy hipotezy zerowej oparte na LRT są zachowawcze, gdy wartość zerowa (np. Σ2 = 0) znajduje się na granicy możliwej przestrzeni; w najprostszym przypadku (wariancja pojedynczego efektu losowego) wartość p jest w przybliżeniu dwa razy większa niż powinna (Pinheiro i Bates 2000).

Istnieje kilka alternatyw:

$\chi^2$
Symuluj prawidłową dystrybucję za pomocą RLRsim(jak opisano również w FAQ).

Przedstawię drugą opcję w następujący sposób:

library("RLRsim")
## reparametrize model so we can get one parameter that we want to be zero:
afex::set_sum_contrasts() ## warning, changes contrasts globally
d <- cbind(d, difference = model.matrix(~condition, d)[,"condition1"])

m_full2 <- lmer(sim_1 ~ condition + (difference | participant_id), data=d, REML = FALSE)
all.equal(deviance(m_full), deviance(m_full2))  ## both full models are identical

## however, we need the full model without correlation!
m_full2b <- lmer(sim_1 ~ condition + (1| participant_id) + 
                   (0 + difference | participant_id), data=d, REML = FALSE)
summary(m_full2b)$varcor
 # Groups           Name        Std.Dev.
 # participant_id   (Intercept) 0.44837 
 # participant_id.1 difference  0.13234 
 # Residual                     1.02677 

## model that only has random effect to be tested
m_red <- update(m_full2b,  . ~ . - (1 | participant_id), data=d, REML = FALSE)
summary(m_red)$varcor
 # Groups         Name       Std.Dev.
 # participant_id difference 0.083262
 # Residual                  1.125116

## Null model 
m_null <- update(m_full2b,  . ~ . - (0 + difference | participant_id), data=d, REML = FALSE)
summary(m_null)$varcor
 # Groups         Name        Std.Dev.
 # participant_id (Intercept) 0.44734 
 # Residual                   1.03571 

exactRLRT(m_red, m_full2b, m_null)
# Using restricted likelihood evaluated at ML estimators.
# Refit with method="REML" for exact results.
# 
#   simulated finite sample distribution of RLRT.
#   
#   (p-value based on 10000 simulated values)
# 
# data:  
# RLRT = 1.9698, p-value = 0.0719

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 RLRsimpozwala na jednoczesne testowanie wielu składników, ale jeśli tak, można to zrobić w ten sam sposób.

Odpowiedź na komentarz:

$\theta_X$ $\theta_0$ $X$

Nie jestem pewien, czy na to pytanie można uzyskać rozsądną odpowiedź.

Losowe przechwytywanie pozwala na idiosynkratyczną różnicę w poziomie ogólnym dla każdego poziomu czynnika grupującego. Na przykład, jeśli zmienną zależną jest czas odpowiedzi, niektórzy uczestnicy są szybsi, a niektórzy wolniej.
Losowe nachylenie umożliwia każdemu poziomowi czynnika grupującego idiosynkratyczny wpływ czynnika, dla którego szacowane są losowe nachylenia. Na przykład, jeśli czynnikiem jest zgodność, to niektórzy uczestnicy mogą mieć większy efekt zgodności niż inni.

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.

— Henrik
źródło

1

Używasz kontrastów leczenia ( 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/…

— Henrik

1

Dzięki! Moje pytania są jednak nieco inne: (1) Twoja odpowiedź wydaje się sugerować, że moje pytanie ogranicza się do „to losowe nachylenie dla warunku różnego od 0”. Czy to prawda? (2) Jeśli odpowiedź na (1) brzmi „tak”, sugeruje to inną interpretację losowego nachylenia dla condition: umożliwia losowe przechwytywanie zmienianie się na różnych poziomach condition. Czy to prawda?

— Patrick S. Forscher,

2

Mój 2 ¢: kontrprzykład @amoeba na proponowaną przez Henrika procedurę jest poprawny. Henrik ma prawie rację, ale porównuje niewłaściwą parę modeli. Porównanie modeli, które odpowiada na pytanie Patricka, to porównanie modeli, które Henrik nazwał m_fullkontra m_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ć ...

— Jake Westfall,

2

To nie jest właściwe wyjaśnienie, ale przestudiowanie mojej odpowiedzi tutaj może rzucić nieco światła na tę sprawę. Zasadniczo parametr korelacji kontroluje, czy linie regresji uczestnika „rozsuwają się w prawo” (korekcja dodatnia), czy „rozkładają się w lewo” (korekta ujemna). Każda z tych implikuje nierówną wariancję średnich warunkowych odpowiedzi uczestników. Kodowanie sumy do zera gwarantuje, że szukamy korelacji we właściwym punkcie na X

— Jake Westfall

2

Rozważę opublikowanie odpowiedzi ze zdjęciami, jeśli znajdę czas ...

— Jake Westfall

5

Twój model

m = lmer(sim_1 ~ condition + (condition | participant_id), data=d)

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

m = lmer(sim_1 ~ 0 + condition + (0 + condition | participant_id), data=d)

Macierz kowariancji losowej ma teraz prostszą interpretację:

Random effects:
 Groups         Name                  Variance Std.Dev. Corr
 participant_id conditioncontrol      0.2464   0.4963       
                conditionexperimental 0.2074   0.4554   0.83

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

delta = as.data.frame(VarCorr(m))[1,4] - as.data.frame(VarCorr(m))[2,4]

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

set.seed(42)
nrep = 100
v = matrix(nrow=nrep, ncol=1)
for (i in 1:nrep)
{
   dp = d
   for (s in unique(d$participant_id)){             
     if (rbinom(1,1,.5)==1){
       dp[p$participant_id==s & d$condition=='control',]$condition = 'experimental'
       dp[p$participant_id==s & d$condition=='experimental',]$condition = 'control'
     }
   }
  m <- lmer(sim_1 ~ 0 + condition + (0 + condition | participant_id), data=dp)
  v[i,] = as.data.frame(VarCorr(m))[1,4] - as.data.frame(VarCorr(m))[2,4]
}
pvalue = sum(abs(v) >= abs(delta)) / nrep

$p=0.7$ nrep

Dokładnie taką samą logikę można zastosować w przypadku bonusu.

— ameba mówi Przywróć Monikę
źródło

Super ciekawe, dziękuję! Muszę się zastanowić, dlaczego działa ta zmiana parametrów, ponieważ wydaje się, że jest to kluczowy wgląd w tę odpowiedź.

— Patrick S. Forscher,

O dziwo, wartości przechwytywania dla grupy w twojej odpowiedzi wydają się różnić od tych w odpowiedzi @MartijnWeterings.

— Patrick S. Forscher,

@ PatrickS.Forscher To dlatego, że myślę, że on generuje inny zestaw danych. Mogę użyć sim_1 ~ 0 + condition + (0 + dummy(condition, "control") + dummy(condition, "experimental") | participant_id)sformułowania i uzyskać taki sam wynik jak w mojej odpowiedzi.

— ameba mówi Przywróć Monikę

1

@ PatrickS.Forscher Nie, użyłem danych wygenerowanych przez twój kod (z twoim ziarnem). Ustawiłem ziarno na 42 tylko podczas przeprowadzania testu permutacji. To Martijn zmienił zestaw danych, a nie ja.

— Ameba mówi Przywróć Monikę

1

Ta propozycja jest zdecydowanie solidna. Jak myślę, że już doświadczyłeś, konfigurowanie testów permutacyjnych dla danych wielopoziomowych nie jest całkiem proste. Podobne podejście, które byłoby nieco łatwiejsze do wdrożenia, to ładowanie parametryczne, które jest dość proste w przypadku lme4 przy użyciu metody simulate () dopasowanych obiektów lmer, tzn. Wywołuje symulację (m) wiele razy, aby zbudować bootstrap dystrybucja. Pomysł do zabawy.

— Jake Westfall,

0

y_{ja jot k} = μ + α_{jot} + {re}_{ja jot} + {mi}_{ja jot k}, {re}_{ja} \sim N. (0, Σ), {mi}_{ja jot k} \sim N. (0, σ^{2)})

$y_{ijk} = \mu + \alpha_j + d_{ij} + e_{ijk}, \quad d_{i} \sim N(0, \Sigma), \quad e_{ijk} \sim N(0, \sigma^2)$

α_{j}

$\alpha_j$

j

$j$

d_{i} = (d_{i 1}, \dots, d_{i J})^{⊤}

$d_i = (d_{i1}, \ldots, d_{iJ})^\top$

i

$i$

j

$j$

y_{i 1 k}

$y_{i1k}$

y_{i 2 k}

$y_{i2k}$

A

$A$

B

$B$

d_{i}

$d_{i}$

Σ = [\begin{matrix} σ_{ZA}^{2)} & σ_{ZA b} \\ σ_{ZA b} & σ_{b}^{2)} \end{matrix}]

$\Sigma = \begin{bmatrix}\sigma^2_A & \sigma_{AB}\\ \sigma_{AB} & \sigma^2_{B} \end{bmatrix}$

$\sigma^2_A$ $\sigma^2_{B}$

$\Sigma$

Odmiana przecięcia, która odpowiada wielkiemu środkowi, jest

σ_{1}^{2)} : = Var (średnia średnia) = Var (\frac{1}{2)} \cdot (ZA + b)) = \frac{1}{4} \cdot (Var (ZA) + Var (b) + 2) \cdot Cov (ZA, b)) .

$\sigma^2_{1} := \text{Var(grand mean)} = \text{Var}(\frac{1}{2} \cdot (A+B)) = \frac{1}{4} \cdot (\text{Var}(A) + \text{Var}(B) + 2 \cdot \text{Cov}(A,B)).$

Wariant kontrastu jest

σ_{2)}^{2)} : = Var (kontrast) = Var (\frac{1}{2)} \cdot (ZA - b)) = \frac{1}{4} \cdot (Var (ZA) + Var (b) - 2) \cdot Cov (ZA, b)) .

$\sigma^2_{2} := \text{Var(contrast)} = \text{Var}(\frac{1}{2} \cdot (A-B)) = \frac{1}{4} \cdot (\text{Var}(A) + \text{Var}(B) - 2 \cdot \text{Cov}(A,B)).$

A kowariancja między przechwyceniem a kontrastem jest

σ_{12} : = Cov (wielki środek, kontrast) = Cov (\frac{1}{2)} \cdot (ZA + b), \frac{1}{2)} \cdot (ZA - b)) = \frac{1}{4} \cdot (Var (ZA) - Var (b)) .

$\sigma_{12} := \text{Cov}(\text{grand mean, contrast}) = \text{Cov}(\frac{1}{2} \cdot (A+B), \frac{1}{2} \cdot (A-B)) \\ = \frac{1}{4} \cdot (\text{Var}(A) - \text{Var}(B)).$

$\Sigma$

Σ = [\begin{matrix} σ_{1}^{2)} + σ_{2)}^{2)} + 2) σ_{12} & σ_{1}^{2)} - σ_{2)}^{2)} \\ σ_{1}^{2)} - σ_{2)}^{2)} & σ_{1}^{2)} + σ_{2)}^{2)} - 2) σ_{12} \end{matrix}] = [\begin{matrix} σ_{ZA}^{2)} & σ_{ZA b} \\ σ_{ZA b} & σ_{b}^{2)} \end{matrix}] .

$\Sigma = \begin{bmatrix}\sigma^2_1 + \sigma^2_2 + 2\sigma_{12} & \sigma^2_1 - \sigma^2_2\\ \sigma^2_1 - \sigma^2_2 & \sigma^2_1 + \sigma^2_2 - 2\sigma_{12} \end{bmatrix} = \begin{bmatrix}\sigma^2_A & \sigma_{AB}\\ \sigma_{AB} & \sigma^2_{B} \end{bmatrix}.$

$\Sigma$

Σ = [\begin{matrix} σ_{1}^{2)} & σ_{1}^{2)} \\ σ_{1}^{2)} & σ_{1}^{2)} \end{matrix}] + [\begin{matrix} σ_{2)}^{2)} & - σ_{2)}^{2)} \\ - σ_{2)}^{2)} & σ_{2)}^{2)} \end{matrix}] + 2) [\begin{matrix} σ_{12} & 0 \\ 0 & - σ_{12} \end{matrix}] .

$\Sigma = \begin{bmatrix}\sigma^2_1 & \sigma^2_1\\ \sigma^2_1 & \sigma^2_1\end{bmatrix} + \begin{bmatrix}\sigma^2_2 & -\sigma^2_2\\ -\sigma^2_2 & \sigma^2_2\end{bmatrix} + 2\begin{bmatrix}\sigma_{12} & 0\\ 0 & -\sigma_{12}\end{bmatrix}.$

$\sigma_{12}$

Σ = [\begin{matrix} σ_{1}^{2)} & σ_{1}^{2)} \\ σ_{1}^{2)} & σ_{1}^{2)} \end{matrix}] + [\begin{matrix} σ_{2)}^{2)} & - σ_{2)}^{2)} \\ - σ_{2)}^{2)} & σ_{2)}^{2)} \end{matrix}] = [\begin{matrix} σ_{1}^{2)} + σ_{2)}^{2)} & σ_{1}^{2)} - σ_{2)}^{2)} \\ σ_{1}^{2)} - σ_{2)}^{2)} & σ_{1}^{2)} + σ_{2)}^{2)} \end{matrix}]

$\Sigma = \begin{bmatrix}\sigma^2_1 & \sigma^2_1\\ \sigma^2_1 & \sigma^2_1\end{bmatrix} + \begin{bmatrix}\sigma^2_2 & -\sigma^2_2\\ -\sigma^2_2 & \sigma^2_2\end{bmatrix} = \begin{bmatrix}\sigma^2_1 + \sigma^2_2 & \sigma^2_1 - \sigma^2_2 \\ \sigma^2_1 - \sigma^2_2 & \sigma^2_1 + \sigma^2_2\end{bmatrix}$

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 = FALSEgdy używasz oszacowania REML), gdzie mod1i mod2są zdefiniowane tak jak @Jake Westfall.

$\sigma_{12}$ $\sigma^2_2$

Σ = [\begin{matrix} σ_{1}^{2)} & σ_{1}^{2)} \\ σ_{1}^{2)} & σ_{1}^{2)} \end{matrix}]

$\Sigma = \begin{bmatrix}\sigma^2_1 & \sigma^2_1\\ \sigma^2_1 & \sigma^2_1\end{bmatrix}$

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

# code snippet from Jake Westfall
d$contrast <- 2*(d$condition == 'experimental') - 1

# new model
mod3 <- lmer(sim_1 ~ contrast + (1 | participant_id) + (1 | contrast:participant_id), 
             data = d, REML = FALSE)

lub (przy użyciu zmiennej / współczynnika kategorialnego condition)

mod4 <- lmer(sim_1 ~ condition + (1 | participant_id) + (1 | condition:participant_id), 
             data = d, REML = FALSE)

z

Σ = [\begin{matrix} σ_{1}^{2)} + σ_{2)}^{2)} & σ_{1}^{2)} \\ σ_{1}^{2)} & σ_{1}^{2)} + σ_{2)}^{2)} \end{matrix}] = [\begin{matrix} σ_{1}^{2)} & σ_{1}^{2)} \\ σ_{1}^{2)} & σ_{1}^{2)} \end{matrix}] + [\begin{matrix} σ_{2)}^{2)} & 0 \\ 0 & σ_{2)}^{2)} \end{matrix}]

$\Sigma = \begin{bmatrix}\sigma^2_1 + \sigma^2_2 & \sigma^2_1\\ \sigma^2_1 & \sigma^2_1 + \sigma^2_2 \end{bmatrix} = \begin{bmatrix}\sigma^2_1 & \sigma^2_1\\ \sigma^2_1 & \sigma^2_1\end{bmatrix} + \begin{bmatrix}\sigma^2_2 & 0\\ 0 & \sigma^2_2\end{bmatrix}$

$\sigma^2_1$ $\sigma^2_2$ $\Sigma$

Poniżej widzimy, że mod1, mod3i mod4otrzymując równoważne pasuje:

# code snippet from Jake Westfall
d$contrast <- 2*(d$condition == 'experimental') - 1

mod1 <- lmer(sim_1 ~ contrast + (contrast || participant_id),
             data = d, REML = FALSE)

mod2 <- lmer(sim_1 ~ contrast + (contrast | participant_id),
             data = d, REML = FALSE)

# new models 
mod3 <- lmer(sim_1 ~ contrast + (1 | participant_id) + (1 | contrast:participant_id), 
             data = d, REML = FALSE) 

mod4 <- lmer(sim_1 ~ condition + (1 | participant_id) + (1 | condition:participant_id), 
             data = d, REML = FALSE)

anova(mod3, mod1)
# Data: d
# Models:
# mod3: sim_1 ~ contrast + (1 | participant_id) + (1 | contrast:participant_id)
# mod1: sim_1 ~ contrast + ((1 | participant_id) + (0 + contrast | participant_id))
#      Df    AIC    BIC  logLik deviance Chisq Chi Df Pr(>Chisq)
# mod3  5 2396.9 2420.3 -1193.5   2386.9                        
# mod1  5 2396.9 2420.3 -1193.5   2386.9     0      0          1

anova(mod4, mod3)
# Data: d
# Models:
# mod4: sim_1 ~ condition + (1 | participant_id) + (1 | condition:participant_id)
# mod3: sim_1 ~ contrast + (1 | participant_id) + (1 | contrast:participant_id)
#      Df    AIC    BIC  logLik deviance Chisq Chi Df Pr(>Chisq)
# mod4  5 2396.9 2420.3 -1193.5   2386.9                        
# mod3  5 2396.9 2420.3 -1193.5   2386.9     0      0          1

$\Sigma$

Σ = [\begin{matrix} σ_{1}^{2)} & σ_{1}^{2)} + σ_{12} \\ σ_{1}^{2)} + σ_{12} & σ_{1}^{2)} + σ_{2)}^{2)} + 2) σ_{12} \end{matrix}] = [\begin{matrix} σ_{1}^{2)} & σ_{1}^{2)} \\ σ_{1}^{2)} & σ_{1}^{2)} \end{matrix}] + [\begin{matrix} 0 & 0 \\ 0 & σ_{2)}^{2)} \end{matrix}] + [\begin{matrix} 0 & σ_{12} \\ σ_{12} & 2) σ_{12} \end{matrix}]

$\Sigma = \begin{bmatrix}\sigma^2_1 & \sigma^2_1 + \sigma_{12}\\ \sigma^2_1 + \sigma_{12} & \sigma^2_1 + \sigma^2_2 + 2\sigma_{12} \end{bmatrix} = \begin{bmatrix}\sigma^2_1 & \sigma^2_1\\ \sigma^2_1 & \sigma^2_1\end{bmatrix} + \begin{bmatrix}0 & 0\\ 0 & \sigma^2_2\end{bmatrix} + \begin{bmatrix}0 & \sigma_{12}\\ \sigma_{12} & 2\sigma_{12}\end{bmatrix}$

$\sigma^2_1$ $A$ $\sigma^2_2$ $A - B$ $\sigma_{12}$

$\sigma_{12}$ $\sigma^2_2$

mod4 $\Sigma$

— statmerkur
źródło

Model efektów mieszanych: Porównaj losową składową wariancji między poziomami zmiennej grupującej

Krótka odpowiedź

Wyjaśnienie wizualne / intuicja

Dowód