Übungszettel Multi-Level-Modelle

M.Psy.205, Dozent: Peter Zezula

German

Links

Übungszettel als PDF-Datei zum Drucken

Übungszettel mit Lösungen

Lösungszettel als PDF-Datei zum Drucken

Der gesamte Übugszettel als .Rmd-Datei (Zum Downloaden: Rechtsklick > Speichern unter…)

Hinweise zur Bearbeitung

1. Bitte beantworten Sie die Fragen in einer .Rmd Datei. Sie können Sie über Datei > Neue Datei > R Markdown... eine neue R Markdown Datei erstellen. Den Text unter dem Setup Chunk (ab Zeile 11) können Sie löschen. Unter diesem Link können Sie auch unsere Vorlage-Datei herunterladen (Rechtsklick > Speichern unter…).
2. Informationen, die Sie für die Bearbeitung benötigen, finden Sie auf der Website der Veranstaltung
3. Zögern Sie nicht, im Internet nach Lösungen zu suchen. Das effektive Suchen nach Lösungen für R-Probleme im Internet ist tatsächlich eine sehr nützliche Fähigkeit, auch Profis arbeiten auf diese Weise. Die beste Anlaufstelle dafür ist der R-Bereich der Programmiererplattform Stackoverflow
4. Auf der Website von R Studio finden Sie sehr hilfreiche Übersichtszettel zu vielen verschiedenen R-bezogenen Themen. Ein guter Anfang ist der Base R Cheat Sheet

Ressourcen

Das Thema Multi-Level-Modelle wird im Field, Kapitel 19 sehr gut verständlich und in adäquater Tiefe behandelt. Starke Empfehlung!

Allgemeines

In Sheet 07 haben Sie Multi-Level-Ansätze bereits im Kontext von Messwiederholungen kennengelernt. Dabei wurden verschiedene Messzeitpunkte auf Level 1 innerhalb verschiedener Personen bzw. Länder auf Level 2 genested/verschachtelt begriffen. Da sich Datenpunkte der gleichen Person etc. vermutlich ähnlicher sind als Datenpunkte unterschiedlicher Personen, ist die Annahme der Unabhängingkeit der Datenpunkte (a.k.a. keine Autokorrelation), die z.B. normale Regressionsverfahren beinhalten, verletzt. Multi-Level-Modelle ermöglichen es, die durch die hierarchische Struktur eines Datensatz entstehende Autokorrelation zu berücksichtigen. Abgesehen von den bereits bekannten Anwendungen auf Personen und Länder, werden in Beispielen häufig Kliniken (wie in den Aufgaben in diesem Sheet) und Schulklassen auf Level 2 berücksichtigt.

1 Daten und Pakete

Laden Sie zunächst das tidyverse und das Paket nlme. Lesen Sie dann den Datensatz https://md.psych.bio.uni-goettingen.de/mv/data/div/cosmetic_surgery.dat (Tabstopp-getrennt) ein. Dieser enthält Daten zu der Lebensqualität von Personen, die sich für eine Schönheitsoperation entschieden haben:

Variable	Erklärung
particnu	Teilnehmernummer
Post_QoL	Lebensqualität nach Eingriff bzw. Wartezeit
Base_QoL	Baseline-Lebensqualität vor Beginn
Clinic	Codenummer dafür, zu welcher der zehn untersuchten Kliniken der Patient gehörte
Surgery	1 = Operation erfolgt, 0 = Wartelistenkontrollgruppe
Reason	1 = Es bestand ein körperlicher Grund für den Eingriff, 0 = Der Eingriff erfolgte zur Aussehensverbesserung
Age	Alter der Probanden
Gender	1 = Männlich, 0 = Weiblich
BDI	Standardisierter Fragebogen zur Erhebung depressiver Störungen
_Text-Variablen	Wiederholungen von Surgery, Reason und Gender in Worten, für etwaige Plotlabels

Uns interessiert vor allem, ob die Lebensqualität nach einer Schönheitsoperation tatsächlich höher ist als nach einer Zeit auf der Warteliste. Baseline-Lebensqualität, Alter, Geschlecht und Depressivität wurden als Kovariate erhoben. Clinic als Level-2-Variable erlaubt uns, für etwaige Variabilität in der OP-Güte zu kontrollieren, und Reason wird uns später eine differenziertere Aussage zum OP-Effekt erlauben.

2 - Fixed Effects only

Erzeugen Sie zunächst ein Modell ohne Random Effects, in dem Sie nur die Kovariaten als Prädiktoren verwenden. Nutzen Sie hierfür anstelle des lm()-Befehls den nlme::gls()-Befehl (Die Syntax ist die gleiche, sie müssen nur zusätzlich method = "ML" spezifizieren). Dies hat den Vorteil, dass die so erzeugten Modelle mittels anova() direkt mit Modellen mit Random Effects verglichen werden können.
Erzeugen Sie nun ein weiteres Modell ohne Random Effects, in dem neben den Kovariaten auch Surgery als Prädiktor verwendet wird. Wenn Sie möchten, versuchen Sie, das neue Modell mithilfe des update()-Befehls aus dem ersten Modell herzuleiten. Wenn das nicht klappt, können Sie das Modell aber auch wie gewohnt aufstellen.
Vergleichen Sie die beiden Modelle mittels anova(). Betrachten Sie zusätzlich die summary() des zweiten Modells, und dort insbesondere den Surgery-Prädiktor.

3 - Random Effects

Es ist gut vorstellbar, dass die Qualität der durchgeführten Operationen zwischen den untersuchten Kliniken schwankt, z.B. aufgrund der persönlichen Fähigkeiten der dort angestellten Chirurgen. Erstellen Sie ein Modell mit den gleichen Prädiktoren wie das zweite Modell aus der Aufgabe zuvor. Erlauben Sie hier jedoch, dass der Effekt von Surgery zwischen den verschiedenen Kliniken schwankt. Versuchen Sie, sich die Syntax hierfür mittels ?nlme::lme() selbst zu erarbeiten. Ansonsten finden Sie in “sheet_repeated_measure” eine Erklärung.
Zur Klärung, ob Random Effects angemessen sind, wird häufig mittels nlme::intervals(Modellname) das 95%-Konfidenzintervall der Streuung der Effekte untersucht. Führen Sie diese Untersuchung durch. Sprechen die Ergebnisse eher für oder gegen Random Effects? Warum?
Nutzen Sie wieder eine anova(), um zu prüfen, ob das neue Modell signifikant mehr Varianz aufklärt als das letzte. Betrachten Sie auch wieder die summary(). Was fällt Ihnen auf?

