Jak uzyskać wartość p (sprawdzić istotność) efektu w mieszanym modelu Lme4?

56

Używam lme4 w R, aby dopasować model mieszany

lmer(value~status+(1|experiment)))

gdzie wartość jest ciągła, status i eksperyment są czynnikami, a ja rozumiem

Linear mixed model fit by REML 
Formula: value ~ status + (1 | experiment) 
  AIC   BIC logLik deviance REMLdev
 29.1 46.98 -9.548    5.911    19.1
Random effects:
 Groups     Name        Variance Std.Dev.
 experiment (Intercept) 0.065526 0.25598 
 Residual               0.053029 0.23028 
Number of obs: 264, groups: experiment, 10

Fixed effects:
            Estimate Std. Error t value
(Intercept)  2.78004    0.08448   32.91
statusD      0.20493    0.03389    6.05
statusR      0.88690    0.03583   24.76

Correlation of Fixed Effects:
        (Intr) statsD
statusD -0.204       
statusR -0.193  0.476

Skąd mam wiedzieć, że wpływ statusu jest znaczący? R zgłasza tylko wartości a nie wartości .ptp

ECII
źródło
1
Przechodząc do odpowiedzi udzielonych na to pytanie, można się zastanawiać, co tak naprawdę interesuje OP: testowanie współczynników na wartości zerowej (test waniliowy przeprowadzany w regularnej regresji liniowej na wartości zerowej ) lub testowanie w celu zminimalizowania wariancji (test otrzymujemy z wielu rodzajów ANOVA). Ci dwaj dążą do różnych rzeczy. Oświecająca odpowiedź, choć nie dotyczy modeli z efektami mieszanymi, znajduje się tutaj . H 0 : β = β null FtH0:β=βnullF
Firebug

Odpowiedzi:

61

Wiele informacji na ten temat można znaleźć w FAQ GLMM . Jednak w twoim konkretnym przypadku sugerowałbym użycie

library(nlme)
m1 <- lme(value~status,random=~1|experiment,data=mydata)
anova(m1)

ponieważ nie potrzebujesz żadnych rzeczy, które lmeroferuje (większa prędkość, obsługa skrzyżowanych efektów losowych, GLMM ...). lmepowinien dać ci dokładnie takie same oszacowania współczynników i wariancji, ale także obliczy dla ciebie wartości df i p (które mają sens w „klasycznym” projekcie, takim, jaki wydajesz się mieć). Możesz także rozważyć losowy termin ~status|experiment(pozwalający na zmianę efektów statusu między blokami lub równoważnie włączając interakcję między stanami). Powyższe plakaty mają również rację, że twoje tstatystyki są tak duże, że twoja wartość p na pewno wyniesie <0,05, ale mogę sobie wyobrazić, że chciałbyś „prawdziwych” wartości p.

Ben Bolker
źródło
3
Nie wiem o tej odpowiedzi. lmerrównie łatwo może zgłaszać te same rodzaje wartości p, ale nie robi tego z ważnych powodów. Myślę, że to komentarz, że są tutaj jakieś „prawdziwe” wartości p, które mnie wkurzą. Można argumentować, że można znaleźć jedną możliwą granicę odcięcia i że wszelkie uzasadnione wartości graniczne zostały przekroczone. Ale nie można argumentować, że istnieje prawdziwa wartość p.
Jan
11
W przypadku klasycznego projektu (zbalansowanego, zagnieżdżonego itp.) Myślę, że rzeczywiście mogę argumentować, że istnieje prawdziwa wartość p, tj. Prawdopodobieństwo uzyskania oszacowania beta o obserwowanej wielkości lub większej, jeśli hipoteza zerowa (beta = 0) były fałszywe ... lme4 nie zapewnia tych mianowników df, jak sądzę, ponieważ trudniej jest ogólnie wykryć strukturę modelu lme4, gdy określony model jest taki, w którym działałaby jakaś heurystyka dla obliczenia klasycznego mianownika df ...
Ben Bolker,
spróbuj summary(m1)zamiast tego (używam tego z pakietem nlme)
jena
36

