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

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

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,

$R^2$

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

$\Omega^{2}_{0}$

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.

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)

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.