Wie erhält man einen "gesamten" p-Wert und eine Effektgröße für einen kategorialen Faktor in einem gemischten Modell (lme4)?

28

Ich möchte einen p-Wert und eine Effektgröße einer unabhängigen kategorialen Variablen (mit mehreren Ebenen) erhalten - das ist "insgesamt" und nicht für jede Ebene separat, wie es die normale Ausgabe von lme4in R ist. Es ist genau wie das, was die Leute berichten, wenn sie eine ANOVA betreiben.

Wie kann ich das bekommen?

user3288202
quelle
Welche Statistiken möchten Sie genau? Sie können die anova()Funktion verwenden, um eine Anova-Tabelle mit linearen gemischten Modellen zu erhalten, genau wie bei linearen Modellen.
Smillig
Ich habe anova () ausprobiert, aber es gibt mir den Wert Df, Sum Sq, Mean Sq und F. Ich sehe keine Effektgröße und keinen p-Wert. Hast du eine Idee dazu?
User3288202
1
Meinen Sie mit Effektgröße so etwas wie ein Äquivalent zu ? In Bezug auf p-Werte gibt es eine lange und umfassende Debatte über ihre Schätzung und ihre Umsetzung in . Schauen Sie sich die Diskussion in dieser Frage für weitere Details an. R2lme4
Smillig
Danke für den Link, Smilig. Bedeutet das, dass aufgrund eines Problems bei der Berechnung des p-Werts auch die Effektgröße des Faktors insgesamt ein Problem ist?
User3288202
Sie sind nicht direkt verwandte Themen. Beachten Sie jedoch, dass sich ein lineares gemischtes Modell nicht genau wie ein lineares Modell ohne zufällige Effekte verhält, sodass ein für das lineare Modell geeignetes Maß nicht unbedingt auf gemischte Modelle verallgemeinert werden muss.
Smillig

Antworten:

48

Beide von Ihnen genannten Konzepte (p-Werte und Effektgrößen von linearen gemischten Modellen) weisen inhärente Probleme auf. In Bezug auf die Effektgröße zitiert Doug Bates, der ursprüngliche Autor von lme4,

Angenommen, man möchte ein R 2 definierenR2 -Maß Erachtens ein Argument für die Behandlung der bestraften Restsumme von Quadraten aus einem linearen gemischten Modell auf die gleiche Weise wie für die Restsumme von Quadraten aus einem linearen Modell angeführt werden. Oder man könnte nur die Restsumme der Quadrate ohne die Strafe oder die minimale Restsumme der Quadrate verwenden, die aus einer gegebenen Menge von Begriffen erhalten werden kann, die einer unendlichen Präzisionsmatrix entspricht. Ich weiß es nicht wirklich. Es hängt davon ab, was Sie zu charakterisieren versuchen.

Weitere Informationen finden Sie in diesem Thread , diesem Thread und dieser Nachricht . Grundsätzlich besteht das Problem darin, dass es keine vereinbarte Methode zur Einbeziehung und Zerlegung der Varianz aus den Zufallseffekten in das Modell gibt. Es werden jedoch einige Standards verwendet. Wenn Sie sich das Wiki ansehen, das für / von der Mailingliste r-sig-mixed-models eingerichtet wurde , werden einige Ansätze aufgelistet.

Eine der vorgeschlagenen Methoden untersucht die Korrelation zwischen den angepassten und den beobachteten Werten. Dies kann in R implementiert werden, wie von Jarrett Byrnes in einem dieser Threads vorgeschlagen:

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

Nehmen wir zum Beispiel an, wir schätzen das folgende lineare Mischmodell:

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

Wir können die Effektgröße mit der oben definierten Funktion berechnen:

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

In Bezug auf die p-Werte ist dies ein viel strittigeres Thema (zumindest in der R / lme4Community). Sehen Sie sich die Diskussionen in den Fragen hier , hier und hier unter vielen anderen an. Wenn Sie erneut auf die Wiki-Seite verweisen, gibt es einige Ansätze, um Hypothesen zu Effekten in linearen gemischten Modellen zu testen. Aufgelistet von "Worst to Best" (laut den Autoren der Wiki-Seite, von denen ich glaube, dass sie Doug Bates sowie Ben Bolker umfasst, der hier viel beiträgt):

  • Wald Z-Tests
  • Für ausgeglichene, verschachtelte LMMs, für die df berechnet werden kann: Wald-t-Tests
  • Likelihood - Ratio - Test, indem das Modell so eingerichtet wird, dass der Parameter isoliert / gelöscht werden kann (über anovaoder)drop1 ), oder Likelihood-Profile berechnet werden
  • MCMC- oder parametrische Bootstrap-Konfidenzintervalle

Sie empfehlen den Markov-Ketten-Monte-Carlo-Stichprobenansatz und listen auch eine Reihe von Möglichkeiten auf, dies ausgehend von den nachstehend aufgeführten pseudo- und vollständig bayesianischen Ansätzen umzusetzen.

