Jak uzyskać „ogólną” wartość p i wielkość efektu dla czynnika jakościowego w modelu mieszanym (Lme4)?

28

Chciałbym uzyskać wartość p i wielkość efektu niezależnej zmiennej kategorialnej (z kilkoma poziomami) - to jest „ogólnie”, a nie dla każdego poziomu osobno, tak jak normalne wyjście z lme4R. To jest jak rzecz, którą ludzie zgłaszają, uruchamiając ANOVA.

Jak mogę to zdobyć?

użytkownik3288202
źródło
Jakich statystyk dokładnie chcesz? Za pomocą tej anova()funkcji można uzyskać tabelę anova z liniowymi modelami mieszanymi, tak jak w przypadku modeli liniowych.
smillig
Próbowałem anova (), ale daje mi wartość Df, Sum Sq, Mean Sq i F. Nie widzę wielkości efektu i wartości p. Czy masz jakiś pomysł na ten temat?
user3288202
1
Według wielkości efektu, masz na myśli coś jak odpowiednik ? Jeśli chodzi o wartości p, trwa długa i merytoryczna debata wokół ich szacunków i ich wdrażania . Więcej informacji można znaleźć w dyskusji w tym pytaniu . R2lme4
smillig
Dzięki za link, Smilig. Czy to oznacza, że ​​ponieważ istnieje problem z obliczaniem wartości p, rozmiar efektu współczynnika w ogóle jest również problemem?
user3288202
Nie są to bezpośrednio powiązane problemy. Należy jednak pamiętać, że liniowy model mieszany nie zachowuje się dokładnie tak, jak model liniowy bez efektów losowych, więc miara, która może być odpowiednia dla modelu liniowego, niekoniecznie uogólnia się na modele mieszane.
smillig

Odpowiedzi:

48

Obie wspomniane koncepcje (wartości p i wielkości efektu liniowych modeli mieszanych) mają nieodłączne problemy. Jeśli chodzi o wielkość efektu , cytując Douga Batesa, oryginalnego autora lme4,

R2

Aby uzyskać więcej informacji, możesz przejrzeć ten wątek , ten wątek i tę wiadomość . Zasadniczo problem polega na tym, że nie ma uzgodnionej metody włączania i dekompozycji wariancji z efektów losowych w modelu. Jednak stosuje się kilka standardów. Jeśli spojrzysz na Wiki skonfigurowaną dla / przez listę mailingową r-sig-mixed-models , wymienionych jest kilka podejść.

Jedna z sugerowanych metod dotyczy korelacji między wartościami dopasowanymi a obserwowanymi. Można to zaimplementować w R, jak sugeruje Jarrett Byrnes w jednym z tych wątków:

r2.corr.mer <- function(m) {
  lmfit <-  lm(model.response(model.frame(m)) ~ fitted(m))
  summary(lmfit)$r.squared
}

Powiedzmy na przykład, że szacujemy następujący liniowy model mieszany:

set.seed(1)
d <- data.frame(y = rnorm(250), x = rnorm(250), z = rnorm(250),
                g = sample(letters[1:4], 250, replace=T)       )
library(lme4)
summary(fm1 <- lmer(y ~ x + (z | g), data=d))
# Linear mixed model fit by REML ['lmerMod']
# Formula: y ~ x + (z | g)
#    Data: d
# REML criterion at convergence: 744.4
# 
# Scaled residuals: 
#     Min      1Q  Median      3Q     Max 
# -2.7808 -0.6123 -0.0244  0.6330  3.5374 
# 
# Random effects:
#  Groups   Name        Variance Std.Dev. Corr 
#  g        (Intercept) 0.006218 0.07885       
#           z           0.001318 0.03631  -1.00
#  Residual             1.121439 1.05898       
# Number of obs: 250, groups: g, 4
# 
# Fixed effects:
#             Estimate Std. Error t value
# (Intercept)  0.02180    0.07795   0.280
# x            0.04446    0.06980   0.637
# 
# Correlation of Fixed Effects:
#   (Intr)
# x -0.005

Możemy obliczyć wielkość efektu za pomocą funkcji zdefiniowanej powyżej:

r2.corr.mer(fm1)
# [1] 0.0160841

Ω02

1-var(residuals(fm1))/(var(model.response(model.frame(fm1))))
# [1] 0.01173721  # Usually, it would be even closer to the value above

W odniesieniu do wartości p jest to kwestia o wiele bardziej kontrowersyjna (przynajmniej w społeczności R / lme4). Zobacz dyskusje w pytaniach tutaj , tutaj i tutaj, wśród wielu innych. Odwołując się ponownie do strony Wiki, istnieje kilka podejść do testowania hipotez dotyczących efektów w liniowych modelach mieszanych. Lista od „najgorszego do najlepszego” (według autorów strony Wiki, która, jak wierzę, obejmuje Douga Batesa, a także Bena Bolkera, który dużo tu wnosi):

  • Wald Z-testy
  • Dla zrównoważonych, zagnieżdżonych LMM, w których można obliczyć df: testy t-Wald
  • Test współczynnika wiarygodności, albo przez skonfigurowanie modelu, aby parametr mógł być izolowany / upuszczany (przez anovalub drop1), lub poprzez obliczanie profili prawdopodobieństwa
  • MCMC lub parametryczne przedziały ufności bootstrapu

Zalecają metodę próbkowania Monte Carlo w łańcuchu Markowa, a także wymieniają szereg możliwości realizacji tego z pseudo i całkowicie bayesowskich metod, wymienionych poniżej.