4 - Interaktion 1

Bisher haben wir die Information, aus welchen Gründen die untersuchten Personen eine Schönheitsoperation anstrebten, ignoriert. Es ist jedoch gut denkbar, dass eine Schönheitsoperation, die ein körperliches Problem behandeln soll (etwa eine Brustverkleinerung zur Behandlung von Rückenschmerzen), einen anderen Einfluss auf die Lebensqualität hat, als eine Schönheitsoperation, die lediglich der Verbesserung des subjektiven Aussehens dient. Erweitern Sie Ihr Modell aus der vorigen Aufgabe daher um eine Interaktion von Surgery und Reason. Führen Sie wie gehabt einen Modellvergleich mit dem vorigen Modell durch, und betrachten Sie die summary des neuen Modells.

5 - Interaktion 2

Das Ergebnis der vorigen Aufgabe legt nahe, dass der Effekt einer Schönheitsoperation von den Gründen für die OP abhängt. Erzeugen Sie zur näheren Untersuchung zwei weitere Modelle mit Random Effects, nur diesmal getrennt für die zwei möglichen Gründe. Die Interaktion von Surgery und Reason sollte nicht im Modell auftauchen, da es ja pro Modell nur einen Wert von Reason gibt. Versuchen Sie dies zunächst, indem Sie sich in die subset-Option von nlme::lme() einlesen. Wenn dies nicht klappt, können Sie auch mit dplyr::filter() zunächst zwei neue Datensätze erstellen. Betrachten Sie beide Modelle, und interpretieren Sie sie inhaltlich.

6 - Random Effects vs Interaktion

Sowohl bei Random Effects, als auch bei Interaktionen schwankt die Wirkung einer Variablen in Abhängigkeit von der Ausprägung einer anderen Variable. In dem hier verwendeten Beispiel war der Effekt von Surgery sowohl unterschiedlich für verschiedene Kliniken, als auch für verschiedene Gründe für die OP. Können Sie sich vorstellen, was der inhaltliche Unterschied ist? Bei welcher Fragestellung wäre eine Interaktion zwischen Klinik und Surgery möglicherweise die interessantere Herangehensweise? Und warum erlauben wir zwar dem Effekt von Surgery, zwischen den Kliniken zu schwanken, nicht aber den Effekten von Alter, Geschlecht, Depressivität etc.?

7 - Plot

Überlegen Sie sich eine angemessene Visualisierung für die Erkenntnisse aus dieser Auswertung!

Lösung

Aufgabe 1

library(tidyverse)
library(nlme)
cosmetic_surgery <- read_delim("https://md.psych.bio.uni-goettingen.de/mv/data/div/cosmetic_surgery.dat", "\t")

## Rows: 276 Columns: 12
## ── Column specification ──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
## Delimiter: "\t"
## chr (3): Surgery_Text, Reason_Text, Gender_Text
## dbl (9): particnu, Post_QoL, Base_QoL, Clinic, Surgery, Reason, Age, Gender, BDI
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

# alternativ
cosmetic_surgery <- readr::read_tsv("https://md.psych.bio.uni-goettingen.de/mv/data/div/cosmetic_surgery.dat")

## Rows: 276 Columns: 12
## ── Column specification ──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
## Delimiter: "\t"
## chr (3): Surgery_Text, Reason_Text, Gender_Text
## dbl (9): particnu, Post_QoL, Base_QoL, Clinic, Surgery, Reason, Age, Gender, BDI
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

Aufgabe 2

m_covariates <- nlme::gls(Post_QoL ~ Base_QoL + Age + Gender + BDI,
                          data = cosmetic_surgery,
                          method = "ML")

m_fixed_effects <- update(m_covariates, .~. + Surgery)
# m_fixed_effects <- nlme::gls(Post_QoL ~ Base_QoL + Age + Gender + BDI + Surgery,
#       data = cosmetic_surgery,
#       method = "ML")

anova(m_covariates, m_fixed_effects)

##                 Model df      AIC      BIC    logLik   Test  L.Ratio p-value
## m_covariates        1  6 1811.000 1832.723 -899.5001                        
## m_fixed_effects     2  7 1807.024 1832.367 -896.5122 1 vs 2 5.975713  0.0145

summary(m_fixed_effects)

## Generalized least squares fit by maximum likelihood
##   Model: Post_QoL ~ Base_QoL + Age + Gender + BDI + Surgery 
##   Data: cosmetic_surgery 
##        AIC      BIC    logLik
##   1807.024 1832.367 -896.5122
## 
## Coefficients:
##                 Value Std.Error   t-value p-value
## (Intercept) 16.674777 2.6650283  6.256885  0.0000
## Base_QoL     0.507214 0.0456353 11.114509  0.0000
## Age          0.231408 0.0537485  4.305381  0.0000
## Gender      -1.016271 1.3183664 -0.770856  0.4415
## BDI          0.130171 0.0367314  3.543866  0.0005
## Surgery     -1.956844 0.8049688 -2.430957  0.0157
## 
##  Correlation: 
##          (Intr) Bas_QL Age    Gender BDI   
## Base_QoL -0.829                            
## Age      -0.170 -0.263                     
## Gender    0.116  0.011 -0.695              
## BDI       0.110 -0.160 -0.522  0.687       
## Surgery  -0.021 -0.022 -0.136 -0.079  0.082
## 
## Standardized residuals:
##         Min          Q1         Med          Q3         Max 
## -2.02107403 -0.75479722 -0.02969699  0.66683065  3.73213879 
## 
## Residual standard error: 6.229488 
## Degrees of freedom: 276 total; 270 residual

#zweites Modell signifikant bessere Varianzaufklärung,
#Surgery signifikanter Prädiktor