Możesz użyć pakietu lmerTest . Po prostu zainstaluj / załaduj, a modele lmerów zostaną przedłużone. Więc np

library(lmerTest)
lmm <- lmer(value~status+(1|experiment)))
summary(lmm)
anova(lmm)

dałby wyniki z wartościami p. Jeśli wartości p są właściwym wskazaniem, jest nieco kwestionowane, ale jeśli chcesz je mieć, to jest sposób na ich uzyskanie.

pbx101
źródło
28

Jeśli potrafisz poradzić sobie z porzucaniem wartości p ( i powinieneś ), możesz obliczyć współczynnik wiarygodności, który reprezentowałby wagę dowodów na wpływ statusu poprzez:

#compute a model where the effect of status is estimated
unrestricted_fit = lmer(
    formula = value ~ (1|experiment) + status
    , REML = F #because we want to compare models on likelihood
)
#next, compute a model where the effect of status is not estimated
restricted_fit = lmer(
    formula = value ~ (1|experiment)
    , REML = F #because we want to compare models on likelihood
)
#compute the AIC-corrected log-base-2 likelihood ratio (a.k.a. "bits" of evidence)
(AIC(restricted_fit)-AIC(unrestricted_fit))*log2(exp(1))
Mike Lawrence
źródło
16
Zauważ, że współczynniki prawdopodobieństwa są asymptotyczne, tzn. Nie uwzględniają niepewności w oszacowaniu wariancji rezydualnej ...
Ben Bolker
5
Interesuje mnie twoja ostatnia linia. Jaka jest interpretacja wyniku? Czy są dostępne źródła?
mguzmann
13

Problem polega na tym, że obliczanie wartości p dla tych modeli nie jest trywialne, patrz omówienie tutaj, więc autorzy lme4pakietu celowo postanowili nie uwzględniać wartości p na wyjściu. Możesz znaleźć metodę ich obliczania, ale niekoniecznie będą one prawidłowe.

Michelle
źródło
9

Zastanów się, o co pytasz. Jeśli chcesz tylko wiedzieć, czy ogólna wartość p dla efektu statusu przekracza jakąś arbitralną wartość odcięcia, na przykład 0,05, to jest to łatwe. Najpierw chcesz poznać ogólny efekt. Możesz to uzyskać od anova.

m <- lmer(...) #just run your lmer command but save the model
anova(m)

Teraz masz wartość F. Możesz to wziąć i sprawdzić w niektórych tabelach F. Wybierz najniższy możliwy nominał. stopnie swobody. Punkt odcięcia będzie około 20. Twoje F może być większe, ale mogę się mylić. Nawet jeśli nie jest, spójrz na liczbę stopni swobody w porównaniu z konwencjonalnymi obliczeniami ANOVA tutaj, korzystając z liczby przeprowadzonych eksperymentów. Utrzymanie tej wartości w granicach około 5 dla granicy. Teraz łatwo przekazujesz go w swoim badaniu. „Prawdziwy” df dla twojego modelu będzie czymś wyższym, ponieważ modelujesz każdy punkt danych w przeciwieństwie do agregowanych wartości, które modelowałaby ANOVA.

Jeśli naprawdę chcesz dokładnej wartości p, nie ma czegoś takiego, chyba że zechcesz wypowiedzieć się na jej temat. Jeśli czytasz Pinheiro i Bates (2001, i być może kilka innych książek na ten temat ... zobacz inne linki w tych odpowiedziach) i wymyślisz argument za konkretnym df, możesz go użyć. Ale tak naprawdę nie szukasz dokładnej wartości p. Wspominam o tym, ponieważ dlatego nie powinieneś zgłaszać dokładnej wartości p, tylko że twój punkt odcięcia został przekroczony.

