Vergleich eines gemischten Modells (Subjekt als Zufallseffekt) mit einem einfachen linearen Modell (Subjekt als fester Effekt)

Ich beende eine Analyse eines großen Datensatzes. Ich möchte das im ersten Teil der Arbeit verwendete lineare Modell verwenden und es mithilfe eines linearen gemischten Modells (LME) neu anpassen. Die LME wäre sehr ähnlich, mit der Ausnahme, dass eine der im Modell verwendeten Variablen als zufälliger Effekt verwendet würde. Diese Daten stammen aus vielen Beobachtungen (> 1000) in einer kleinen Gruppe von Probanden (~ 10), und ich weiß, dass die Modellierung des Effekts eines Probanden besser als zufälliger Effekt erfolgt (dies ist eine Variable, die ich verschieben möchte). Der R-Code würde folgendermaßen aussehen:

my_modelB <- lm(formula = A ~ B + C + D)    
lme_model <- lme(fixed=A ~ B + C, random=~1|D, data=my_data, method='REML')

Alles läuft gut und die Ergebnisse sind sehr ähnlich. Es wäre schön, wenn ich so etwas wie RLRsim oder einen AIC / BIC verwenden könnte, um diese beiden Modelle zu vergleichen und zu entscheiden, welches am besten geeignet ist. Meine Kollegen möchten die LME nicht melden, da es keine leicht zugängliche Möglichkeit gibt, die "bessere" auszuwählen, obwohl ich denke, dass die LME das geeignetere Modell ist. Irgendwelche Vorschläge?

Dies soll die Antwort von @ ocram ergänzen, da es zu lang ist, um als Kommentar zu posten. Ich würde es A ~ B + Cals Ihr Nullmodell behandeln, damit Sie die statistische Signifikanz eines Dzufälligen Abschnitts auf einer Ebene in einem verschachtelten Modellaufbau beurteilen können . Wie ocram hervorhob, werden Regelmäßigkeitsbedingungen verletzt, wenn , und die Likelihood-Ratio-Teststatistik (LRT) wird nicht notwendigerweise asymptotisch verteilt sein . Die Lösung, die mir beigebracht wurde, bestand darin, das LRT (dessen Bootstrap-Verteilung wahrscheinlich nicht ) parametrisch zu booten und einen Bootstrap-p-Wert wie folgt zu berechnen:χ 2 χ $H_0: \sigma^2 = 0$$\chi^2$$\chi^2$

library(lme4)
my_modelB <- lm(formula = A ~ B + C)
lme_model <- lmer(y ~ B + C + (1|D), data=my_data, REML=F)
lrt.observed <- as.numeric(2*(logLik(lme_model) - logLik(my_modelB)))
nsim <- 999
lrt.sim <- numeric(nsim)
for (i in 1:nsim) {
    y <- unlist(simulate(mymodlB))
    nullmod <- lm(y ~ B + C)
    altmod <- lmer(y ~ B + C + (1|D), data=my_data, REML=F)
    lrt.sim[i] <- as.numeric(2*(logLik(altmod) - logLik(nullmod)))
}
mean(lrt.sim > lrt.observed) #pvalue

Der Anteil der Bootstrap-LRTs, der extremer ist als der beobachtete LRT, ist der p-Wert.

Ich bin mir nicht ganz sicher, welches Modell passt, wenn Sie die lme-Funktion verwenden. (Ich denke, der zufällige Effekt soll einer Normalverteilung mit dem Mittelwert Null folgen?). Das lineare Modell ist jedoch ein Sonderfall des gemischten Modells, wenn die Varianz des Zufallseffekts Null ist. Obwohl einige technische Schwierigkeiten bestehen (da in der Grenze des Parameterraums für die Varianz liegt), sollte es möglich sein, zu testen vs ...H 0 : v a r i a n c e = 0 H 1 : v a r i a n c e > $0$$H_0:variance = 0$$H_1: variance > 0$

BEARBEITEN

Um Verwirrung zu vermeiden: Der oben erwähnte Test wird manchmal verwendet, um zu entscheiden, ob der zufällige Effekt signifikant ist oder nicht ... aber nicht, um zu entscheiden, ob er in einen festen Effekt umgewandelt werden soll oder nicht.