Aufgabe 3

m_random <- nlme::lme(Post_QoL ~ Base_QoL + Age + Gender + BDI + Surgery,
                      data = cosmetic_surgery,
                      random = ~Surgery | Clinic,
                      method = "ML")

nlme::intervals(m_random)

## Approximate 95% confidence intervals
## 
##  Fixed effects:
##                   lower       est.      upper
## (Intercept) 27.33892415 34.1257908 40.9126574
## Base_QoL     0.10270956  0.2031965  0.3036834
## Age          0.18480652  0.2835955  0.3823844
## Gender      -2.98124310 -0.7258946  1.5294539
## BDI          0.03362661  0.1012212  0.1688157
## Surgery     -3.14643520 -0.9008826  1.3446700
## 
##  Random Effects:
##   Level: Clinic 
##                               lower       est.       upper
## sd((Intercept))           2.9765969  4.8911584  8.03717504
## sd(Surgery)               1.2816426  2.7914417  6.07981234
## cor((Intercept),Surgery) -0.9811474 -0.8268695 -0.02891568
## 
##  Within-group standard error:
##    lower     est.    upper 
## 4.911125 5.356450 5.842157

# Weder das Konfidenzinterval von sd((Intercept)), noch von sd(Surgery) beinhaltet
# die Null, insofern kann von signifikanter Schwankung der Effekte
# und somit gerechtfertigten Random Effects ausgegangen werden.

anova(m_fixed_effects, m_random)

##                 Model df      AIC      BIC    logLik   Test  L.Ratio p-value
## m_fixed_effects     1  7 1807.024 1832.367 -896.5122                        
## m_random            2 10 1763.483 1799.687 -871.7414 1 vs 2 49.54154  <.0001

summary(m_random)

## Linear mixed-effects model fit by maximum likelihood
##   Data: cosmetic_surgery 
##        AIC      BIC    logLik
##   1763.483 1799.687 -871.7414
## 
## Random effects:
##  Formula: ~Surgery | Clinic
##  Structure: General positive-definite, Log-Cholesky parametrization
##             StdDev   Corr  
## (Intercept) 4.891158 (Intr)
## Surgery     2.791442 -0.827
## Residual    5.356450       
## 
## Fixed effects:  Post_QoL ~ Base_QoL + Age + Gender + BDI + Surgery 
##                Value Std.Error  DF   t-value p-value
## (Intercept) 34.12579  3.484780 261  9.792811  0.0000
## Base_QoL     0.20320  0.051596 261  3.938226  0.0001
## Age          0.28360  0.050724 261  5.590941  0.0000
## Gender      -0.72589  1.158030 261 -0.626836  0.5313
## BDI          0.10122  0.034707 261  2.916443  0.0038
## Surgery     -0.90088  1.153000 261 -0.781338  0.4353
##  Correlation: 
##          (Intr) Bas_QL Age    Gender BDI   
## Base_QoL -0.782                            
## Age      -0.109 -0.259                     
## Gender    0.086 -0.017 -0.619              
## BDI      -0.050 -0.004 -0.485  0.631       
## Surgery  -0.231 -0.074 -0.094 -0.041  0.045
## 
## Standardized Within-Group Residuals:
##        Min         Q1        Med         Q3        Max 
## -2.0491309 -0.7468661 -0.1504611  0.6704928  2.9847839 
## 
## Number of Observations: 276
## Number of Groups: 10

#Das neue Modell klärt signifikant mehr Varianz auf, jedoch ist
#der Fixed Effect von Surgery nun nicht mehr signifikant!

Aufgabe 4

m_interaction <- update(m_random, .~. + Surgery:Reason)
# m_interaction <- nlme::lme(Post_QoL ~ Base_QoL +
#                              Age +
#                              Gender +
#                              BDI +
#                              Surgery +
#                              Surgery:Reason,
#                            data = cosmetic_surgery,
#                            random = ~Surgery | Clinic,
#                            method = "ML")
anova(m_random, m_interaction)

##               Model df      AIC      BIC    logLik   Test  L.Ratio p-value
## m_random          1 10 1763.483 1799.687 -871.7414                        
## m_interaction     2 11 1748.482 1788.306 -863.2410 1 vs 2 17.00082  <.0001

summary(m_interaction)

## Linear mixed-effects model fit by maximum likelihood
##   Data: cosmetic_surgery 
##        AIC      BIC   logLik
##   1748.482 1788.306 -863.241
## 
## Random effects:
##  Formula: ~Surgery | Clinic
##  Structure: General positive-definite, Log-Cholesky parametrization
##             StdDev   Corr  
## (Intercept) 4.463083 (Intr)
## Surgery     1.229488 -0.812
## Residual    5.231446       
## 
## Fixed effects:  Post_QoL ~ Base_QoL + Age + Gender + BDI + Surgery + Surgery:Reason 
##                    Value Std.Error  DF   t-value p-value
## (Intercept)    31.324071  3.273622 260  9.568629  0.0000
## Base_QoL        0.222262  0.050330 260  4.416098  0.0000
## Age             0.294238  0.047655 260  6.174378  0.0000
## Gender         -1.248325  1.135014 260 -1.099832  0.2724
## BDI             0.155575  0.035984 260  4.323450  0.0000
## Surgery        -4.484673  1.152615 260 -3.890868  0.0001
## Surgery:Reason  5.819623  1.312833 260  4.432874  0.0000
##  Correlation: 
##                (Intr) Bas_QL Age    Gender BDI    Surgry
## Base_QoL       -0.798                                   
## Age            -0.033 -0.304                            
## Gender          0.114 -0.029 -0.639                     
## BDI            -0.078  0.017 -0.505  0.540              
## Surgery         0.007 -0.120 -0.034  0.053 -0.201       
## Surgery:Reason -0.114  0.067 -0.068 -0.124  0.353 -0.709
## 
## Standardized Within-Group Residuals:
##        Min         Q1        Med         Q3        Max 
## -2.2788267 -0.7482315 -0.1017572  0.6956910  3.0657904 
## 
## Number of Observations: 276
## Number of Groups: 10

