---
title: 'Übungszettel Multi-Level-Modelle'
author: "Joshua Driesen (joshua.driesen@stud.uni-goettingen.de)"
subtitle: "M.Psy.205, Dozent: Peter Zezula"
output:
  html_document:
    theme: paper
    toc: yes
    toc_depth: 3
    toc_float: yes
  pdf_document:
    latex_engine: xelatex
header-includes:
 - \usepackage{comment}
params:
 soln: TRUE
mainfont: Helvetica
linkcolor: blue
---

<!-- Hier und im Titel immer auf die aktuelle Version des Übungszettels aktualisieren -->
# German

## Links

[Übungszettel als PDF-Datei zum Drucken](https://md.psych.bio.uni-goettingen.de/mv/sheets/sheet_mlm.pdf)

`r if(!params$soln) {"\\begin{comment}"}`
**Übungszettel mit Lösungen**

[Lösungszettel als PDF-Datei zum Drucken](https://md.psych.bio.uni-goettingen.de/mv/sheets/sheet_mlm_solutions.pdf)

[Der gesamte Übugszettel als .Rmd-Datei](https://md.psych.bio.uni-goettingen.de/mv/sheets/sheet_mlm.Rmd) (Zum Downloaden: Rechtsklick > Speichern unter...)

`r if(!params$soln) {"\\end{comment}"}`

## 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](https://md.psych.bio.uni-goettingen.de/mv/sheets/students_template.Rmd) 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](https://www.psych.uni-goettingen.de/de/it/team/zezula/courses/multivariate)
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](https://stackoverflow.com/questions/tagged/r)
4. Auf der Website von R Studio finden Sie sehr [hilfreiche Übersichtszettel](https://www.rstudio.com/resources/cheatsheets/) zu vielen verschiedenen R-bezogenen Themen. Ein guter Anfang ist der [Base R Cheat Sheet](http://github.com/rstudio/cheatsheets/raw/master/base-r.pdf)

## 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!

`r if(!params$soln) {"\\begin{comment}"}`
## Lösung
### Aufgabe 1

```{r}
library(tidyverse)
library(nlme)
cosmetic_surgery <- read_delim("https://md.psych.bio.uni-goettingen.de/mv/data/div/cosmetic_surgery.dat", "\t")
# alternativ
cosmetic_surgery <- readr::read_tsv("https://md.psych.bio.uni-goettingen.de/mv/data/div/cosmetic_surgery.dat")

```

### Aufgabe 2

```{r}
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)
summary(m_fixed_effects)
#zweites Modell signifikant bessere Varianzaufklärung,
#Surgery signifikanter Prädiktor
```

### Aufgabe 3

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

nlme::intervals(m_random)
# 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)
summary(m_random)
#Das neue Modell klärt signifikant mehr Varianz auf, jedoch ist
#der Fixed Effect von Surgery nun nicht mehr signifikant!
```

### Aufgabe 4

```{r}
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)
summary(m_interaction)
#Neues Modell signifikant besser, Interaktion signifikanter Prädiktor.
#Interessanterweise wird auch Surgery selbst wieder signifikant!
```

### Aufgabe 5

```{r}
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)
summary(m_appearance)
```

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

```{r}
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)))


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)))

# 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)))
 

#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)

```

`r if(!params$soln) {"\\end{comment}"}`

## 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](https://md.psych.bio.uni-goettingen.de/mv/sheets/sheet_mlm.pdf)

`r if(!params$soln) {"\\begin{comment}"}`
**Exercise sheet with solutions included**

[Exercise sheet with solutions included as PDF](https://md.psych.bio.uni-goettingen.de/mv/sheets/sheet_mlm.pdf)

[The source code of this sheet as .Rmd](https://md.psych.bio.uni-goettingen.de/mv/sheets/sheet_mlm.Rmd) 
(Right click and "store as" to download ...)


`r if(!params$soln) {"\\end{comment}"}`


## 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](https://md.psych.bio.uni-goettingen.de/mv/sheets/students_template.Rmd) by downloading it.

2. You may find the informations useful that you can find on the
   [front page of this course](https://www.psych.uni-goettingen.de/de/it/team/zezula/courses/multivariate).
   
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](https://stackoverflow.com/questions/tagged/r)
4. You can find very useful [cheat sheets](https://www.rstudio.com/resources/cheatsheets/) 
   for various R-related topics.
   A good starting point is the [Base R Cheat Sheet](http://github.com/rstudio/cheatsheets/raw/master/base-r.pdf).


## 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!

`r if(!params$soln) {"\\begin{comment}"}`
## solution
### task 1

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

```

### Aufgabe 2

```{r}
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)
summary(m_fixed_effects)
#second model significantly better
#surgery significant predictor
```

### task 3

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

nlme::intervals(m_random)
# 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)
summary(m_random)
#new model explains significantl more variance, but
#surgery's no longer a significant predictor.
```

### task 4

```{r}
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)
summary(m_interaction)
#new model significantly better, the interaction reaches significance
#interestingly, surgery's significant again, too!
```

### task 5

```{r}
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)
summary(m_appearance)
```

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

```{r}
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)))


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)))

# 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)))
 

#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)

```

`r if(!params$soln) {"\\end{comment}"}`

## 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.
