Welche praktischen Auswirkungen hat das Schätzen in einem mehrstufigen Modell im Vergleich zum Nichtschätzen von Korrelationsparametern für zufällige Effekte?

Welche praktischen und interpretationsbezogenen Auswirkungen hat das Schätzen in einem mehrstufigen Modell im Vergleich zum Nichtschätzen von Korrelationsparametern für zufällige Effekte? Der praktische Grund, dies zu erfragen, ist, dass es im früheren Framework in R keine implementierte Methode zum Schätzen von p-Werten über MCMC-Techniken gibt, wenn Schätzungen im Modell der Korrelationen zwischen Parametern vorgenommen werden.

Wenn Sie sich zum Beispiel dieses Beispiel (Abschnitte unten) ansehen, welche praktischen Auswirkungen haben M2 und M3? Offensichtlich wird in einem Fall P5 nicht geschätzt und im anderen Fall.

Fragen

1. Aus praktischen Gründen (dem Wunsch, durch MCMC-Techniken einen p-Wert zu erhalten) kann es sinnvoll sein, ein Modell ohne Korrelationen zwischen zufälligen Effekten anzupassen, auch wenn P5 im Wesentlichen nicht Null ist. Wenn man dies tut und dann die p-Werte mit der MCMC-Technik schätzt, sind die Ergebnisse interpretierbar? (Ich weiß @ Ben Bolker zuvor erwähnt , dass „ die Kombination von Signifikanztests mit MCMC ist ein wenig zusammenhanglos, statistisch gesehen , obwohl ich den Drang zu verstehen , dies zu tun (immer Konfidenzintervall ist mehr erträglich)“ , also , wenn es Sie besser schlafen nachts so tun, als hätte ich vertrauensintervalle angegeben.)
2. Wenn man P5 nicht abschätzen kann, ist das dasselbe, als wenn man behauptet, es sei 0?
3. Wenn P5 wirklich nicht Null ist, wie werden dann die geschätzten Werte von P1-P4 beeinflusst?
4. Wenn P5 wirklich ungleich Null ist, wie werden dann die Fehlerschätzungen für P1-P4 beeinflusst?
5. Wenn P5 wirklich nicht Null ist, inwiefern sind Interpretationen eines Modells, die P5 nicht enthalten, fehlerhaft?

Ausgeborgt von der Antwort von @Mike Lawrence (diejenigen, die mehr wissen als ich, können dies durch die vollständige Modellnotation ersetzen, ich bin nicht ganz sicher, dass ich dies mit angemessener Genauigkeit tun kann):

M2: V1 ~ (1|V2) + V3 + (0+V3|V2)(Schätzungen P1 - P4)

M3: V1 ~ (1+V3|V2) + V3(Schätzungen P1-P5)

Parameter, die geschätzt werden könnten:

P1 : Ein globaler Schnittpunkt

P2 : Zufällige Effektabschnitte für V2 (dh für jede Ebene von V2 die Abweichung des Abschnitts dieser Ebene vom globalen Abschnitt)

P3 : Eine einzelne globale Schätzung für den Effekt (Steigung) von V3

P4 : Der Effekt von V3 in jeder Ebene von V2 (genauer gesagt, der Grad, in dem der V3-Effekt in einer bestimmten Ebene vom globalen Effekt von V3 abweicht), während eine Null-Korrelation zwischen den Intercept-Abweichungen und den V3-Effekt-Abweichungen über Ebenen hinweg erzwungen wird von V2.

P5 : Die Korrelation zwischen Schnittabweichungen und V3-Abweichungen über Ebenen von V2

Antworten, die aus einer ausreichend großen und umfassenden Simulation zusammen mit dem zugehörigen Code in R unter Verwendung von lmer abgeleitet wurden, wären akzeptabel.

Berücksichtigen Sie die Schlafstudiendaten, die in lme4 enthalten sind. Bates diskutiert dies in seinem Online- Buch über lme4. In Kapitel 3 betrachtet er zwei Modelle für die Daten.

und