#Neues Modell signifikant besser, Interaktion signifikanter Prädiktor.
#Interessanterweise wird auch Surgery selbst wieder signifikant!

Aufgabe 5

physical_subset <- cosmetic_surgery$Reason == 1
appearance_subset <- !physical_subset #! bedeutet "nicht"

m_physical <- update(m_random, subset = physical_subset)
m_appearance <- update(m_random, subset = appearance_subset)

# Alternative ohne subset
# df_physical <- dplyr::filter(cosmetic_surgery, Reason == 1)
# df_appearance <- dplyr::filter(cosmetic_surgery, Reason == 0)
# 
# m_physical_2 <- update(m_random, data = df_physical)
# m_appearance_2 <- update(m_random, data = df_appearance)

summary(m_physical)

## Linear mixed-effects model fit by maximum likelihood
##   Data: cosmetic_surgery 
##   Subset: physical_subset 
##        AIC      BIC    logLik
##   1154.884 1186.702 -567.4421
## 
## Random effects:
##  Formula: ~Surgery | Clinic
##  Structure: General positive-definite, Log-Cholesky parametrization
##             StdDev   Corr  
## (Intercept) 4.679271 (Intr)
## Surgery     1.668346 -0.794
## Residual    5.471991       
## 
## Fixed effects:  Post_QoL ~ Base_QoL + Age + Gender + BDI + Surgery 
##                 Value Std.Error  DF   t-value p-value
## (Intercept) 29.893046  4.411360 163  6.776379  0.0000
## Base_QoL     0.265651  0.069800 163  3.805898  0.0002
## Age          0.274836  0.066637 163  4.124374  0.0001
## Gender      -0.460955  1.502015 163 -0.306891  0.7593
## BDI          0.118640  0.064379 163  1.842838  0.0672
## Surgery      0.666083  1.162912 163  0.572772  0.5676
##  Correlation: 
##          (Intr) Bas_QL Age    Gender BDI   
## Base_QoL -0.834                            
## Age      -0.067 -0.283                     
## Gender    0.147 -0.079 -0.613              
## BDI      -0.184  0.087 -0.381  0.518       
## Surgery  -0.036 -0.055 -0.209 -0.111  0.005
## 
## Standardized Within-Group Residuals:
##        Min         Q1        Med         Q3        Max 
## -2.3720422 -0.7824301 -0.1276846  0.6989338  2.9496881 
## 
## Number of Observations: 178
## Number of Groups: 10

summary(m_appearance)

## Linear mixed-effects model fit by maximum likelihood
##   Data: cosmetic_surgery 
##   Subset: appearance_subset 
##        AIC      BIC    logLik
##   573.6497 599.4993 -276.8248
## 
## Random effects:
##  Formula: ~Surgery | Clinic
##  Structure: General positive-definite, Log-Cholesky parametrization
##             StdDev   Corr  
## (Intercept) 5.766445 (Intr)
## Surgery     2.514861 -0.772
## Residual    3.495341       
## 
## Fixed effects:  Post_QoL ~ Base_QoL + Age + Gender + BDI + Surgery 
##                 Value Std.Error DF   t-value p-value
## (Intercept) 28.394897  4.225895 84  6.719262  0.0000
## Base_QoL     0.147063  0.055898 84  2.630933  0.0101
## Age          0.198532  0.060153 84  3.300452  0.0014
## Gender      -4.696973  1.523294 84 -3.083433  0.0028
## BDI          0.472555  0.059681 84  7.918080  0.0000
## Surgery     -3.163418  1.240325 84 -2.550476  0.0126
##  Correlation: 
##          (Intr) Bas_QL Age    Gender BDI   
## Base_QoL -0.709                            
## Age      -0.020 -0.252                     
## Gender    0.051  0.048 -0.643              
## BDI      -0.258 -0.016 -0.527  0.381       
## Surgery  -0.248 -0.107 -0.060  0.096  0.112
## 
## Standardized Within-Group Residuals:
##        Min         Q1        Med         Q3        Max 
## -2.4115387 -0.5446037 -0.1070166  0.5074123  3.4520485 
## 
## Number of Observations: 98
## Number of Groups: 9

Bei körperlichen Gründen führt die OP zu einer leichten, jedoch nicht signifikanten Verbesserung der Lebensqualität im Vergleich zur Wartekontrollgruppe. Wird die OP jedoch nur zur Aussehensverbesserung durchgeführt, senkt sie die Lebensqualität sogar signifikant.

Aufgabe 6

Zunächst müssen die Daten für Random Effects eine hierarchische Struktur aufweisen, was für eine Interaktion nicht nötig ist. Der Hauptunterschied liegt jedoch darin, ob wir uns für den Einfluss der einen Variable auf die Wirkung der anderen interessieren, oder ob wir sie eher als störende Autokorrelationsquelle betrachten. So wollten wir in diesem Beispiel ja explizit wissen, welche Wirkung eine Schönheits-OP bei welchem dahinterstehenden Grund hat. Uns interessierte dagegen nicht, bei welcher Klinik welcher OP-Effekt vorlag. Der Random Effect diente hier also primär zur Reduktion von Fehlervarianz.
Wären wir jedoch auf der Suche nach dem besten Schönheitschirurgen der Stadt, etwa weil wir selbst eine solche Operation anstreben, dann wäre ein besseres Vorgehen, Clinic als Prädiktor ins Modell aufzunehmen, obwohl wir aufgrund der Dummykodierung hierdurch neun zusätzliche Prädiktoren, und durch die Interaktion mit Surgery nochmal neun weitere Prädiktoren hätten. Wir würden uns dann für diejenige Klinik entscheiden, bei der der Estimate von KlinikDummy:Surgery am höchsten ist, da die OPs an dieser Klinik offensichtlich die Lebensqualität am meisten steigern.
Grundsätzlich sollte man nur Random Effects von Variablen erlauben, bei denen es plausible Gründe für die Annahme gibt, dass sich ihr Effekt zwischen den verschiedenen Level-2-Variablen unterscheidet. Dass die Qualität der durchgeführten Operationen von Klinik zu Klinik schwanken, ist in unserem Beispiel deutlich plausibler, als dass der Effekt von Alter von Klinik zu Klinik schwankt. Wenn wir aber beispielsweise wüssten, dass sich die Kliniken hinsichtlich der psychologischen Betreuung unterscheiden, könnten wir beispielsweise auch für die Depressivität (BDI) einen Random Effect vermuten. Hier besteht also teils erheblicher Interpretationsspielraum, für den es nicht immer eine eindeutig beste Lösung gibt.