Pseudo-Bayesian:

  • Post-hoc-Stichprobe, wobei normalerweise (1) flache Prioritäten angenommen werden und (2) von der MLE ausgegangen wird, wobei möglicherweise die ungefähre Varianz-Kovarianz-Schätzung verwendet wird, um eine Kandidatenverteilung zu wählen
  • Via mcmcsamp(falls für Ihr Problem verfügbar: zB LMMs mit einfachen Zufallseffekten - keine GLMMs oder komplexen Zufallseffekten)
    Via pvals.fncim languageRPaket, ein Wrapper fürmcmcsamp )
  • In AD Model Builder möglicherweise über das glmmADMBPaket (verwenden Sie die mcmc=TRUEOption) oder das R2admbPaket (schreiben Sie Ihre eigene Modelldefinition in AD Model Builder) oder außerhalb von R
  • Über die simFunktion aus dem armPaket (simuliert das Posterior nur für die Beta (Fixed-Effect) -Koeffizienten)

Voll Bayesianische Ansätze:

  • Über das MCMCglmmPaket
  • Verwenden glmmBUGS(eine WinBUGS-Wrapper / R- Schnittstelle)
  • Verwenden von JAGS / WinBUGS / OpenBUGS usw. über die Pakete rjags/ r2jags/ R2WinBUGS/BRugs

Zur Veranschaulichung, um zu zeigen, wie dies aussehen könnte, wird im Folgenden eine MCMCglmmSchätzung unter Verwendung des MCMCglmmPakets angegeben, das ähnliche Ergebnisse wie das obige Modell liefert und eine Art von Bayes'schen p-Werten aufweist:

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

Ich hoffe das hilft etwas. Ich denke, der beste Rat für jemanden, der mit linearen gemischten Modellen beginnt und versucht, sie in R zu schätzen, ist, die Wiki-FAQs zu lesen, aus denen die meisten dieser Informationen stammen. Es ist eine hervorragende Ressource für alle Arten von Themen mit gemischten Effekten, von einfach bis fortgeschritten und von der Modellierung bis zum Plotten.

smillig
quelle
Vielen Dank smilig. Daher kann ich möglicherweise keine Effektgröße für die Gesamtparameter angeben.
User3288202
r2
3
+6, beeindruckend klar, umfassend und gründlich kommentiert.
gung - Wiedereinsetzung von Monica
1
ausserdem kann man sich das afex-paket und vor allem die mixed-funktion ansehen. siehe hier
beginneR
6

Zur Berechnung der Signifikanz ( p ) -Werte gibt Luke (2016) an , dass die optimale Methode für die Berechnung der Signifikanz in linearen Mischeffektmodellen in R entweder die Kenward-Roger- oder die Satterthwaite-Näherung für Freiheitsgrade ist (verfügbar in R mit Paketen wie lmerTestoder afex).

Abstrakt

Modelle mit gemischten Effekten werden bei der Analyse experimenteller Daten immer häufiger verwendet. Im lme4-Paket in R sind die Standards für die Bewertung der Signifikanz fester Effekte in diesen Modellen (dh das Erhalten von p-Werten) jedoch etwas vage. Es gibt gute Gründe dafür, aber da Forscher, die diese Modelle verwenden, in vielen Fällen verpflichtet sind, p-Werte zu melden, ist eine Methode zur Bewertung der Signifikanz der Modellausgabe erforderlich. Dieser Aufsatz berichtet über die Ergebnisse von Simulationen, aus denen hervorgeht, dass die beiden gängigsten Methoden zur Bewertung der Signifikanz unter Verwendung von Likelihood-Ratio-Tests und Anwendung der z-Verteilung auf die Wald-t-Werte aus der Modellausgabe (t-as-z) etwas anti-konservativ sind. speziell für kleinere Probengrößen. Andere Methoden zur Bewertung der Signifikanz,Die Ergebnisse dieser Simulationen legen nahe, dass die Fehlerraten von Typ 1 am nächsten bei 0,05 liegen, wenn Modelle unter Verwendung von REML angepasst werden und p-Werte unter Verwendung der Kenward-Roger- oder Satterthwaite-Näherungen abgeleitet werden, da beide Näherungen akzeptable Fehlerraten von Typ 1 auch für kleinere ergaben Proben.

(Betonung hinzugefügt)

Pablo Bernabeu
quelle
4
+1 Danke, dass du diesen Link geteilt hast. Ich möchte nur kurz darauf hinweisen, dass die Kenward-Roger-Näherung im lmerTestPaket enthalten ist.
Amöbe sagt Reinstate Monica
5

Ich benutze das lmerTestPaket. Dies beinhaltet zweckmäßigerweise eine Schätzung des p-Werts in der anova()Ausgabe für meine MLM-Analysen, gibt jedoch aus den in anderen Posts hier angegebenen Gründen keine Effektgröße an.

Bruna
quelle
1
In meinem Fall bevorzuge ich den paarweisen Vergleich mit lsmeans, da ich alle Kontrastpaare einschließlich der p-Werte erhalte. Wenn ich lmerTest verwende, muss ich das Modell sechs Mal mit verschiedenen Basislinien ausführen, um alle Kontrastpaare zu sehen.
User3288202