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

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 .$t$$p$

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.

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.

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

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.

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

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 .

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.

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

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.