Aufgabe 7

cosmetic_surgery$Surgery_Text <- factor(cosmetic_surgery$Surgery_Text,
                                        levels = c("Waiting List",
                                                   "Cosmetic Surgery"))
cosmetic_surgery$Reason_Text <- factor(cosmetic_surgery$Reason_Text,
                                       levels = c("Change Appearance",
                                                  "Physical reason"))

plot_cosmetic_surgery <- ggplot(data = cosmetic_surgery,
                                aes(
                                  x = Surgery_Text,
                                  y = Post_QoL,
                                  group = Reason_Text,
                                  color = Reason_Text))

plot_cosmetic_surgery + stat_summary(fun.y = mean, geom = "line") +
  scale_x_discrete(rev(levels(cosmetic_surgery$Surgery_Text)))

## Warning: `fun.y` is deprecated. Use `fun` instead.

plot_cosmetic_surgery + 
  stat_summary(fun.data = mean_se, geom = "errorbar", width=0.2) +
  stat_summary(fun.y = mean, geom = "line") +
  scale_x_discrete(rev(levels(cosmetic_surgery$Surgery_Text)))

## Warning: `fun.y` is deprecated. Use `fun` instead.

# oder etwas schöner nicht überlappend
plot_cosmetic_surgery + 
  stat_summary(fun.data = mean_se, position=position_dodge(width=0.2), geom = "errorbar", width=0.2) +
  stat_summary(fun.y = mean, geom = "line", position=position_dodge(width=0.2)) +
  scale_x_discrete(rev(levels(cosmetic_surgery$Surgery_Text)))

## Warning: `fun.y` is deprecated. Use `fun` instead.

#zur Veranschaulichung der Random Effects:

plot_cosmetic_surgery + stat_summary(fun.y = mean, geom = "line", position=position_dodge(width=0.2)) +
  scale_x_discrete(rev(levels(cosmetic_surgery$Surgery_Text))) +
  stat_summary(fun.data = mean_se, position=position_dodge(width=0.2), geom = "errorbar", width=0.2) +
  facet_wrap(~Clinic, ncol = 3)

## Warning: `fun.y` is deprecated. Use `fun` instead.

Literatur

Anmerkung: Diese Übungszettel basieren zum Teil auf Aufgaben aus dem Lehrbuch Dicovering Statistics Using R (Field, Miles & Field, 2012). Sie wurden für den Zweck dieser Übung modifiziert, und der verwendete R-Code wurde aktualisiert.

Field, A., Miles, J., & Field, Z. (2012). Discovering Statistics Using R. London: SAGE Publications Ltd.

English

Links

Exercise sheet in PDF

Exercise sheet with solutions included

Exercise sheet with solutions included as PDF

The source code of this sheet as .Rmd (Right click and “store as” to download …)

Some hints

1. Please give your answers in a .Rmd file. You may generate one from scratch using the file menu: ‘File > new file > R Markdown …’ Delete the text below Setup Chunk (starting from line 11). Alternatively you may use this sample Rmd by downloading it.

2. You may find the informations useful that you can find on the front page of this course.

3. Don’t hesitate to google for solutions. Effective web searches to find solutions for R-problems is a very useful ability, professionals do that too … A really good starting point might be the R area of the programmers platform Stackoverflow

4. You can find very useful cheat sheets for various R-related topics. A good starting point is the Base R Cheat Sheet.

Ressources

Multi-Level-Models are explained in Field’s chapter 19, we highly recommend his accessible work!

general remarks

In Sheet 07, you’ve already encountered multi-level concepts in the context of repeated measures designs. There, we interpreted different time points on level 1 as nested within people (or countries) on level 2. Since different data points of one person are probobaly more similar to each other than to data points of a different person, the simple-regression assumption of independency is most likely violated (i.e., there’s autocorrelation). Multi-level-models allow you to include the dependency structure resulting from the hierarchical structure of your data in your analysis. In examples, apart from people and countries, frequent level 2 entities you’ll encounter are clinics or schools .

1 Data and packages

Load the tidyverse and the nlme-package. Read in the data at https://md.psych.bio.uni-goettingen.de/mv/data/div/cosmetic_surgery.dat (separated by tabs). This dataframe contains data on the quality of life of people who’ve decided to undergo cosmetic surgery:

Variable	explanation
particnu	participant number
Post_QoL	Quality of life after surgery/after the waiting time
Base_QoL	Baseline quality of life at the outset of the study
Clinic	Code number identifying in which of the ten clinics studied the surgery was performed
Surgery	1 = surgery performed, 0 = waiting list control
Reason	1 = there was a physicla reason for the surgery, 0 = the surgery was performed for purely aesthetical reason
Age	age of the subject
Gender	1 = male, 0 = female
BDI	standardized depression questionnaire
_Text-Variables	Repetition of surgery, reason and gender as strings, for use as plot labels

We’re mostly interested in whether quality of life really is higher after cosmetic surgery as compared to a waiting list control. Baseline quality of life, age, gender and depression were measured as covariates. Using clinic as a level 2 variable allows us, to control for possible differences in the quality of surgeries. The variable Reason will enable a more detailed analysis of surgery effects later on.

2 - Fixed Effects only

First, create a model without random effects, using only the covariates as predictors. Instead of lm(), use the nlme::gls() command - the syntax is the same, you’ll just have to add method = "ML". Creating a fixed-effects-model this way has the advantage of being able to compare this model with a random-effects-model using the anova() command. Now, create another model without random effects, using surgery as well as the covariates as predictors. If you want, try creating this new model from the old one using the update() command. Should that not work, you can just as well create it the old-fashioned way. Compare these two models with anova(). Additionally, check out the summary() of the second model, especially the surgery predictor.

3 - Random Effects

