Dlaczego otrzymuję zerową wariancję efektu losowego w moim modelu mieszanym, pomimo pewnych różnic w danych?

Przeprowadziliśmy regresję logistyczną efektów mieszanych przy użyciu następującej składni;

# fit model
fm0 <- glmer(GoalEncoding ~ 1 + Group + (1|Subject) + (1|Item), exp0,
             family = binomial(link="logit"))
# model output
summary(fm0)

Przedmiot i przedmiot to efekty losowe. Otrzymujemy nieparzysty wynik, który jest współczynnikiem, a odchylenie standardowe dla przedmiotowego terminu wynosi zero;

Generalized linear mixed model fit by maximum likelihood (Laplace
Approximation) [glmerMod]
Family: binomial  ( logit )
Formula: GoalEncoding ~ 1 + Group + (1 | Subject) + (1 | Item)
Data: exp0

AIC      BIC      logLik deviance df.resid 
449.8    465.3   -220.9    441.8      356 

Scaled residuals: 
Min     1Q Median     3Q    Max 
-2.115 -0.785 -0.376  0.805  2.663 

Random effects:
Groups  Name        Variance Std.Dev.
Subject (Intercept) 0.000    0.000   
Item    (Intercept) 0.801    0.895   
Number of obs: 360, groups:  Subject, 30; Item, 12

Fixed effects:
                Estimate Std. Error z value Pr(>|z|)    
 (Intercept)     -0.0275     0.2843    -0.1     0.92    
 GroupGeMo.EnMo   1.2060     0.2411     5.0  5.7e-07 ***
 ---
 Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

 Correlation of Fixed Effects:
             (Intr)
 GroupGM.EnM -0.002

To nie powinno się zdarzyć, ponieważ oczywiście istnieją różnice między podmiotami. Kiedy przeprowadzamy tę samą analizę w stacie

xtmelogit goal group_num || _all:R.subject || _all:R.item

Note: factor variables specified; option laplace assumed

Refining starting values: 

Iteration 0:   log likelihood = -260.60631  
Iteration 1:   log likelihood = -252.13724  
Iteration 2:   log likelihood = -249.87663  

Performing gradient-based optimization: 

Iteration 0:   log likelihood = -249.87663  
Iteration 1:   log likelihood = -246.38421  
Iteration 2:   log likelihood =  -245.2231  
Iteration 3:   log likelihood = -240.28537  
Iteration 4:   log likelihood = -238.67047  
Iteration 5:   log likelihood = -238.65943  
Iteration 6:   log likelihood = -238.65942  

Mixed-effects logistic regression               Number of obs      =       450
Group variable: _all                            Number of groups   =         1

                                                Obs per group: min =       450
                                                               avg =     450.0
                                                               max =       450

Integration points =   1                        Wald chi2(1)       =     22.62
Log likelihood = -238.65942                     Prob > chi2        =    0.0000

------------------------------------------------------------------------------
        goal |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
   group_num |   1.186594    .249484     4.76   0.000     .6976147    1.675574
       _cons |  -3.419815   .8008212    -4.27   0.000    -4.989396   -1.850234
------------------------------------------------------------------------------

------------------------------------------------------------------------------
  Random-effects Parameters  |   Estimate   Std. Err.     [95% Conf. Interval]
-----------------------------+------------------------------------------------
_all: Identity               |
               sd(R.subject) |   7.18e-07   .3783434             0           .
-----------------------------+------------------------------------------------
_all: Identity               |
                 sd(R.trial) |   2.462568   .6226966      1.500201    4.042286
------------------------------------------------------------------------------
LR test vs. logistic regression:     chi2(2) =   126.75   Prob > chi2 = 0.0000

Note: LR test is conservative and provided only for reference.
Note: log-likelihood calculations are based on the Laplacian approximation.

wyniki są zgodne z oczekiwaniami przy niezerowym współczynniku / se dla terminu podmiotu.