Die Studie umfasste 18 Probanden, die über einen Zeitraum von 10 Tagen ohne Schlaf untersucht wurden. Die Reaktionszeiten wurden zu Beginn und an den folgenden Tagen berechnet. Es gibt einen deutlichen Effekt zwischen der Reaktionszeit und der Dauer des Schlafentzugs. Es gibt auch signifikante Unterschiede zwischen den Probanden. Modell A erlaubt die Möglichkeit einer Wechselwirkung zwischen zufälligen Intercept- und Slope-Effekten: Stellen Sie sich vor, dass Menschen mit schlechten Reaktionszeiten stärker unter den Auswirkungen von Schlafentzug leiden. Dies würde eine positive Korrelation in den Zufallseffekten implizieren.

In Bates 'Beispiel gab es keine offensichtliche Korrelation aus dem Gitterplot und keinen signifikanten Unterschied zwischen den Modellen. Um die oben gestellte Frage zu untersuchen, entschloss ich mich jedoch, die angepassten Werte der Schlafstudie zu verwenden, die Korrelation zu verbessern und die Leistung der beiden Modelle zu untersuchen.

Wie Sie auf dem Bild sehen können, sind lange Reaktionszeiten mit einem größeren Leistungsverlust verbunden. Die für die Simulation verwendete Korrelation betrug 0,58

Ich habe 1000 Proben mit der Simulationsmethode in lme4 simuliert, basierend auf den angepassten Werten meiner künstlichen Daten. Ich passte M0 und Ma aneinander an und schaute mir die Ergebnisse an. Der ursprüngliche Datensatz enthielt 180 Beobachtungen (10 für jede von 18 Personen), und die simulierten Daten hatten dieselbe Struktur.

Das Fazit ist, dass es sehr wenig Unterschiede gibt.

1. Die festen Parameter haben bei beiden Modellen exakt die gleichen Werte.
2. Die zufälligen Effekte sind leicht unterschiedlich. Für jede simulierte Stichprobe gibt es 18 zufällige Intercept- und 18 Slope-Effekte. Für jede Stichprobe müssen diese Effekte zu 0 addiert werden, was bedeutet, dass der mittlere Unterschied zwischen den beiden Modellen (künstlich) 0 beträgt. Die Varianzen und Kovarianzen unterscheiden sich jedoch. Die mittlere Kovarianz unter MA betrug 104 gegenüber 84 unter M0 (tatsächlicher Wert 112). Die Varianzen der Steigungen und Abschnitte waren unter M0 größer als unter MA, vermutlich um den zusätzlichen Spielraum für Wackelbewegungen zu erhalten, der ohne einen freien Kovarianzparameter benötigt wird.
3. Die ANOVA-Methode für lmer gibt eine F-Statistik zum Vergleichen des Slope-Modells mit einem Modell mit nur einem zufälligen Achsenabschnitt an (keine Auswirkung aufgrund von Schlafentzug). Offensichtlich war dieser Wert bei beiden Modellen sehr groß, aber er war bei MA typischerweise (aber nicht immer) größer (Mittelwert 62 gegenüber Mittelwert 55).
4. Die Kovarianz und Varianz der festen Effekte sind unterschiedlich.
5. Etwa die Hälfte der Zeit weiß es, dass MA richtig ist. Der mittlere p-Wert zum Vergleichen von M0 mit MA beträgt 0,0442. Trotz des Vorhandenseins einer aussagekräftigen Korrelation und 180 ausgewogener Beobachtungen würde das richtige Modell nur in etwa der Hälfte der Fälle ausgewählt.
6. Die vorhergesagten Werte unterscheiden sich bei beiden Modellen, jedoch nur geringfügig. Die mittlere Differenz zwischen den Vorhersagen beträgt 0 mit einem SD von 2,7. Der sd der vorhergesagten Werte selbst beträgt 60,9

