Ich habe einen unsymmetrischen Datensatz mit wiederholten Messungen zur Analyse und ich habe gelesen, dass die Art und Weise, wie die meisten statistischen Pakete mit ANOVA umgehen (dh die Summe der Quadrate des Typs III), falsch ist. Daher würde ich gerne ein gemischtes Effektmodell verwenden, um diese Daten zu analysieren. Ich habe viel über gemischte Modelle in gelesen R, aber ich bin immer noch sehr neu in RModellen mit gemischten Effekten und nicht sehr zuversichtlich, dass ich die Dinge richtig mache. Beachten Sie, dass ich mich noch nicht vollständig von "traditionellen" Methoden scheiden lassen kann und weiterhin Werte und Post-Hoc-Tests benötige .$p$

Ich würde gerne wissen, ob der folgende Ansatz Sinn macht oder ob ich etwas schrecklich Falsches tue. Hier ist mein Code:

# load packages
library(lme4)
library(languageR)
library(LMERConvenienceFunctions)
library(coda)
library(pbkrtest)

# import data
my.data <- read.csv("data.csv")

# create separate data frames for each DV & remove NAs
region.data <- na.omit(data.frame(time=my.data$time, subject=my.data$subject, dv=my.data$dv1))

# output summary of data
data.summary <- summary(region.data)

# fit model
# "time" is a factor with three levels ("t1", "t2", "t3")
region.lmer <- lmer(dv ~ time + (1|subject), data=region.data)

# check model assumptions
mcp.fnc(region.lmer)

# remove outliers (over 2.5 standard deviations)
rm.outliers <- romr.fnc(region.lmer, region.data, trim=2.5)
region.data <- rm.outliers$data
region.lmer <- update(region.lmer)

# re-check model assumptions
mcp.fnc(region.lmer)

# compare model to null model
region.lmer.null <- lmer(dv ~ 1 + (1|subject), data=region.data)
region.krtest <- KRmodcomp(region.lmer, region.lmer.null)

# output lmer summary
region.lmer.summary <- summary(region.lmer)

# run post hoc tests
t1.pvals <- pvals.fnc(region.lmer, ndigits=10, withMCMC=TRUE)

region.lmer <- lmer(dv ~ relevel(time,ref="t2") + (1|subject), data=region.data)
t2.pvals <- pvals.fnc(region.lmer, ndigits=10, withMCMC=TRUE)

region.lmer <- lmer(dv ~ relevel(time,ref="t3") + (1|subject), data=region.data)
t3.pvals <- pvals.fnc(region.lmer, ndigits=10, withMCMC=TRUE)

# Get mcmc mean and 50/95% HPD confidence intervals for graphs
# repeated three times and stored in a matrix (not shown here for brevity)
as.numeric(t1.pvals$fixed$MCMCmean)
as.numeric(t1.pvals$fixed$HPD95lower)
as.numeric(t1.pvals$fixed$HPD95upper)
HPDinterval(as.mcmc(t1.pvals$mcmc),prob=0.5)
    HPDinterval(as.mcmc(t1.pvals$mcmc),prob=0.5)

Einige spezifische Fragen, die ich habe:

1. Ist dies eine gültige Methode zur Analyse von Mixed-Effects-Modellen? Wenn nicht, was soll ich stattdessen tun?
2. Sind die von mcp.fnc ausgegebenen Kritik-Plots gut genug, um Modellannahmen zu überprüfen, oder sollte ich zusätzliche Schritte unternehmen?
3. Ich verstehe, dass gemischte Modelle nur gültig sind, wenn die Daten die Annahmen von Normalität und Homoskedastizität berücksichtigen. Wie kann ich beurteilen, was "ungefähr normal" ist und was nicht, indem ich die von mcp.fnc generierten Kritik-Plots betrachte? Muss ich nur ein Gefühl dafür bekommen, oder ist ihre Vorgehensweise vorgeschrieben? Wie robust sind gemischte Modelle in Bezug auf diese Annahmen?
4. Ich muss Unterschiede zwischen den drei Zeitpunkten für ~ 20 Merkmale (Biomarker) der Probanden in meiner Stichprobe bewerten. Ist das Anpassen und Testen von separaten Modellen für jedes akzeptabel, solange ich alle durchgeführten Tests berichte (signifikant oder nicht), oder benötige ich irgendeine Form der Korrektur für mehrere Vergleiche.

Um das Experiment etwas genauer zu gestalten, hier einige Details. Wir verfolgten eine Reihe von Teilnehmern in Längsrichtung, während sie sich einer Behandlung unterzogen. Wir haben vor Beginn der Behandlung und zu zwei Zeitpunkten danach eine Reihe von Biomarkern gemessen. Was ich sehen möchte, ist, ob zwischen den drei Zeitpunkten Unterschiede in diesen Biomarkern bestehen.

Das meiste, was ich hier mache, stütze ich auf dieses Tutorial , aber ich habe einige Änderungen vorgenommen, die auf meinen Bedürfnissen und den Dingen basieren, die ich lese. Die Änderungen, die ich vorgenommen habe, sind:

1. Relevel des "Zeit" -Faktors, um Vergleiche von t1-t2, t2-t3 und t1-t3 mit pvals.fnc (aus dem languageR-Paket) zu erhalten
2. Vergleiche mein gemischtes Modell mit dem Null-Modell unter Verwendung eines ungefähren F-Tests, der auf einem Kenward-Roger-Ansatz basiert (unter Verwendung des Pakets pbkrtest), anstatt eines Likelihood-Ratio-Tests (weil ich gelesen habe, dass Kenward-Roger momentan besser betrachtet wird).
3. Verwenden Sie das LMERConvenienceFunctions-Paket, um Annahmen zu überprüfen und Ausreißer zu entfernen (weil ich gelesen habe, dass gemischte Modelle sehr empfindlich gegenüber Ausreißern sind).

Ihre Frage (n) ist ein bisschen "groß", daher beginne ich mit einigen allgemeinen Kommentaren und Tipps.

Einige Hintergrundinformationen und nützliche Pakete

Sie sollten sich wahrscheinlich einige der Tutorial-Einführungen zur Verwendung gemischter Modelle ansehen. Ich würde empfehlen, mit Baayen et al. (2008) zu beginnen - Baayen ist der Autor von languageR. Barr et al. (2013) diskutieren einige Probleme mit der Zufallseffektstruktur, und Ben Bolker ist einer der aktuellen Entwickler von lme4. Baayen et al. Sind besonders gut für Ihre Fragen geeignet, da sie viel Zeit damit verbringen, über die Verwendung von Bootstrapping- / Permutationstests (das Zeug dahinter mcp.fnc) zu diskutieren .

Literatur

• Baayen, RH, D. Davidson und D. Bates (2008). Mixed-Effects-Modellierung mit gekreuzten Zufallseffekten für Objekte und Objekte . Journal of Memory and Language, 59 (390–412).
• Barr, Dale J., R. Levy, C. Scheepers und HJ Tily (2013). Random-Effects-Struktur für das Testen von Bestätigungshypothesen: Halten Sie sie maximal . Journal of Memory and Language, 68: 255– 278.
• Bolker, Benjamin M., ME Brooks, CJ Clark, SW Geange, JR Poulsen, MHH Stevens und J.-SS White (2009). Verallgemeinerte lineare Mischmodelle: ein praktischer Leitfaden für Ökologie und Evolution . Trends Ecol Evol, 24 (3): 127–35.

Florian Jaeger hat auch ein paar Sachen auf die praktische Seite des Mischmodelle auf seinem Labor des Blog .

Es gibt auch eine Reihe nützlicher R-Pakete zum Visualisieren und Testen gemischter Modelle wie lmerTestund effects. Das effectsPaket ist besonders schön, weil Sie damit die linearen Regressions- und Konfidenzintervalle zeichnen können, die sich hinter den Kulissen abspielen.

Passgenauigkeit und Modellvergleich

$p$lmerTest

anova()mer$\chi^2$$\chi^2$$p$-Werte zum direkten Vergleich zweier Modelle. Der Nachteil dabei ist, dass nicht sofort klar ist, wie gut Ihre Passform ist.

$t$summary()$|t|>2$fixef()

Sie sollten auch sicherstellen, dass keiner Ihrer festen Effekte zu stark korreliert ist - Multikollinearität verletzt die Modellannahmen. Florian Jaeger hat dazu einiges geschrieben , sowie einige mögliche Lösungen.

Mehrere Vergleiche

Ich werde Ihre Frage 4 etwas direkter beantworten. Die Antwort ist im Grunde die gleiche wie bei der herkömmlichen ANOVA. Leider scheint dies ein Punkt zu sein, an dem für die meisten Forscher eine große Unsicherheit besteht. (Es ist dasselbe wie bei der herkömmlichen ANOVA, da sowohl lineare Modelle mit gemischten Effekten als auch die ANOVA auf dem allgemeinen linearen Modell basieren. Gemischte Modelle haben lediglich einen zusätzlichen Term für die zufälligen Effekte.) Wenn Sie davon ausgehen, dass die Zeitfenster a ergeben Wenn Sie die Auswirkungen der Zeit vergleichen möchten, sollten Sie die Zeit in Ihr Modell einbeziehen. Dies bietet Ihnen übrigens auch eine bequeme Möglichkeit zu beurteilen, ob die Zeit einen Unterschied gemacht hat: Gibt es einen (festen) Haupteffekt für die Zeit? Der Nachteil dieser Route ist, dass Ihr Modell VIEL komplexer wird und die Single "super" Modell mit Zeit als Parameter wird wahrscheinlich länger dauern als drei kleinere Modelle ohne Zeit als Parameter. In der Tat ist das klassische Tutorial-Beispiel für gemischte Modellesleepstudy die Zeit als Parameter verwendet.

$t$foreachlme4$\chi^2$

Wenn Ihre Merkmale die abhängige Variable sind, müssen Sie ohnehin verschiedene Modelle berechnen, und dann können Sie die Ergebnisse mit dem AIC und dem BIC vergleichen.