Początkowo sądziliśmy, że może to mieć coś wspólnego z kodowaniem terminu podmiotu, ale zmiana tego ciągu z ciągu na liczbę całkowitą nie zrobiła żadnej różnicy.

Oczywiście analiza nie działa poprawnie, ale nie jesteśmy w stanie określić źródła trudności. (Uwaga: ktoś inny na tym forum ma podobny problem, ale ten wątek pozostaje bez odpowiedzi link do pytania )

Omówiono to nieco na stronie https://bbolker.github.io/mixedmodels-misc/glmmFAQ.html (wyszukaj „pojedyncze modele”); jest to powszechne, zwłaszcza gdy istnieje niewielka liczba grup (chociaż 30 nie jest szczególnie małe w tym kontekście).

Jedną różnicą między lme4wieloma innymi pakietami jest to, że wiele pakietów, w tym lme4poprzednik nlme, radzi sobie z faktem, że oszacowania wariancji muszą być nieujemne poprzez dopasowanie wariancji w skali logarytmicznej: oznacza to, że szacunki wariancji nie mogą być dokładnie zerowe, tylko bardzo bardzo mały. lme4, natomiast wykorzystuje ograniczoną optymalizację, dzięki czemu może zwracać wartości dokładnie zerowe ( więcej dyskusji można znaleźć na stronie http://arxiv.org/abs/1406.5823 str. 24). http://rpubs.com/bbolker/6226 daje przykład.

W szczególności, patrząc uważnie na wyniki wariancji między pacjentami ze Staty, masz szacunkową wartość 7,18e-07 (w stosunku do przecięcia -3,4) z odchyleniem standardowym Wald wynoszącym 0,3783434 (w tym przypadku zasadniczo bezużytecznym!) I 95% CI wymienione jako „0”; jest to technicznie „niezerowe”, ale jest tak bliskie zeru, jak program zgłosi ...

Jest dobrze znane i teoretycznie możliwe do udowodnienia (np. Stram i Lee Biometrics 1994), że rozkład zerowy dla składników wariancyjnych jest mieszaniną masy punktowej („spike”) przy zerowym i rozkładem chi-kwadrat od zera. Nic dziwnego (ale nie wiem, czy jest to udowodnione / dobrze znane), rozkład próbkowania oszacowań komponentu wariancji często ma skok zerowy, nawet gdy prawdziwa wartość nie jest zerowa - patrz np. Http://rpubs.com/ bbolker / 4187 na przykład lub ostatni przykład na ?bootMerstronie:

library(lme4)
library(boot)
## Check stored values from a longer (1000-replicate) run:
load(system.file("testdata","boo01L.RData",package="lme4"))
plot(boo01L,index=3) 

Nie sądzę, że jest problem. Lekcja na podstawie wyników modelu jest taka, że ​​chociaż istnieje „oczywista” zmienność wyników podmiotu, zakres tej zmienności podmiotu można w pełni lub praktycznie całkowicie wyjaśnić jedynie przez sam warunek rezydualny wariancji. Nie ma wystarczającej dodatkowej zmienności na poziomie podmiotu, aby uzasadnić dodanie dodatkowego losowego efektu na poziomie podmiotu w celu wyjaśnienia wszystkich zaobserwowanych zmian.

Pomyśl o tym w ten sposób. Wyobraź sobie, że symulujemy dane eksperymentalne w ramach tego samego paradygmatu. Ustawiamy parametry tak, aby istniała zmienność resztkowa na zasadzie prób po próbie, ale 0 zmienność na poziomie podmiotu (tj. Wszyscy badani mają tę samą „prawdziwą średnią” plus błąd). Teraz za każdym razem, gdy symulujemy dane z tego zestawu parametrów, oczywiście stwierdzimy, że badani nie mają dokładnie takiej samej wydajności. Niektóre kończą się niskimi wynikami, niektóre z wysokimi. Ale wszystko to tylko ze względu na resztkową zmienność na poziomie próby. „Wiemy” (dzięki określeniu parametrów symulacji), że tak naprawdę nie ma żadnej zmiany na poziomie podmiotu.