uzyskiwanie stopni swobody od lmer

11

Dopasowałem lmer model z następującymi (aczkolwiek wymyślonymi danymi wyjściowymi):

Random effects:
 Groups        Name        Std.Dev.
 day:sample (Intercept)    0.09
 sample        (Intercept) 0.42
 Residual                  0.023 

Naprawdę chciałbym zbudować przedział ufności dla każdego efektu, używając następującej formuły:

(n1)s2χα/2,n12,(n1)s2χ1α/2,n12

Czy istnieje sposób, aby wygodnie wydostać się ze stopni swobody?

użytkownik1357015
źródło
1
Myślę, że chcesz sprawdzić lmerTest . Istnieje wiele przybliżeń przybliżenia df w modelu efektów mieszanych dla efektów stałych (np. Satterthwaite , Kenward-Roger itp.). W twoim przypadku wydaje mi się, że nadmiernie komplikujesz swoje życie. Zakładasz, że każdy efekt jest gaussowski. Wystarczy użyć odchylenia standardowego, aby uzyskać wybrany przedział ufności.
usεr11852
3
@ usεr11852 W modelu z efektami mieszanymi zakłada się, że każdy efekt jest gaussowski, ale parametrem jest wariancja rozkładu Gaussa, a nie średnia. Tak więc rozkład jego estymatora będzie bardzo wypaczony, a normalny przedział ufności odchyleń standardowych ± 2 nie będzie odpowiedni.
Karl Ove Hufthammer 18.04.15
1
@KarlOveHufthammer: Masz rację, zwracając na to uwagę; Rozumiem, co masz na myśli (i prawdopodobnie OP). Myślałem, że był zaniepokojony środkami i / lub realizacją efektów losowych, kiedy wspominał o stopniach swobody.
usεr11852 18.04.15
stopnie swobody są „problematyczne” dla modeli mieszanych, patrz: stat.ethz.ch/pipermail/r-help/2006- maj 094765.html i stats.stackexchange.com/questions/84268/…
Tim

Odpowiedzi:

17

Zamiast tego po prostu utworzę przedziały ufności prawdopodobieństwa profilu . Są niezawodne i bardzo łatwe do obliczenia przy użyciu pakietu „lme4”. Przykład:

> library(lme4)
> fm = lmer(Reaction ~ Days + (Days | Subject),
            data=sleepstudy)
> summary(fm)
[]
Random effects:
 Groups   Name        Variance Std.Dev. Corr
 Subject  (Intercept) 612.09   24.740       
          Days         35.07    5.922   0.07
 Residual             654.94   25.592       

Możesz teraz obliczyć przedziały ufności prawdopodobieństwa profilu za pomocą confint()funkcji:

> confint(fm, oldNames=FALSE)
Computing profile confidence intervals ...
                               2.5 %  97.5 %
sd_(Intercept)|Subject        14.381  37.716
cor_Days.(Intercept)|Subject  -0.482   0.685
sd_Days|Subject                3.801   8.753
sigma                         22.898  28.858
(Intercept)                  237.681 265.130
Days                           7.359  13.576

Możesz również użyć parametrycznego ładowania początkowego, aby obliczyć przedziały ufności. Oto składnia R (używając parmargumentu, aby ograniczyć parametry, dla których chcemy mieć przedziały ufności):

> confint(fm, method="boot", nsim=1000, parm=1:3)
Computing bootstrap confidence intervals ...
                              2.5 % 97.5 %
sd_(Intercept)|Subject       11.886 35.390
cor_Days.(Intercept)|Subject -0.504  0.929
sd_Days|Subject               3.347  8.283

Wyniki będą się naturalnie różnić nieco dla każdego przebiegu. Możesz zwiększyć, nsimaby zmniejszyć tę zmienność, ale zwiększy to również czas potrzebny do oszacowania przedziałów ufności.