Warum passiert das? @gung vermutete vernünftigerweise, dass das Versäumnis, die Möglichkeit einer Korrelation zu berücksichtigen, dazu führt, dass die zufälligen Effekte unkorreliert bleiben. Vielleicht sollte es so sein; Bei dieser Implementierung können die zufälligen Effekte jedoch korreliert werden, was bedeutet, dass die Daten die Parameter unabhängig vom Modell in die richtige Richtung ziehen können. Die Unrichtigkeit des falschen Modells zeigt sich in der Wahrscheinlichkeit, weshalb Sie die beiden Modelle auf dieser Ebene (manchmal) unterscheiden können. Das Mixed-Effects-Modell passt im Grunde genommen lineare Regressionen an jedes Subjekt an, die davon abhängen, wie das Modell sie für angebracht hält. Das falsche Modell erzwingt die Anpassung weniger plausibler Werte als das richtige Modell. Letztendlich werden die Parameter jedoch von der Anpassung an die tatsächlichen Daten bestimmt.

Hier ist mein etwas klobiger Code. Die Idee war, die Schlafstudiendaten anzupassen und dann einen simulierten Datensatz mit den gleichen Parametern, aber einer größeren Korrelation für die zufälligen Effekte zu erstellen. Dieser Datensatz wurde in simulate.lmer () eingespeist, um 1000 Proben zu simulieren, von denen jede in beide Richtungen angepasst wurde. Sobald ich angepasste Objekte gepaart hatte, konnte ich mithilfe von T-Tests oder was auch immer verschiedene Merkmale der Anpassung herausziehen und vergleichen.

    # Fit a model to the sleep study data, allowing non-zero correlation
fm01 <- lmer(Reaction ~ 1 + Days +(1+Days|Subject), data=sleepstudy, REML=FALSE)
# Now use this to build a similar data set with a correlation = 0.9
# Here is the covariance function for the random effects
# The variances come from the sleep study. The covariance is chosen to give a larger correlation
sigma.Subjects <- matrix(c(565.5,122,122,32.68),2,2) 
# Simulate 18 pairs of random effects
ranef.sim <- mvrnorm(18,mu=c(0,0),Sigma=sigma.Subjects)
# Pull out the pattern of days and subjects.
XXM <- model.frame(fm01) 
n <- nrow(XXM) # Sample size
# Add an intercept to the model matrix.
XX.f <- cbind(rep(1,n),XXM[,2])
# Calculate the fixed effects, using the parameters from the sleep study. 
yhat <- XX.f %*%  fixef(fm01 )
# Simulate a random intercept for each subject
intercept.r <- rep(ranef.sim[,1], each=10) 
# Now build the random slopes
slope.r <- XXM[,2]*rep(ranef.sim[,2],each=10)
# Add the slopes to the random intercepts and fixed effects
yhat2 <- yhat+intercept.r+slope.r
# And finally, add some noise, using the variance from the sleep study
y <- yhat2 + rnorm(n,mean=0,sd=sigma(fm01))
# Here is new "sleep study" data, with a stronger correlation.
new.data <- data.frame(Reaction=y,Days=XXM$Days,Subject=XXM$Subject)
# Fit the new data with its correct model
fm.sim <- lmer(Reaction ~ 1 + Days +(1+Days|Subject), data=new.data, REML=FALSE)
# Have a look at it
xyplot(Reaction ~ Days | Subject, data=new.data, layout=c(6,3), type=c("p","r"))
# Now simulate 1000 new data sets like new.data and fit each one
# using the right model and zero correlation model.
# For each simulation, output a list containing the fit from each and
# the ANOVA comparing them.
n.sim <- 1000
    sim.data <- vector(mode="list",)
    tempReaction <- simulate(fm.sim, nsim=n.sim)
    tempdata <- model.frame(fm.sim)
    for (i in 1:n.sim){
        tempdata$Reaction <- tempReaction[,i]
			output0 <- lmer(Reaction ~ 1 + Days +(1|Subject)+(0+Days|Subject), data = tempdata, REML=FALSE)
			output1 <- lmer(Reaction ~ 1 + Days +(Days|Subject), data=tempdata, REML=FALSE)
			temp <- anova(output0,output1)
			pval <- temp$`Pr(>Chisq)`[2]
        sim.data[[i]] <- list(model0=output0,modelA=output1, pvalue=pval)
    }

Placidia hat bereits eine gründliche Antwort mit simulierten Daten basierend auf dem sleepstudyDatensatz gegeben. Hier ist eine andere (weniger strenge) Antwort, die auch die sleepstudyDaten verwendet.

