Visualisierung gemischter Modellergebnisse

15

Eines der Probleme, die ich bei gemischten Modellen immer hatte, ist das Herausfinden von Datenvisualisierungen, wie sie auf einem Papier oder Poster landen können, sobald die Ergebnisse vorliegen.

Im Moment arbeite ich an einem Poisson-Mischeffektmodell mit einer Formel, die ungefähr so ​​aussieht:

a <- glmer(counts ~ X + Y + Time + (Y + Time | Site) + offset(log(people))

Mit etwas, das in glm () gepasst ist, könnte man leicht die predict () verwenden, um Vorhersagen für einen neuen Datensatz zu erhalten und daraus etwas aufzubauen. Aber wie würden Sie bei einer Ausgabe wie dieser so etwas wie eine grafische Darstellung der Rate über die Zeit mit den Verschiebungen von X (und wahrscheinlich mit einem festgelegten Wert von Y) erstellen? Ich denke, man könnte die Anpassung gut genug aus den Schätzungen der festen Effekte vorhersagen, aber was ist mit dem 95% -KI?

Gibt es noch etwas, das jemandem einfallen könnte, um die Ergebnisse zu visualisieren? Die Ergebnisse des Modells sind unten:

Random effects:
 Groups     Name        Variance   Std.Dev.  Corr          
 Site       (Intercept) 5.3678e-01 0.7326513               
            time        2.4173e-05 0.0049167  0.250        
            Y           4.9378e-05 0.0070270 -0.911  0.172 

Fixed effects:
              Estimate Std. Error z value Pr(>|z|)    
(Intercept) -8.1679391  0.1479849  -55.19  < 2e-16
X            0.4130639  0.1013899    4.07 4.62e-05
time         0.0009053  0.0012980    0.70    0.486    
Y            0.0187977  0.0023531    7.99 1.37e-15

Correlation of Fixed Effects:
         (Intr) Y    time
X      -0.178              
time    0.387 -0.305       
Y      -0.589  0.009  0.085
r data-visualization mixed-model Fomite
quelle
1
(+1) @EpiGrad: Warum machen Sie sich Sorgen über den CI (dh über den Standardfehler) der Vorhersagen aus dem Fixeffekt-Teil Ihres Modells?
Boscovich
1
@andrea Eine intellektuelle und eine praktische Antwort: Intellektuell bevorzuge ich im Allgemeinen die Quantifizierung und Visualisierung von Unsicherheit, wenn ich kann. Praktisch, weil ich mir ziemlich sicher bin, dass ein Rezensent danach fragen wird.
Fomite
Ja klar, aber ich meinte etwas anderes. Mein Kommentar war nicht klar genug, sorry. Sie schreiben in Ihrer Frage "aber was ist mit dem 95% CI?". Mein Kommentar lautet: Warum berechnen Sie nicht den Standardfehler der Vorhersage aus dem Teil des Modells mit festem Effekt? Wenn Sie in der Lage sind, die vorhergesagten Werte aus dem Festeffektteil zu berechnen, können Sie auch die SE und damit den CI berechnen. @EpiGrad
boscovich
@andrea Ah. Die Sorge ist, dass eines der Dinge, die ich vorhersagen möchte, Zeit, auch einen zufälligen Effekt hat, womit ich keine Ahnung habe, was ich tun soll.
Fomite
Nun, du willst es vorhersagen counts, nicht time. Sie beheben Werte X, Yund timeund mit dem festen Effekten Teil Ihres Modells Sie vorhersagen counts. Es ist wahr, dass timein Ihrem Modell auch ein Zufallseffekt enthalten ist (genau wie das Abfangen und Y), aber das spielt hier keine Rolle, da das Verwenden nur des Teils Ihres Modells mit festem Effekt für die Vorhersage so ist, als würden Sie die Zufallseffekte auf 0 setzen @EpiGrad
boscovich

Antworten:

4

Wenn Sie countsmit dem Fixeffekt-Teil Ihres Modells vorhersagen, setzen Sie die zufälligen Effekte auf Null (dh ihren Mittelwert). Dies bedeutet, dass Sie sie "vergessen" und Standardmaschinen verwenden können, um die Vorhersagen und die Standardfehler der Vorhersagen zu berechnen (mit denen Sie die Konfidenzintervalle berechnen können).

Dies ist ein Beispiel mit Stata, aber ich nehme an, es kann leicht in die R-Sprache "übersetzt" werden:

webuse epilepsy, clear
xtmepoisson seizures treat visit || subject: visit
predict log_seiz, xb
gen pred_seiz = exp(log_seiz)
predict std_log_seiz, stdp
gen ub = exp(log_seiz+invnorm(.975)*std_log_seiz)
gen lb = exp(log_seiz-invnorm(.975)*std_log_seiz)

tw (line pred_seiz ub lb visit if treat == 0, sort lc(black black black) ///
 lp(l - -)), scheme(s1mono) legend(off) ytitle("Predicted Seizures") ///
 xtitle("Visit")

Das Diagramm bezieht sich auf treat == 0und soll ein Beispiel sein (es visitist keine wirklich kontinuierliche Variable, aber es dient nur der Veranschaulichung). Die gestrichelten Linien geben 95% -Konfidenzintervalle an.

Bildbeschreibung hier eingeben

boscovich
quelle