It’s possible that there’s considerable variability in the quality of the surgeries performed depending on the clinic, e.g. due to differing skills of the surgeons.
Create a model with the same predictors as the second model from the task before. For this new model however, allow the surgery effect to differ between the different clinics. Try figuring out the syntax for this yourself, using ?nlme::lme(). If that doesn’t work, you’ll find an explanation in “sheet_repeated_measure”. To make sure that random effects are warranted, there’s often an analysis of the 95% confidence intervals of the standard deviation of the effects. You get these using nlme::intervals(model_name). Perform this command. Do the results support allowing random effects of not? Why? Use another anova() to analyse whether the new model’s able to explain significantly more variance than the last one. Also check out the summar() again. What do you notice?

4 - Interaction 1

Up until now, we haven’t considered the reasons for getting cosmetic surgery. However, it’s definitely possible that cosmetic surgery that is supposed to help with a physical problem (e.g., a breast reduction to help with back pain) might have a different effect on quality of life than a surgery for purely aesthetic reasons. Expand your model from the last ask to include the interaction between surgery and reason. Again, perform an ANOVA and look at the summary.

5 - Interaction 2

The last task’s result indicates that a surgery’s effect depends on the reasons for the surgery. To investigate this further, create two new models with random effects, one for each reason. These models should not include the interaction term from before, because there’s always only one reason per model. First, try to solve this task using the subset()-option of nlme::lme(). If that doesn’t work, you can create two separate dataframes using dplyr::filter(). Look at both models and interpret them.

6 - Random Effects vs Interaction

Both with random effects and interactions, the effect of one variable varies depending on the value of another variable. In our example, surgery’s effect was affected both by different clinics and by different reasons. Can you explain the theoretical difference between these? What kind of question might have warranted looking at an interaction between clinic and surgery? And why did we allow surgery’s effect to vary betwen clincs, but not the effects of age, gender or depression?

7 - Plot

Create a fitting visualization for the conclusions of this analysis!

solution

task 1

library(tidyverse)
library(nlme)
cosmetic_surgery <- read_delim("https://md.psych.bio.uni-goettingen.de/mv/data/div/cosmetic_surgery.dat", "\t")

## Rows: 276 Columns: 12
## ── Column specification ──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
## Delimiter: "\t"
## chr (3): Surgery_Text, Reason_Text, Gender_Text
## dbl (9): particnu, Post_QoL, Base_QoL, Clinic, Surgery, Reason, Age, Gender, BDI
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

Aufgabe 2

m_covariates <- nlme::gls(Post_QoL ~ Base_QoL + Age + Gender + BDI,
                          data = cosmetic_surgery,
                          method = "ML")

m_fixed_effects <- update(m_covariates, .~. + Surgery)
# m_fixed_effects <- nlme::gls(Post_QoL ~ Base_QoL + Age + Gender + BDI + Surgery,
#       data = cosmetic_surgery,
#       method = "ML")

anova(m_covariates, m_fixed_effects)

##                 Model df      AIC      BIC    logLik   Test  L.Ratio p-value
## m_covariates        1  6 1811.000 1832.723 -899.5001                        
## m_fixed_effects     2  7 1807.024 1832.367 -896.5122 1 vs 2 5.975713  0.0145

summary(m_fixed_effects)

## Generalized least squares fit by maximum likelihood
##   Model: Post_QoL ~ Base_QoL + Age + Gender + BDI + Surgery 
##   Data: cosmetic_surgery 
##        AIC      BIC    logLik
##   1807.024 1832.367 -896.5122
## 
## Coefficients:
##                 Value Std.Error   t-value p-value
## (Intercept) 16.674777 2.6650283  6.256885  0.0000
## Base_QoL     0.507214 0.0456353 11.114509  0.0000
## Age          0.231408 0.0537485  4.305381  0.0000
## Gender      -1.016271 1.3183664 -0.770856  0.4415
## BDI          0.130171 0.0367314  3.543866  0.0005
## Surgery     -1.956844 0.8049688 -2.430957  0.0157
## 
##  Correlation: 
##          (Intr) Bas_QL Age    Gender BDI   
## Base_QoL -0.829                            
## Age      -0.170 -0.263                     
## Gender    0.116  0.011 -0.695              
## BDI       0.110 -0.160 -0.522  0.687       
## Surgery  -0.021 -0.022 -0.136 -0.079  0.082
## 
## Standardized residuals:
##         Min          Q1         Med          Q3         Max 
## -2.02107403 -0.75479722 -0.02969699  0.66683065  3.73213879 
## 
## Residual standard error: 6.229488 
## Degrees of freedom: 276 total; 270 residual

#second model significantly better
#surgery significant predictor

task 3

m_random <- nlme::lme(Post_QoL ~ Base_QoL + Age + Gender + BDI + Surgery,
                      data = cosmetic_surgery,
                      random = ~Surgery | Clinic,
                      method = "ML")

nlme::intervals(m_random)

## Approximate 95% confidence intervals
## 
##  Fixed effects:
##                   lower       est.      upper
## (Intercept) 27.33892415 34.1257908 40.9126574
## Base_QoL     0.10270956  0.2031965  0.3036834
## Age          0.18480652  0.2835955  0.3823844
## Gender      -2.98124310 -0.7258946  1.5294539
## BDI          0.03362661  0.1012212  0.1688157
## Surgery     -3.14643520 -0.9008826  1.3446700
## 
##  Random Effects:
##   Level: Clinic 
##                               lower       est.       upper
## sd((Intercept))           2.9765969  4.8911584  8.03717504
## sd(Surgery)               1.2816426  2.7914417  6.07981234
## cor((Intercept),Surgery) -0.9811474 -0.8268695 -0.02891568
## 
##  Within-group standard error:
##    lower     est.    upper 
## 4.911125 5.356450 5.842157

# neither the confidence interval of sd((Intercept)), nor the confidence 
# interval of sd(Surgery) include zero, so that we can assume a significant
# variability of effects. This justifies the use of random effects.


anova(m_fixed_effects, m_random)