Wir sehen, dass man die geschätzte Korrelation zwischen dem zufälligen Schnittpunkt und der zufälligen Steigung beeinflussen kann, indem man die zufällige Prädiktorvariable "verschiebt". Schauen Sie sich die Ergebnisse von Modellen fm1und fm2unten an:

library(lmer)

#Fit Models
fm1 <- lmer(Reaction ~ Days + (Days | Subject), sleepstudy)
k <- 3 # Shift "Days" by an arbitrary amount
fm2 <- lmer(Reaction ~ I(Days + k) + (I(Days + k)| Subject), sleepstudy)

fm1 # Model Output
# Linear mixed model fit by REML ['lmerMod']
# Formula: Reaction ~ Days + (Days | Subject)
# Data: sleepstudy
# REML criterion at convergence: 1743.628
# Random effects:
#   Groups   Name        Std.Dev. Corr
# Subject  (Intercept) 24.740       
# Days         5.922   0.07
# Residual             25.592       
# Number of obs: 180, groups:  Subject, 18
# Fixed Effects:
#   (Intercept)         Days  
# 251.41        10.47

fm2 # Model Output
# Linear mixed model fit by REML ['lmerMod']
# Formula: Reaction ~ I(Days + k) + (I(Days + k) | Subject)
# Data: sleepstudy
# REML criterion at convergence: 1743.628
# Random effects:
#   Groups   Name        Std.Dev. Corr 
# Subject  (Intercept) 29.498        
# I(Days + k)  5.922   -0.55
# Residual             25.592        
# Number of obs: 180, groups:  Subject, 18
# Fixed Effects:
#   (Intercept)  I(Days + k)  
# 220.00        10.47

# Random effects from both models
cbind(ranef(fm1)$Subject,ranef(fm2)$Subject)
# (Intercept)        Days (Intercept) I(Days + k)
# 308   2.2585654   9.1989719 -25.3383538   9.1989727
# 309 -40.3985769  -8.6197032 -14.5394628  -8.6197043
# 310 -38.9602458  -5.4488799 -22.6136027  -5.4488807
# 330  23.6904985  -4.8143313  38.1334933  -4.8143315
# 331  22.2602027  -3.0698946  31.4698868  -3.0698946
# 332   9.0395259  -0.2721707   9.8560377  -0.2721706
# 333  16.8404311  -0.2236244  17.5113040  -0.2236243
# 334  -7.2325792   1.0745761 -10.4563076   1.0745761
# 335  -0.3336958 -10.7521591  31.9227854 -10.7521600
# 337  34.8903508   8.6282840   9.0054946   8.6282850
# 349 -25.2101104   1.1734142 -28.7303527   1.1734141
# 350 -13.0699567   6.6142050 -32.9125736   6.6142054
# 351   4.5778352  -3.0152572  13.6236077  -3.0152574
# 352  20.8635924   3.5360133  10.2555505   3.5360138
# 369   3.2754530   0.8722166   0.6588028   0.8722167
# 370 -25.6128694   4.8224646 -40.0802641   4.8224648
# 371   0.8070397  -0.9881551   3.7715053  -0.9881552
# 372  12.3145393   1.2840297   8.4624492   1.2840300

Aus der Modellausgabe sehen wir, dass sich die zufällige Varianzkorrelation geändert hat. Die Steigungen (fest und zufällig) blieben jedoch ebenso wie die Restvarianzschätzung gleich. Die (festen und zufälligen) Abschnittsschätzungen änderten sich als Reaktion auf die verschobene Variable.

De-Korrelation Zufalls Intercept-Slope - Kovarianz für LMMs ist Dr. Jack Weiss Skriptum diskutiert hier . Weiss merkt an, dass das Reduzieren der Varianzkorrelation auf diese Weise unter anderem manchmal zur Modellkonvergenz beitragen kann.

Das obige Beispiel variiert die zufällige Korrelation (Parameter "P5"). Aus der obigen Ausgabe geht hervor, dass der dritte Teil des OPs wie folgt behandelt wird:

#   Parameter           Status
=================================
P1  Fixed Intercept     Affected
P2  Random Intercepts   Affected
P3  Fixed Slope         Not Affected
P4  Random Slopes       Not Affected
P5  Random Correlation  Affected