Naprawdę powinieneś rozważyć odpowiedź Mike'a Lawrence'a, ponieważ cały pomysł trzymania się punktu przejścia dla wartości p jako ostatecznej i najważniejszej informacji do wyciągnięcia z twoich danych jest ogólnie mylny (ale może nie być w twoim przypadku, ponieważ nie „ naprawdę mam wystarczającą ilość informacji, aby wiedzieć). Mike używa ciekawej wersji kalkulacji LR dla zwierząt domowych, ale może być trudno znaleźć na jej temat wiele dokumentacji. Jeśli spojrzysz na wybór i interpretację modelu za pomocą AIC, może ci się spodobać.

Jan
źródło
9

Edycja: Ta metoda nie jest już obsługiwana w nowszych wersjach Lme4. Skorzystaj z pakietu lmerTest, zgodnie z sugestią pbx101 w tej odpowiedzi .

Na liście R znajduje się post autora lme4, dlaczego wartości p nie są wyświetlane. Sugeruje zamiast tego użycie próbek MCMC, które robisz za pomocą pvals.fnc z pakietu languageR:

library("lme4")
library("languageR")
model=lmer(...)
pvals.fnc(model)

Przykład i szczegóły można znaleźć na stronie http://www2.hawaii.edu/~kdrager/MixedEffectsModels.pdf .

Jeff
źródło
3
lme4 już tego nie obsługuje. Ten post można zaktualizować, aby oszczędzić ludziom konieczności dowiedzenia się tego tak, jak właśnie to zrobiłem.
timothy.s.lau,
5

Czy chcesz wiedzieć, czy połączony efekt statusma znaczący wpływ na value? Jeśli tak, możesz użyć Anovafunkcji w carpakiecie (nie mylić z anovafunkcją w bazie R).

dat <- data.frame(
  experiment = sample(c("A","B","C","D"), 264, replace=TRUE), 
  status = sample(c("D","R","A"), 264, replace=TRUE), 
  value = runif(264)   
)
require(lme4)
(fm <- lmer(value~status+(1|experiment), data=dat))

require(car)
Anova(fm)

Zobacz ?Anovapo załadowaniu carpaczki.

smillig
źródło
Jakiś pomysł, jak car::Anova()uniknąć lepkich problemów związanych z obliczaniem wartości p, które łączy Michelle?
Mike Lawrence
Nie wiem, ale domyślam się, że pozwala uniknąć lepkich problemów, ignorując je! Po ponownym przeczytaniu oryginalnego postu czuję, że mogłem źle zrozumieć pytanie. Jeśli OP chce dokładnych wartości p dla stałych parametrów efektów, ma kłopoty. Ale jeśli PO chce tylko wiedzieć, czy są znaczące, myślę, że wartości t są większe niż jakakolwiek niepewność w sposobie obliczenia dokładnej wartości p. (Innymi słowy, są znaczące.)
smillig
1
Myślę, że zdecydowanie dobrym pomysłem było przekierowanie do obliczeń ANOVA, aby dowiedzieć się, jaki jest ogólny efekt statystyk, ale nie jestem pewien, czy ustalanie wartości p jest dobre. Zwykłe anovapolecenie da ci F.
Jan
Myślę, że to trochę trudniejsze niż pozorne. Uruchamianie ANOVA jest poprawne, gdy chcesz zminimalizować wariancję, ale myślę, że z treści pytania wynika, że ​​OP chce ustalić krańcowy efekt zmiennych, tj. Współczynniki testowe względem wartości zerowej.
Firebug,
0

Funkcja pvals.fncnie jest już obsługiwana przez lme4. Korzystając z pakietu pakiet lmerTest, można użyć innej metody do obliczenia wartości p, takiej jak przybliżenia Kenwarda-Rogera

model=lmer(value~status+1|experiment)
anova(model, ddf="Kenward-Roger")
Stefano Cacciatore
źródło
0

Po prostu załadowanie pakietu afex wydrukuje wartości p na wyjściu funkcji lmer z pakietu lme4 (nie musisz używać afexa, po prostu go załaduj):

library(lme4)  #for mixed model
library(afex)  #for p-values

Spowoduje to automatyczne dodanie kolumny wartości p do danych wyjściowych lmera (twój model) dla ustalonych efektów.

Maryam Nasseri
źródło