Multilevel Models 1: Hierarchische Datenstrukturen
R-Klaus
1 Juli 2019
Szenario: Kosmetische Chirurgie und Lebensqualität
In einer medizinisch-psychologischen Studie wurde untersucht, ob Menschen, die sich einer nicht medizinisch notwendigen Operation unterzogen haben, zu einem bestimmten Zeitpunkt nach der OP eine höhere Lebensqualität/-zufriedenheit aufweisen als Personen auf der Warteliste (Kontrollgruppe).Dabei soll dafür kontrolliert werden, in welcher von 10 Kliniken die Personen ihre OP hatten bzw. auf der Warteliste standen.
Siehe Andy Field, Kapitel 19.
Folgende Hypothese wurden aufgestellt:
Ein Modell, das die Durchführung einer Operation und die durchführende Klinik beinhaltet, sagt die spätere Lebensqualität (PostQoL) von Personen vorher.
Datenstruktur erkennen und Multilevel-Modell aufstellen
Es liegt eine hierarchische Datenstruktur vor. Die Patient*innen und ihre PostQoL sind zehn Kliniken untergeordnet. Die Personen innerhalb einer Klinik werden sich in ihren Outcomes möglicherweise ähnlicher sein als sie es zu den Personen in anderen Kliniken sind. Unser Modell sollte diese Struktur berücksichtigen.
Ob eine Operation durchgeführt wurde oder nicht, wird als ganz normaler Prädiktor (fixed effect) ins Modell aufgenommen - genau so wie bisher bekannt. Der Einfluss der Klinik kann mit einem random intercept und/oder einem random slope eingefangen werden. Wir werden prüfen, welche dieser Komponenten nötig sind, um die Daten bestmöglich zu beschreiben.
Random intercept = der y-Achsenabschnitt der 10 Kliniken unterscheidet sich. Intercept = Wert des Kriteriums, wenn alle Prädiktoren gleich Null sind. Ein signifikantes random intercept würde hier also bedeuten, dass die PostQoL der Patient*innen sich auch zwischen den Kliniken unterscheidet, wenn gar keine Operation durchgeführt wurde. (Das könnte beispielsweise daran liegen, dass die Warteliste-Gruppe in manchen Kliniken besser betreut oder beraten wird als in anderen. Oder es liegt gar kein kausaler Zusammenhang vor, sondern manche Kliniken liegen geographisch einfach in wohlhabenderen Gegenden, in denen die Lebensqualität generell höher ist. Immerhin wurden die Personen den Kliniken ja vermutlich nicht randomisiert zugeteilt.)
Random slope = der Einfluss des fixed factors Operation unterscheidet sich zwischen den Kliniken: In manchen Kliniken unterscheiden sich die behandelte Gruppe und die Kontrollgruppe in ihrer späteren Lebensqualität vielleicht stärker als in anderen, in manchen vielleicht auch gar nicht. Oder: In manchen Kliniken geht die Veränderung in eine positive Richtung, in anderen sogar in eine negative. Das random slope könnte also die Qualität der verschiedenen Kliniken einfangen (wobei die möglicherweise resultierenden Unterschiede natürlich auch durch andere Faktoren erklärbar sein könnten, beispielsweise könnten manche Kliniken eher von Menschen besucht werden, die ohnehin gesünder sind oder ähnliches).
Ein Modell kann beide oder nur einen dieser random factors beinhalten. Bedeutung random intercept ohne random slope = Ohne Intervention ist die PostQoL zwischen den Kliniken unterschiedlich. Das Durchführen der OP führt dann aber in jeder Klinik in gleichem Maße zu einer Veränderung der PostQoL. Bedeutung random slope ohne random intercept = Ohne Intervention gibt es keine Unterschiede zwischen den Kliniken. Wenn eine OP durchgeführt wird, sind die Auswirkungen der OP auf die PostQoL aber je nach Klinik unterschiedlich. Welche random factors man ins Modell aufnimmt, hängt von den inhaltlichen Hypothesen ab. In unserem Fall sind beide Effekte plausibel und wir sollten daher prüfen, welche Komponenten unser Modell beinhalten sollte, um die Daten bestmöglich zu beschreiben.
Vorgehen:
Wir werden mittels mehrerer, schrittweiser Modellvergleiche prüfen, ob ein Modell mit random intercept und random slope unsere Daten vorhersagen kann:
Personenindex i = Patient*in
Gruppenindex j = Klinik
beta0,j: mittleres Intercept über Kliniken hinweg
u0,j: Standardabweichung der Intercepts für die verschiedenen Kliniken um das mittlere Intercept
beta1,j: mittlere Steigung für den Faktor Operation über die Kliniken hinweg
u1,j: Standardabweichung der Steigungen für die verschiedenen Kliniken
Statistische Hypothesen über Modell
Siehe Formelsammlung:
Schritte beim Vorgehen:
VL, Folie 26 (siehe darauffolgende Folien für noch mehr Details):
Einlesen der Daten, Laden der nötigen Pakete
data <- read.csv("surgery_data.csv",sep=";")neue Pakete ggf. installieren:
library(Rcmdr)
#install.packages("nlme",dependencies=T)
library(nlme) #für Berechnen des Modells
#install.packages("lmtest",dependencies=T)
library(lmtest) #für likelihood ratio tests
#install.packages("rcompanion",dependencies=T)
library(rcompanion) #für EffektgrößenGraphiken erstellen
Graphik ohne Krankenhäuser:
Erstmal könnte man sich ansehen, wie der Zusammenhang aussieht, wenn man die 10 Krankenhäuser nicht mit einbezieht.
RCmdr: Grafiken, Streudiagramm:
data$Surgery_label <- factor(data$Surgery_label)
with(data, plotMeans(Post_QoL, Surgery_label, error.bars="conf.int",
level=0.95, connect=TRUE))
Weiterhin könnten wir uns anschauen, wie die Mittelwerte aussehen, wenn man Clinic mit einbezieht:
data$Clinic_label <- factor(data$Clinic_label)
with(data, plotMeans(Post_QoL, Clinic_label, Surgery_label,
error.bars="conf.int", level=0.95, connect=TRUE, legend.pos="farright"))
Um eine Grafik zu erhalten, die ähnlich zu der Visualisierung in der VL ist (und optisch näher am linearen Modell), können wir auch ein Streudiagramm mit den folgenden Spezifikationen erstellen:
- x-Achse: surgery
- y-Achse PostQoL
- Gruppierungsfaktor Clinic
- Kleinst-Quadratlinie bei Optionen auswählen
scatterplot(Post_QoL~Surgery_num | Clinic_label, regLine=TRUE, smooth=FALSE,
boxplots=FALSE, by.groups=TRUE, data=data)## Warning in scatterplot.default(X[, 2], X[, 1], groups = X[, 3], xlab = xlab, : number of groups exceeds number of available colors
## colors are recycled
Hier sieht man, dass die Steigungen (zwischen Kontroll- und OP-Gruppe) für unterschiedliche Kliniken durchaus unterschiedlich aussehen. Weiterhin sieht man, dass auch in der Kontrollgruppe bereits Unterschiede zwischen Kliniken zu bestehen scheinen.
Schrittweiser Modellvergleich
Strategie:
Nullmodell aufstellen: Vergleich mit einem Modell mit random intercept. Zeigt, ob es baseline-Unterschiede zwischen Kliniken gibt und ein random intercept demnach angemessen ist.
Wenn Modell mit RI besser als Nullmodell: Aufstellen eines Modells mit random intercept plus fixed effect für den Prädiktor Operation, Vergleich mit dem Modell mit nur random intercept. Zeigt uns, ob der Faktor Operation über die Kliniken hinweg einen Einfluss hat.
Wenn Modell mit RI und fixed factor OP besser als Modell mit nur RI: Aufstellen eines Modells mit random intercept, Prädiktor plus random slope für Prädiktor. Vergleich mit dem vorherigen Modell. Zeigt uns, ob der Einfluss der Operation in den verschiedenen Kliniken unterschiedlich ist.
Schritt1: Nullmodell vs. Random Intercept:
Null <- gls(Post_QoL ~ 1, data = data, method="ML")
summary(Null)## Generalized least squares fit by maximum likelihood
## Model: Post_QoL ~ 1
## Data: data
## AIC BIC logLik
## 2017.124 2024.365 -1006.562
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 59.60978 0.5596972 106.5036 0
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -2.1127754 -0.7875625 -0.1734394 0.7962286 3.0803354
##
## Residual standard error: 9.281527
## Degrees of freedom: 276 total; 275 residualGeneralized least squares = auf kleinsten Quadraten basierende Methode, Regressionsparameter auch dann effizient zu schätzen, wenn Abhängigkeiten in den Daten vorliegen (zb. hierarchische Struktur).
Maximum likelihood: Schätzmethode für die Parameter. Es werden diejenigen Parameter gesucht, unter denen die erhaltenen Daten am wahrscheinlichsten sind.
Random intercept-Modell (Achtung: jetzt mit lme-Funktion):
(Es ist für die Berechnung hier übrigens egal, ob man die numerische Clinic-Variable oder den Clinic-Faktor (mit labels) einsetzt.)
RI <- lme(Post_QoL ~ 1, data= data, random = ~1|Clinic_num, method="ML")
summary(RI)## Linear mixed-effects model fit by maximum likelihood
## Data: data
## AIC BIC logLik
## 1911.473 1922.334 -952.7364
##
## Random effects:
## Formula: ~1 | Clinic_num
## (Intercept) Residual
## StdDev: 5.909691 7.238677
##
## Fixed effects: Post_QoL ~ 1
## Value Std.Error DF t-value p-value
## (Intercept) 60.08377 1.923283 266 31.24022 0
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -1.8828507 -0.7606631 -0.1378732 0.7075242 2.8607949
##
## Number of Observations: 276
## Number of Groups: 10Unter “random effects” erhalten wir jetzt die Standardabweichung des Intercepts für die verschiedenen Kliniken. (Und die der Residuen.)
Wie tippt man den vertikalen Strich ein (der ist links neben “Y”)?
Das AIC des zweiten Modells ist kleiner als das des ersten, was eine bessere Passung auf die Daten anzeigt. Wir prüfen nun, ob diese Verbesserung signifikant ist:
lrtest(Null,RI)## Warning in modelUpdate(objects[[i - 1]], objects[[i]]): original model was of
## class "gls", updated model is of class "lme"## Likelihood ratio test
##
## Model 1: Post_QoL ~ 1
## Model 2: Post_QoL ~ 1
## #Df LogLik Df Chisq Pr(>Chisq)
## 1 2 -1006.56
## 2 3 -952.74 1 107.65 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Die Verbesserung ist signifikant: wir sollten den random intercept beibehalten.
Die Richtung des Unterschieds (welches Modell jetzt besser ist) geht aus dem lrtest-Output selbst nicht hervor. Das müssen wir am AIC der Modelle ablesen (niedriger = besser).
Schritt 2: RI-Modell gegen RI-Modell plus fixed factor (unser Prädiktor):
RI_pred <- lme(Post_QoL ~ Surgery_num, data = data, random = ~1|Clinic_num, method="ML")
summary(RI_pred)## Linear mixed-effects model fit by maximum likelihood
## Data: data
## AIC BIC logLik
## 1910.137 1924.619 -951.0686
##
## Random effects:
## Formula: ~1 | Clinic_num
## (Intercept) Residual
## StdDev: 6.099513 7.18542
##
## Fixed effects: Post_QoL ~ Surgery_num
## Value Std.Error DF t-value p-value
## (Intercept) 59.30517 2.0299632 265 29.21490 0.000
## Surgery_num 1.66583 0.9091314 265 1.83233 0.068
## Correlation:
## (Intr)
## Surgery_num -0.21
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -1.8904290 -0.7191399 -0.1420998 0.7177762 2.8644538
##
## Number of Observations: 276
## Number of Groups: 10Das AIC ist im Vergleich zum vorheigen Modell leicht gesunken, der Koeffizient für den fixed factor Surgery zudem (knapp) nicht signifikant.
Test der beiden Modelle gegeneinander:
lrtest(RI,RI_pred)## Likelihood ratio test
##
## Model 1: Post_QoL ~ 1
## Model 2: Post_QoL ~ Surgery_num
## #Df LogLik Df Chisq Pr(>Chisq)
## 1 3 -952.74
## 2 4 -951.07 1 3.3356 0.06779 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Sie unterscheiden sich nicht signifikant. Über die verschiedenen Kliniken hinweg scheint sich die Durchführung einer Operation also nicht auf die spätere Lebensqualität auszuwirken.
Es könnte aber sein, dass sich die Effekte zwischen den Kliniken so unterscheiden, dass sie sich, wenn man über alle Kliniken mittelt, ausgleichen. Ein solches Muster können wir nur über ein random slope einfangen.
Schritt 3: Modell mit RI und fixed factor Operation gegen Modell mit RI, fixed factor Operation und random slope für Operation:
RI_pred_RS <- lme(Post_QoL ~ Surgery_num, data = data, random = ~Surgery_num|Clinic_num, method="ML")
summary(RI_pred_RS)## Linear mixed-effects model fit by maximum likelihood
## Data: data
## AIC BIC logLik
## 1838.972 1860.694 -913.4859
##
## Random effects:
## Formula: ~Surgery_num | Clinic_num
## Structure: General positive-definite, Log-Cholesky parametrization
## StdDev Corr
## (Intercept) 8.645944 (Intr)
## Surgery_num 8.063630 -0.947
## Residual 6.101323
##
## Fixed effects: Post_QoL ~ Surgery_num
## Value Std.Error DF t-value p-value
## (Intercept) 59.38746 2.794293 265 21.25312 0.0000
## Surgery_num 0.36775 2.675365 265 0.13746 0.8908
## Correlation:
## (Intr)
## Surgery_num -0.927
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -2.94484564 -0.66035532 -0.05116588 0.66293618 2.58513956
##
## Number of Observations: 276
## Number of Groups: 10Unter “random effects” finden sich nun die geschätzten Standardabweichungen für Intercept und Steigung.
Das AIC ist deutlich gesunken.
Wir prüfen, ob das Modell mit random slope signifikant besser ist als das vorherige:
lrtest(RI_pred,RI_pred_RS)## Likelihood ratio test
##
## Model 1: Post_QoL ~ Surgery_num
## Model 2: Post_QoL ~ Surgery_num
## #Df LogLik Df Chisq Pr(>Chisq)
## 1 4 -951.07
## 2 6 -913.49 2 75.165 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Das Modell mit random intercept und random slope beschreibt unsere Daten also am besten.
Die statistische Hypothese kann angenommen werden.
Effektgröße
Hier kann eine Annäherung an R-Squared berechnet werden, bei der das final gewählte Modell mit dem Nullmodell verglichen wird:
nagelkerke(RI_pred_RS,Null)## $Models
##
## Model: "lme.formula, Post_QoL ~ Surgery_num, data, ~Surgery_num | Clinic_num, ML"
## Null: "gls, Post_QoL ~ 1, data, ML"
##
## $Pseudo.R.squared.for.model.vs.null
## Pseudo.R.squared
## McFadden 0.0924695
## Cox and Snell (ML) 0.4905720
## Nagelkerke (Cragg and Uhler) 0.4909050
##
## $Likelihood.ratio.test
## Df.diff LogLik.diff Chisq p.value
## -4 -93.076 186.15 3.5558e-39
##
## $Number.of.observations
##
## Model: 276
## Null: 276
##
## $Messages
## [1] "Note: For models fit with REML, these statistics are based on refitting with ML"
##
## $Warnings
## [1] "None"Cox & Snell oder Nagelkerke können als R^2-Annäherung interpretiert werden.
Exkurs: Was passiert, wenn man mit den Daten einfach eine Anova oder eine Regression rechnet?
Anova
Zweifaktoriell über Commander
AnovaModel.1 <- lm(Post_QoL ~ Clinic_label*Surgery_label, data=data,
contrasts=list(Clinic_label ="contr.Sum", Surgery_label
="contr.Sum"))
Anova(AnovaModel.1)## Anova Table (Type II tests)
##
## Response: Post_QoL
## Sum Sq Df F value Pr(>F)
## Clinic_label 10016.6 9 29.9434 < 2e-16 ***
## Surgery_label 205.8 1 5.5361 0.01939 *
## Clinic_label:Surgery_label 4216.1 9 12.6037 < 2e-16 ***
## Residuals 9515.1 256
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Wir würden hier also zu einem ähnlichen Schluss kommen. Was dahintersteckt, ist aber aufwändiger als ein Multilevelmodell, da im Hintergrund ja Dummy-Variablen für den Faktor Klinik angelegt werden. Hier wären das 9 Dummy-Variablen, aber wenn es noch mehr Gruppen gäbe, könnten es auch deutlich mehr werden. Bei einem Multilevelmodell würde bei diesem Beispiel immer dieselbe Anzahl von Parametern geschätzt (die Mittelwerte und Streuungen für Intercept und Slope, wenn man sie als random effects einsetzt), auch wenn wir 500 Kliniken vergleichen würden.
Regression
Achtung: hier muss zumindest für Clinic der Faktor als Prädiktor gewählt werden, sonst macht der Output keinen Sinn (die Zahlen 1-10 würden als kontinuierliche Variable verstanden. Bei surgery, wo es nur 0 und 1 gibt, macht das keinen Unterschied, der gleiche Gruppenunterschied würde in beiden Fällen ausgegeben. Der Konsistenz halber wird hier allerdings ebenfalls der Faktor verwendet.)
summary(lm(Post_QoL ~ Surgery_label*Clinic_label, data = data))##
## Call:
## lm(formula = Post_QoL ~ Surgery_label * Clinic_label, data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -18.2941 -3.9765 -0.2529 3.8766 16.6133
##
## Coefficients:
## Estimate Std. Error t value
## (Intercept) 60.4000 3.0483 19.814
## Surgery_labelWaiting List 8.7235 3.3880 2.575
## Clinic_labelClinic10 -0.5182 3.5596 -0.146
## Clinic_labelClinic2 8.2000 4.0897 2.005
## Clinic_labelClinic3 -3.1133 3.4307 -0.907
## Clinic_labelClinic4 -2.4615 3.4859 -0.706
## Clinic_labelClinic5 -3.0400 3.4307 -0.886
## Clinic_labelClinic6 2.6688 3.4081 0.783
## Clinic_labelClinic7 -1.4833 3.3700 -0.440
## Clinic_labelClinic8 -5.4235 3.3880 -1.601
## Clinic_labelClinic9 1.0176 3.3880 0.300
## Surgery_labelWaiting List:Clinic_labelClinic10 -19.6103 4.0885 -4.796
## Surgery_labelWaiting List:Clinic_labelClinic2 -3.3185 4.5575 -0.728
## Surgery_labelWaiting List:Clinic_labelClinic3 -0.7402 4.2040 -0.176
## Surgery_labelWaiting List:Clinic_labelClinic4 1.6321 4.0650 0.402
## Surgery_labelWaiting List:Clinic_labelClinic5 -6.9335 4.2040 -1.649
## Surgery_labelWaiting List:Clinic_labelClinic6 -15.6692 4.0817 -3.839
## Surgery_labelWaiting List:Clinic_labelClinic7 -19.3873 3.9661 -4.888
## Surgery_labelWaiting List:Clinic_labelClinic8 -10.5200 4.1691 -2.523
## Surgery_labelWaiting List:Clinic_labelClinic9 -18.7048 4.1284 -4.531
## Pr(>|t|)
## (Intercept) < 2e-16 ***
## Surgery_labelWaiting List 0.010591 *
## Clinic_labelClinic10 0.884375
## Clinic_labelClinic2 0.046013 *
## Clinic_labelClinic3 0.365007
## Clinic_labelClinic4 0.480738
## Clinic_labelClinic5 0.376393
## Clinic_labelClinic6 0.434316
## Clinic_labelClinic7 0.660196
## Clinic_labelClinic8 0.110653
## Clinic_labelClinic9 0.764140
## Surgery_labelWaiting List:Clinic_labelClinic10 2.74e-06 ***
## Surgery_labelWaiting List:Clinic_labelClinic2 0.467188
## Surgery_labelWaiting List:Clinic_labelClinic3 0.860377
## Surgery_labelWaiting List:Clinic_labelClinic4 0.688380
## Surgery_labelWaiting List:Clinic_labelClinic5 0.100316
## Surgery_labelWaiting List:Clinic_labelClinic6 0.000156 ***
## Surgery_labelWaiting List:Clinic_labelClinic7 1.80e-06 ***
## Surgery_labelWaiting List:Clinic_labelClinic8 0.012232 *
## Surgery_labelWaiting List:Clinic_labelClinic9 9.02e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 6.097 on 256 degrees of freedom
## Multiple R-squared: 0.5998, Adjusted R-squared: 0.5701
## F-statistic: 20.19 on 19 and 256 DF, p-value: < 2.2e-16Die dummykodierten Kliniken und vor allem die vielen Interaktionsvariablen sind nicht besonders leicht zu interpretieren. Daher funktioniert das Vorgehen zwar, aber ist an dieser Stelle nicht besonders informativ (und davon abgesehen gilt dasselbe Argument bzgl. der potentiell vielen Dummy-Variablen wie oben.)
PS.:
Nicolaas Jan Dirk Nagelkerke wünscht alles Gute für die Klausur!