##                 Model df      AIC      BIC    logLik   Test  L.Ratio p-value
## m_fixed_effects     1  7 1807.024 1832.367 -896.5122                        
## m_random            2 10 1763.483 1799.687 -871.7414 1 vs 2 49.54154  <.0001

summary(m_random)

## Linear mixed-effects model fit by maximum likelihood
##   Data: cosmetic_surgery 
##        AIC      BIC    logLik
##   1763.483 1799.687 -871.7414
## 
## Random effects:
##  Formula: ~Surgery | Clinic
##  Structure: General positive-definite, Log-Cholesky parametrization
##             StdDev   Corr  
## (Intercept) 4.891158 (Intr)
## Surgery     2.791442 -0.827
## Residual    5.356450       
## 
## Fixed effects:  Post_QoL ~ Base_QoL + Age + Gender + BDI + Surgery 
##                Value Std.Error  DF   t-value p-value
## (Intercept) 34.12579  3.484780 261  9.792811  0.0000
## Base_QoL     0.20320  0.051596 261  3.938226  0.0001
## Age          0.28360  0.050724 261  5.590941  0.0000
## Gender      -0.72589  1.158030 261 -0.626836  0.5313
## BDI          0.10122  0.034707 261  2.916443  0.0038
## Surgery     -0.90088  1.153000 261 -0.781338  0.4353
##  Correlation: 
##          (Intr) Bas_QL Age    Gender BDI   
## Base_QoL -0.782                            
## Age      -0.109 -0.259                     
## Gender    0.086 -0.017 -0.619              
## BDI      -0.050 -0.004 -0.485  0.631       
## Surgery  -0.231 -0.074 -0.094 -0.041  0.045
## 
## Standardized Within-Group Residuals:
##        Min         Q1        Med         Q3        Max 
## -2.0491309 -0.7468661 -0.1504611  0.6704928  2.9847839 
## 
## Number of Observations: 276
## Number of Groups: 10

#new model explains significantl more variance, but
#surgery's no longer a significant predictor.

task 4

m_interaction <- update(m_random, .~. + Surgery:Reason)
# m_interaction <- nlme::lme(Post_QoL ~ Base_QoL +
#                              Age +
#                              Gender +
#                              BDI +
#                              Surgery +
#                              Surgery:Reason,
#                            data = cosmetic_surgery,
#                            random = ~Surgery | Clinic,
#                            method = "ML")
anova(m_random, m_interaction)

##               Model df      AIC      BIC    logLik   Test  L.Ratio p-value
## m_random          1 10 1763.483 1799.687 -871.7414                        
## m_interaction     2 11 1748.482 1788.306 -863.2410 1 vs 2 17.00082  <.0001

summary(m_interaction)

## Linear mixed-effects model fit by maximum likelihood
##   Data: cosmetic_surgery 
##        AIC      BIC   logLik
##   1748.482 1788.306 -863.241
## 
## Random effects:
##  Formula: ~Surgery | Clinic
##  Structure: General positive-definite, Log-Cholesky parametrization
##             StdDev   Corr  
## (Intercept) 4.463083 (Intr)
## Surgery     1.229488 -0.812
## Residual    5.231446       
## 
## Fixed effects:  Post_QoL ~ Base_QoL + Age + Gender + BDI + Surgery + Surgery:Reason 
##                    Value Std.Error  DF   t-value p-value
## (Intercept)    31.324071  3.273622 260  9.568629  0.0000
## Base_QoL        0.222262  0.050330 260  4.416098  0.0000
## Age             0.294238  0.047655 260  6.174378  0.0000
## Gender         -1.248325  1.135014 260 -1.099832  0.2724
## BDI             0.155575  0.035984 260  4.323450  0.0000
## Surgery        -4.484673  1.152615 260 -3.890868  0.0001
## Surgery:Reason  5.819623  1.312833 260  4.432874  0.0000
##  Correlation: 
##                (Intr) Bas_QL Age    Gender BDI    Surgry
## Base_QoL       -0.798                                   
## Age            -0.033 -0.304                            
## Gender          0.114 -0.029 -0.639                     
## BDI            -0.078  0.017 -0.505  0.540              
## Surgery         0.007 -0.120 -0.034  0.053 -0.201       
## Surgery:Reason -0.114  0.067 -0.068 -0.124  0.353 -0.709
## 
## Standardized Within-Group Residuals:
##        Min         Q1        Med         Q3        Max 
## -2.2788267 -0.7482315 -0.1017572  0.6956910  3.0657904 
## 
## Number of Observations: 276
## Number of Groups: 10

#new model significantly better, the interaction reaches significance
#interestingly, surgery's significant again, too!

task 5

physical_subset <- cosmetic_surgery$Reason == 1
appearance_subset <- !physical_subset #! means "not"

m_physical <- update(m_random, subset = physical_subset)
m_appearance <- update(m_random, subset = appearance_subset)

# Alternative without subset
# df_physical <- dplyr::filter(cosmetic_surgery, Reason == 1)
# df_appearance <- dplyr::filter(cosmetic_surgery, Reason == 0)
# 
# m_physical_2 <- update(m_random, data = df_physical)
# m_appearance_2 <- update(m_random, data = df_appearance)

summary(m_physical)

## Linear mixed-effects model fit by maximum likelihood
##   Data: cosmetic_surgery 
##   Subset: physical_subset 
##        AIC      BIC    logLik
##   1154.884 1186.702 -567.4421
## 
## Random effects:
##  Formula: ~Surgery | Clinic
##  Structure: General positive-definite, Log-Cholesky parametrization
##             StdDev   Corr  
## (Intercept) 4.679271 (Intr)
## Surgery     1.668346 -0.794
## Residual    5.471991       
## 
## Fixed effects:  Post_QoL ~ Base_QoL + Age + Gender + BDI + Surgery 
##                 Value Std.Error  DF   t-value p-value
## (Intercept) 29.893046  4.411360 163  6.776379  0.0000
## Base_QoL     0.265651  0.069800 163  3.805898  0.0002
## Age          0.274836  0.066637 163  4.124374  0.0001
## Gender      -0.460955  1.502015 163 -0.306891  0.7593
## BDI          0.118640  0.064379 163  1.842838  0.0672
## Surgery      0.666083  1.162912 163  0.572772  0.5676
##  Correlation: 
##          (Intr) Bas_QL Age    Gender BDI   
## Base_QoL -0.834                            
## Age      -0.067 -0.283                     
## Gender    0.147 -0.079 -0.613              
## BDI      -0.184  0.087 -0.381  0.518       
## Surgery  -0.036 -0.055 -0.209 -0.111  0.005
## 
## Standardized Within-Group Residuals:
##        Min         Q1        Med         Q3        Max 
## -2.3720422 -0.7824301 -0.1276846  0.6989338  2.9496881 
## 
## Number of Observations: 178
## Number of Groups: 10