Pseudobayesian:

  • Próbkowanie post-hoc, zazwyczaj (1) przy założeniu płaskich priorytetów i (2) rozpoczynając od MLE, być może przy użyciu przybliżonego oszacowania wariancji-kowariancji w celu wybrania rozkładu kandydata
  • Via mcmcsamp(jeśli dostępne dla twojego problemu: tj. LMM z prostymi efektami losowymi - nie GLMM lub złożonymi efektami losowymi)
    Via pvals.fncw languageRpakiecie, opakowanie dla mcmcsamp)
  • W AD Model Builder, ewentualnie za pośrednictwem glmmADMBpakietu (użyj mcmc=TRUEopcji) lub R2admbpakietu (napisz własną definicję modelu w AD Model Builder), lub poza R
  • Poprzez simfunkcję z armpakietu (symuluje tylną tylko dla współczynników beta (efektu stałego)

Podejścia w pełni bayesowskie:

  • Poprzez MCMCglmmpaczkę
  • Za pomocą glmmBUGS(a WinBUGS owijki / R Interface)
  • Używanie JAGS / WinBUGS / OpenBUGS itp. Za pośrednictwem pakietów rjags/ r2jags/ R2WinBUGS/BRugs

Dla ilustracji, aby pokazać, jak to może wyglądać, poniżej znajduje się MCMCglmmszacunek przy użyciu MCMCglmmpakietu, który zobaczysz daje podobne wyniki jak w powyższym modelu i ma pewne bayesowskie wartości p:

library(MCMCglmm)
summary(fm2 <- MCMCglmm(y ~ x, random=~us(z):g, data=d))
# Iterations = 3001:12991
# Thinning interval  = 10
#  Sample size  = 1000 
# 
#  DIC: 697.7438 
# 
#  G-structure:  ~us(z):g
# 
#       post.mean  l-95% CI u-95% CI eff.samp
# z:z.g 0.0004363 1.586e-17 0.001268    397.6
# 
#  R-structure:  ~units
# 
#       post.mean l-95% CI u-95% CI eff.samp
# units    0.9466   0.7926    1.123     1000
# 
#  Location effects: y ~ x 
# 
#             post.mean l-95% CI u-95% CI eff.samp pMCMC
# (Intercept)  -0.04936 -0.17176  0.07502     1000 0.424
# x            -0.07955 -0.19648  0.05811     1000 0.214

Mam nadzieję, że to trochę pomaga. Myślę, że najlepszą radą dla kogoś, kto zaczyna od liniowych modeli mieszanych i próbuje oszacować je w R, jest przeczytanie często zadawanych pytań na Wiki, skąd pochodzi większość tych informacji. Jest to doskonałe źródło wszelkiego rodzaju motywów z efektami mieszanymi od podstawowego do zaawansowanego oraz od modelowania do kreślenia.

smillig
źródło
Dziękuję bardzo smilig. Może więc nie zgłosić wielkości efektu dla ogólnych parametrów.
user3288202
r2
3
+6, imponująco jasne, kompleksowe i dokładnie opatrzone adnotacjami.
gung - Przywróć Monikę
1
dodatkowo możesz rzucić okiem na pakiet afex, a zwłaszcza na funkcję mieszaną. patrz tutaj
początek
6

Jeśli chodzi o obliczanie wartości istotności ( p ), Luke (2016) Oceniając istotność w liniowych modelach efektów mieszanych w R podaje, że optymalną metodą jest przybliżenie Kenwarda-Rogera lub Satterthwaite dla stopni swobody (dostępne w R z pakietami takimi jak lmerTestlub afex).

Abstrakcyjny

Modele z efektami mieszanymi są coraz częściej wykorzystywane w analizie danych eksperymentalnych. Jednak w pakiecie lme4 w R standardy oceny istotności efektów stałych w tych modelach (tj. Uzyskiwanie wartości p) są nieco niejasne. Istnieją ku temu dobre powody, ale ponieważ w wielu przypadkach badacze korzystający z tych modeli są zobowiązani do zgłaszania wartości p, potrzebna jest metoda oceny istotności wyników modelu. W niniejszym artykule przedstawiono wyniki symulacji pokazujących, że dwie najczęstsze metody oceny istotności, wykorzystujące testy współczynnika wiarygodności i zastosowanie rozkładu Z do wartości Wald t z wyników modelu (t-as-z), są nieco antykonserwatywne, szczególnie dla mniejszych próbek. Inne metody oceny istotności,Wyniki tych symulacji sugerują, że wskaźniki błędów Typu 1 są najbliższe 0,05, gdy modele są montowane za pomocą REML, a wartości p są uzyskiwane za pomocą przybliżeń Kenwarda-Rogera lub Satterthwaite, ponieważ te przybliżenia dały akceptowalne poziomy błędu Typu 1 nawet dla mniejszych próbki.

(podkreślenie dodane)

Pablo Bernabeu
źródło
4
+1 Dziękujemy za udostępnienie tego linku. Powiem krótko, że przybliżenie Kenwarda-Rogera jest dostępne w lmerTestpakiecie.
ameba mówi Przywróć Monikę
5

Korzystam z lmerTestpaczki. To dogodnie obejmuje oszacowanie wartości p na anova()wyjściu dla moich analiz MLM, ale nie daje wielkości efektu z powodów podanych w innych postach tutaj.

Bruna
źródło
1
W moim przypadku wolę porównywanie par przy użyciu lsmeans, ponieważ daje mi wszystkie pary kontrastów, w tym wartości p. Jeśli użyję lmerTest, będę musiał uruchomić model sześć razy z różnymi liniami bazowymi, aby zobaczyć wszystkie pary kontrastów.
user3288202