Karl Ove Hufthammer
źródło
1
Dobra odpowiedź (+1). Chciałbym również wspomnieć o tym, że w tym przypadku można również uzyskać CI z parametrycznego bootstrapu . Ten wątek zawiera bardzo interesującą dyskusję na ten temat.
usεr11852 18.04.15
@ usεr11852 Dzięki za sugestię. Dodałem teraz przykład wykorzystujący parametryczny bootstrap.
Karl Ove Hufthammer 18.04.15
6

Stopnie swobody dla modeli mieszanych są „problematyczne”. Aby przeczytać więcej na ten temat, możesz sprawdzić lmer, wartości p i cały ten post napisany przez Douglasa Batesa. Również R-sig-mieszane modele FAQ podsumowuje powody dlaczego jest to uciążliwe:

  • Zasadniczo nie jest jasne, czy rozkład zerowy obliczonego stosunku sum kwadratów jest tak naprawdę rozkładem F dla dowolnego wyboru mianownika stopni swobody. Chociaż dotyczy to specjalnych przypadków, które odpowiadają klasycznym projektom eksperymentalnym (zagnieżdżone, podzielone, losowy blok itp.), Najwyraźniej nie jest to prawdą w przypadku bardziej złożonych projektów (niezrównoważone, GLMM, korelacja czasowa lub przestrzenna itp.).
  • Dla każdego zaproponowanego prostego przepisu stopni swobody (ślad matrycy kapelusza itp.) Wydaje się, że istnieje co najmniej jeden dość prosty kontrprzykład, w którym przepis źle się zawodzi.
  • Inne sugerowane schematy aproksymacji df (Satterthwaite, Kenward-Roger, itp.) Byłyby najwyraźniej dość trudne do wdrożenia w lme4 / nlme,
    (...)
  • Ponieważ pierwotni autorzy lme4 nie są przekonani o użyteczności ogólnego podejścia do testowania w odniesieniu do przybliżonego rozkładu zerowego, a także z powodu narzutu związanego z kopaniem kodu w celu włączenia odpowiedniej funkcjonalności (jako łatki lub dodatku -on), jest mało prawdopodobne, aby sytuacja ta uległa zmianie w przyszłości.

FAQ podaje również kilka alternatyw

  • użyj MASS :: glmmPQL (używa starych reguł nlme w przybliżeniu równoważnych regułom „wewnętrzny-zewnętrzny” SAS) dla GLMM lub (n) lme dla LMM
  • Odgadnij mianownik df ze standardowych reguł (dla standardowych projektów) i zastosuj je do testów t lub F.
  • Uruchom model w lme (jeśli to możliwe) i użyj podanego tam mianownika df (który jest zgodny z prostą regułą „wewnętrzna-zewnętrzna”, która powinna odpowiadać odpowiedzi kanonicznej dla prostych / ortogonalnych projektów), zastosowanej do testów t lub F. Aby zapoznać się z wyraźną specyfikacją reguł, których używam, patrz strona 91 Pinheiro i Bates - ta strona jest dostępna w Google Books
  • używać SAS, Genstat (AS-REML), Stata?
  • Załóżmy nieskończony mianownik df (tj. Test Z / chi-kwadrat zamiast t / F), jeśli liczba grup jest duża (> 45? Wprowadzono różne ogólne zasady określające, jak duża jest „w przybliżeniu nieskończona”, w tym [w Angrist i Pischke's „Ekonometria w większości nieszkodliwa”], 42 (w hołdzie Douglasowi Adamsowi)

Ale jeśli interesują Cię przedziały ufności, istnieją lepsze podejścia, np. Oparte na bootstrap, jak sugeruje Karl Ove Hufthammer w swojej odpowiedzi lub te zaproponowane w FAQ.

Tim
źródło
„Odgadnij mianownik df ze standardowych zasad (dla standardowych projektów) i zastosuj je do testów t lub F”; Naprawdę chciałbym, żeby ktoś mógł to rozwinąć. Na przykład w przypadku projektu zagnieżdżonego (np. Pacjenci kontra kontrole, kilka próbek na pacjenta; przy czym identyfikator pacjenta jest efektem losowym), w jaki sposób uzyskujemy stopnie swobody dla takiego projektu?
Arnaud A