summary(m_appearance)

## Linear mixed-effects model fit by maximum likelihood
##   Data: cosmetic_surgery 
##   Subset: appearance_subset 
##        AIC      BIC    logLik
##   573.6497 599.4993 -276.8248
## 
## Random effects:
##  Formula: ~Surgery | Clinic
##  Structure: General positive-definite, Log-Cholesky parametrization
##             StdDev   Corr  
## (Intercept) 5.766445 (Intr)
## Surgery     2.514861 -0.772
## Residual    3.495341       
## 
## Fixed effects:  Post_QoL ~ Base_QoL + Age + Gender + BDI + Surgery 
##                 Value Std.Error DF   t-value p-value
## (Intercept) 28.394897  4.225895 84  6.719262  0.0000
## Base_QoL     0.147063  0.055898 84  2.630933  0.0101
## Age          0.198532  0.060153 84  3.300452  0.0014
## Gender      -4.696973  1.523294 84 -3.083433  0.0028
## BDI          0.472555  0.059681 84  7.918080  0.0000
## Surgery     -3.163418  1.240325 84 -2.550476  0.0126
##  Correlation: 
##          (Intr) Bas_QL Age    Gender BDI   
## Base_QoL -0.709                            
## Age      -0.020 -0.252                     
## Gender    0.051  0.048 -0.643              
## BDI      -0.258 -0.016 -0.527  0.381       
## Surgery  -0.248 -0.107 -0.060  0.096  0.112
## 
## Standardized Within-Group Residuals:
##        Min         Q1        Med         Q3        Max 
## -2.4115387 -0.5446037 -0.1070166  0.5074123  3.4520485 
## 
## Number of Observations: 98
## Number of Groups: 9

For physical reasons, surgery leads to a slight, non-significant increase in quality of life compared to the control group. A surgery for aesthetic reasons even significantly decreases quality of life!

task 6

First, for random effects to be possible there needs to be a hierarchical structure to the data. This isn’t neccessary for an interaction. However, the main difference is whether we’re theoretically interested in the effect of one variable on the effect of another, or whether we consider them an annoying source of autocorrelation. In our example, we were explicitly interested in the dependency of surgery’s effect on the reason for the surgery. In comparison, we were not interested in the question which clinic lead to which surgery effect. We included the random effects primarily to reduce the residual variance. However, had we been searching for the best cosmetic surgery clinic in the city, we would have included Clinic as predictor, even though that would’ve led to nine additional dummy predictors. After that, we would’ve chosen the clinic with the highest estimate for clinic_dummy:Surgery - that’s the clinic whose surgeries increase quality of life the most.

As a rule, you should only include random effects of variables where there are plausible reasons for the assumption that their effects differ between different level 2 entities. In our example, it was highly plausible that the quality of surgeries varied from clinic to clinic. In comparison, there was no reason to assume the age’s effect would differ depending on the clinic. However, if we had known that the clinic differ regarding their psychological counselling, we might have included random effects for depression as well. There is some room for interpretation when deciding which random effects to include, no one solution is neccessarily the “correct” one.

task 7

cosmetic_surgery$Surgery_Text <- factor(cosmetic_surgery$Surgery_Text,
                                        levels = c("Waiting List",
                                                   "Cosmetic Surgery"))
cosmetic_surgery$Reason_Text <- factor(cosmetic_surgery$Reason_Text,
                                       levels = c("Change Appearance",
                                                  "Physical reason"))

plot_cosmetic_surgery <- ggplot(data = cosmetic_surgery,
                                aes(
                                  x = Surgery_Text,
                                  y = Post_QoL,
                                  group = Reason_Text,
                                  color = Reason_Text))

plot_cosmetic_surgery + stat_summary(fun.y = mean, geom = "line") +
  scale_x_discrete(rev(levels(cosmetic_surgery$Surgery_Text)))

## Warning: `fun.y` is deprecated. Use `fun` instead.

plot_cosmetic_surgery + 
  stat_summary(fun.data = mean_se, geom = "errorbar", width=0.2) +
  stat_summary(fun.y = mean, geom = "line") +
  scale_x_discrete(rev(levels(cosmetic_surgery$Surgery_Text)))

## Warning: `fun.y` is deprecated. Use `fun` instead.

# or a little prettier without overlaying bars:
plot_cosmetic_surgery + 
  stat_summary(fun.data = mean_se, position=position_dodge(width=0.2), geom = "errorbar", width=0.2) +
  stat_summary(fun.y = mean, geom = "line", position=position_dodge(width=0.2)) +
  scale_x_discrete(rev(levels(cosmetic_surgery$Surgery_Text)))

## Warning: `fun.y` is deprecated. Use `fun` instead.

#visualizing the random effects

plot_cosmetic_surgery + stat_summary(fun.y = mean, geom = "line", position=position_dodge(width=0.2)) +
  scale_x_discrete(rev(levels(cosmetic_surgery$Surgery_Text))) +
  stat_summary(fun.data = mean_se, position=position_dodge(width=0.2), geom = "errorbar", width=0.2) +
  facet_wrap(~Clinic, ncol = 3)

## Warning: `fun.y` is deprecated. Use `fun` instead.

Literature

Annotation: This exercise sheet bases in part on exercises, that you can find in the textbook Dicovering Statistics Using R (Field, Miles & Field, 2012). They were modified for the purpose of this sheet and the R-code was actualized.

Field, A., Miles, J., & Field, Z. (2012). Discovering Statistics Using R. London: SAGE Publications Ltd.