--- title: 'Übungszettel MANOVA' author: "Johannes Brachem (johannes.brachem@stud.uni-goettingen.de)" subtitle: "M.Psy.205, Dozent: Dr. 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 --- # Deutsche Version ## Links [Übungszettel als PDF-Datei zum Drucken](https://pzezula.pages.gwdg.de/sheet_manova.pdf) `r if(!params$soln) {"\\begin{comment}"}` **Übungszettel mit Lösungen** [Lösungszettel als PDF-Datei zum Drucken](https://pzezula.pages.gwdg.de/sheet_manova_solutions.pdf) [Der gesamte Übugszettel als .Rmd-Datei](https://pzezula.pages.gwdg.de/sheet_manova.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://pzezula.pages.gwdg.de/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 Da es sich um eine praktische Übung handelt, können wir Ihnen nicht alle neuen Befehle einzeln vorstellen. Stattdessen finden Sie hier Verweise auf sinnvolle Ressourcen, in denen Sie für die Bearbeitung unserer Aufgaben nachschlagen können. Ressource | Beschreibung -------------|-------------------------------- Field, Kapitel 16 | Buchkapitel, das Schritt für Schritt erklärt, worum es geht, und wie man Multivariate Varianzanalysen in R durchführt. **Große Empfehlung!** ## Tipp der Woche Hinweis: Der Tipp ist diese Woche etwas länger. ### 1) Befehle nutzen, ohne Pakete zu laden Normalerweise müssen Sie ein Paket mit `library()` laden, damit Sie Funktionen aus diesem Paket nutzen können. Wenn Sie eine Funktion aber nur sehr selten benötigen, können Sie einen Trick verwenden, nämlich den doppelten Doppelpunkt `::`. Das funktioniert nach dem Schema `package::function()` und wird von uns auch benutzt, um Funktionen mit häufigen Namen richtig einzusetzen. Wenn Sie diese Schreibweise verwenden, muss das jeweilige Paket nicht geladen sein. Beispiel: `psych::describe()`. Das Paket muss allerdings installiert sein, damit diese Schreibweise funktioniert. ### 2) Wie genau funktioniert eigentlich die Pipe? `%>%` Die Pipe kommt aus dem Paket `magrittr` und wird im `tidyverse` häufig verwendet. Sie kennen die Pipe vor allem aus `dplyr`, aber wenn Sie `magrittr` oder `dplyr` geladen haben (d.h. mit `library()` aktiviert), dann können Sie die Pipe für fast alle Befehle in R benutzen, wenn Sie möchten. Das funktioniert so: Der Code, der auf der *linken* Seite der Pipe steht, wird von R "unter der Haube" auf der *rechten* Seite der Pipe eingesetzt, und zwar standardmäßig immer als erstes *Argument* einer *Funktion*, so dass andere Argumente der Funktion mit einem Komma dahinter gestellt werden. Wenn Sie den Code auf der linken Seite der Pipe an eine andere Stelle auf der rechten Seite der Funktion einsetzen möchten, dann können Sie R mit einem Punkt `.` sagen, wo der Code der linken Seite eingesetzt werden soll (siehe Beispiel 4). #### Beispiel 1 Nehmen wir mal an, wir wollen nur VP über 18 auswählen. ```{r, eval = FALSE} # ohne Pipe filter(example_data, age > 18) # mit Pipe example_data %>% filter(age > 18) ``` #### Beispiel 2 Das ist besonders nützlich, wenn Befehle verschachtelt sind. Nehmen wir mal an, wir wollen nur VP über 18 auswählen und uns nur die Gruppe und unsere abhängige Variable anzeigen lassen. ```{r, eval = FALSE} # ohne Pipe select(filter(example_data, age > 18), group, dependent_variable) # mit Pipe example_data %>% filter(age > 18) %>% select(group, dependent_variable) ``` #### Beispiel 3 Das funktioniert für fast alle Funktionen. Hinweis: `na.rm = TRUE` sorgt dafür, dass Missing Values bei der Berechnung des Mittelwerts weggelassen werden. ```{r, eval = FALSE} # ohne Pipe mean(age, na.rm = TRUE) # mit Pipe age %>% mean(na.rm = TRUE) ``` #### Beispiel 4 ```{r, eval = FALSE} # ohne Pipe lm(dependent_variable ~ group, data = example_data) # mit Pipe example_data %>% lm(dependent_variable ~ group, data = .) ``` ## Vorgehen Hinweis: Screenshot aus Field, Miles & Field (2012). Für fast alle Verfahren finden Sie in diesem Lehrbuch ähnliche Übersichten. Diese sind sehr wertvoll für die Prüfungsvorbereitung und darüber hinaus. ![](data/screenshot_manova.png) ## 1) Daten einlesen 1. Laden Sie die nötigen Pakete und setzen Sie ein sinnvolles Arbeitsverzeichnis. 2. Laden Sie den Datensatz `ocd_data.dat` über den Link [https://pzezula.pages.gwdg.de/data/ocd_data.dat](https://pzezula.pages.gwdg.de/data/ocd_data.dat) in R ein. 3. Kodieren Sie die Variable `group` in `ocd_data` als Faktor mit sinnvoller Baseline. Geben Sie der Gruppe "No Treatment Control" das Label "NT". `r if(!params$soln) {"\\begin{comment}"}` ### Lösung #### Unteraufgabe 1 ```{r} library(tidyverse) ``` ```{r, eval = FALSE} # nur als Beispiel ... setwd("~/ownCloud/_Arbeit/Hiwi Peter/gitlab_sheets") ``` #### Unteraufgabe 2 ```{r} ocd_data <- read_delim("https://pzezula.pages.gwdg.de/data/ocd_data.dat", delim = "\t") ``` #### Unteraufgabe 3 ```{r} ocd_data$group <- ocd_data$group %>% factor(levels = c("No Treatment Control", "BT", "CBT"), labels = c("NT", "BT", "CBT")) ``` `r if(!params$soln) {"\\end{comment}"}` ## 2) Überblick über die Daten ### Bedeutung der Variablen Unser Beispieldatensatz enthält hypothetische Daten zur Evaluation von Therapieprogrammen bei Zwangsstörungen. Variable | Beschreibung -------------|-------------------------------- group | Faktor, der angibt, welche Art von Therapie die Versuchsperson erhalten hat. NT = Keine Therapie, BT = Verhaltenstherapie, CBT = Kognitive Verhaltenstherapie actions | Häufigkeit von Zwangshandlungen nach der Behandlung thoughts | Häufigkeit von Zwangsgedanken nach der Behandlung 1. Erstellen Sie zunächst einen einfachen Scatterplot, der den Zusammenhang zwischen zwanghaftem Verhalten auf der x-Achse und zwanghaften Gedanken auf der y-Achse darstellt. 2. Fügen Sie dem Plot eine Regressionslinie hinzu. 3. Nutzen Sie den Befehl `facet_wrap()`, um den Plot aus Aufgabe 2.2 nach Gruppen getrennt darzustellen. 4. Nutzen Sie den Code `ocd_data %>% select(actions, thoughts) %>% by(ocd_data$group, cov)`, um sich die Varianz-Kovarianz-Matrizen für jede Gruppe anzeigen zu lassen. 5. Nutzen Sie den Code `ocd_data %>% by(ocd_data$group, psych::describe)`, um sich für jede Gruppe deskriptive Daten ausgeben zu lassen. Können Sie den Befehl verstehen? *Falls Sie eine Fehlermeldung bekommen, installieren Sie das Paket `psych`.* 6. Führen Sie den unten stehenden Code aus, um die Annahme der multivariaten Normalverteilung an Ihren Daten zu überprüfen. Können Sie den Code verstehen? Sind die Daten in jeder Gruppe multivariat normalverteilt? *Falls Sie eine Fehlermeldung bekommen, installieren Sie das Paket `mvnormtest`.* 7. Der Box's M-Test prüft die Gleichheit der Varianz-Covarianz-Matrizen. Führen Sie den Befehl `boxM()` aus der `library(heplots)` durch, um diese Voraussetzung zu überprüfen. Können wir von Gleichheit der Varianz-Covarianz-Matrizen ausgehen? **Falls Sie eine Fehlermeldung bekommen, installieren Sie das Paket `heplots`.** ```{r, eval = FALSE} # Daten vorbereiten nt <- ocd_data %>% filter(group == "NT") %>% select(2:3) %>% t() bt <- ocd_data %>% filter(group == "BT") %>% select(2:3) %>% t() cbt <- ocd_data %>% filter(group == "CBT") %>% select(2:3) %>% t() # Tests durchführen mvnormtest::mshapiro.test(nt) mvnormtest::mshapiro.test(bt) mvnormtest::mshapiro.test(cbt) ``` `r if(!params$soln) {"\\begin{comment}"}` ### Lösung #### Unteraufgabe 1 ```{r} # Objekt erstellen baseplot <- ggplot(ocd_data, aes(x = actions, y = thoughts)) # Punkte hinzufügen scatterplot <- baseplot + geom_point() # Plot anzeigen scatterplot ``` #### Unteraufgabe 2 ```{r} # Regressionslinie hinzufügen lineplot <- scatterplot + geom_smooth(method = "lm") # Plot anzeigen lineplot ``` #### Unteraufgabe 3 ```{r} # Plot aufteilen groupplot <- lineplot + facet_wrap(~ group) # Plot anzeigen groupplot ``` #### Unteraufgabe 4 ```{r} ocd_data %>% select(actions, thoughts) %>% by(ocd_data$group, cov) ``` Hier können wir die Covarianzen in den verschiedenen Untergruppen direkt descriptiv vergleichen. #### Unteraufgabe 5 ```{r} ocd_data %>% by(ocd_data$group, psych::describe) ``` #### Unteraufgabe 6 Mit `t()` transponieren wir hier die Matrizen der Daten pro Gruppe. Vorher waren die Informationen von oben nach unten gespeichert (in Spalten), jetzt sind sie von links nach rechts gespeichert (in Zeilen). Das ist ungewöhnlich, aber nötig für den Test. ```{r} # Daten vorbereiten nt <- ocd_data %>% filter(group == "NT") %>% select(2:3) %>% t() bt <- ocd_data %>% filter(group == "BT") %>% select(2:3) %>% t() cbt <- ocd_data %>% filter(group == "CBT") %>% select(2:3) %>% t() # Tests durchführen mvnormtest::mshapiro.test(nt) mvnormtest::mshapiro.test(bt) mvnormtest::mshapiro.test(cbt) ``` Für die erste Gruppe haben wir ein signifikantes Ergebnis, d.h. dort sind die Daten nicht multivariat normalverteilt. Für den Zweck dieses Übungszettels führen wir die Analyse trotzdem fort. #### Unteraufgabe 7 ```{r} res.boxm <- heplots::boxM(ocd_data[,c('actions', 'thoughts')], group=ocd_data$group) res.boxm # summary(res.boxm) # for details ``` Der p-Wert des BoxM-Test ist nicht unter 0.05, also halten wir noch an der H0 fest, dass sich die Varianz-Covarianz-Matrizen nicht signifikant unterscheiden. `r if(!params$soln) {"\\end{comment}"}` ## 3) MANOVA durchführen ### Vorab: Erklärung MANOVAS können in R mit dem Befehl `manova()` durchgeführt werden. Dieser Befehl funktioniert genau so wie `lm()` und `aov()` in der Form: `manova(outcome ~ predictor, data = data)`. Der Unterschied ist, dass im Vorfeld alle verwendeten Outcome-Variablen mit dem `cbind()`-Befehl zu einem Objekt "zusammengeschnürt" werden. 1. Setzen Sie Kontraste für die Interpretation der Analyse. Das funktioniert genau so wie bei ANOVAs und Regressionen. Schauen Sie zur Not in Ihren Aufzeichnung von früheren Übungszetteln nach. a) Erster Kontrast: Vergleich von BT und NT b) Zweiter Kontrast: Vergleich von CBT und NT *Hinweis: Hier handelt es sich um nicht-orthogonale Kontraste. Das ist an dieser Stelle in Ordnung, weil wir nur eine Prädiktorvariable haben. (siehe Field, Kap. 16.6.6: Setting Contrasts)* *Bemerkung: Wir müssen keine Kontraste setzen. Bei den Default-Einstellungen würde `BT` zur Referenzgruppe und wir würden den Unterschied, auch der Gruppe `CBT` als Effekt prüfen.* 2. Erstellen Sie das Outcome-Objekt, indem Sie `ocd_data$thoughts` und `ocd_data$actions` mit Hilfe von `cbind()` verbinden. 3. Nutzen Sie den Befehl `manova()`, um die Analyse durchzuführen. Speichern Sie das Ergebnis in einem Objekt. 4. Wenden Sie auf das Objekt aus 3.3 den Befehl `summary()` mit dem zusätzlichen Argument `intercept = TRUE` an. 5. Welchen Schluss können Sie aus dem Output ziehen? `r if(!params$soln) {"\\begin{comment}"}` ### Lösung #### Unteraufgabe 1 ```{r} .bt.vs.nt <- c(-1,1,0) .cbt.vs.nt <- c(-1,0,1) # Kontraste an den Faktor binden contrasts(ocd_data$group) <- cbind(.bt.vs.nt, .cbt.vs.nt) mean(ocd_data$actions) ocd_data %>% group_by(group) %>% summarize(mean(actions)) ``` #### Unteraufgabe 2 ```{r} outcome <- cbind(ocd_data$actions, ocd_data$thoughts) ``` #### Unteraufgabe 3 ```{r} model1 <- manova(outcome ~ group, data = ocd_data) ``` #### Unteraufgabe 4 ```{r} summary(model1, intercept = TRUE) ``` Wir fordern den Intercept mit an, um zu zeigen, dass auch die Manova im Prinzp als Teil des allgemeinen linearen Modells verstanden werden kann. In den Linearkombinationen der Erklärvariablen gibt es ja immer auch einen Koeffizienten für den Intercept. #### Unteraufgabe 5 Der F-Test für die Gruppen wird signifikant (F(4,54) = 2.56, p = .049). Das bedeutet, dass die Art der Therapie die Zwanghaftigkeit der untersuchten Personen, gemessen durch Handlungen und Gedanken, beeinflusst hat. Mehr Schlüsse können wir erst einmal nicht ziehen. `r if(!params$soln) {"\\end{comment}"}` ## 4) MANOVA interpretieren 1. Wenden Sie den Befehl `summary.aov()` auf ihr MANOVA-Modell an. 2. Interpretieren Sie die Ergebnisse kurz. a) Dürfen wir auf Grundlage dieser Ergebnisse die geplanten Kontraste untersuchen? 3. Erstellen Sie für jedes einzelne Outcome ein eigenes ANOVA-Modell und betrachten Sie den Output, um Ihre Kontraste zu interpretieren. 4. Was folgern Sie aus den Ergebnissen? `r if(!params$soln) {"\\begin{comment}"}` ### Lösung #### Unteraufgabe 1 ```{r} summary.aov(model1) ``` #### Unteraufgabe 2 Beide einzelnen ANOVAs liefern keine signifikanten Ergebnisse, d.h. wir haben keinen Hinweis darauf, dass besonders die Zwangsgedanken oder besonders die Zwangshandlungen von der Art der Therapie beeinflusst werden. Das ist sehr interessant, da die Therapie die Kombination beider Komponenten ja scheinbar durchaus beeinflusst. Auf dieser Grundlage haben wir eigentlich keinen Anlass, die geplanten Kontraste zu untersuchen. Wir tun es dennoch, zu Demonstrationszwecken. #### Unteraufgabe 3 ```{r} out1 <- lm(actions ~ group, data = ocd_data) out2 <- lm(thoughts ~ group, data = ocd_data) summary(out1) summary(out2) ``` #### Unteraufgabe 4 Für Zwangshandlungen erweist sich die Verhaltenstherapie als signifikant besser als die Kontrollgruppe, für Zwangsgedanken dagegen die kognitive Verhaltenstherapie. `r if(!params$soln) {"\\end{comment}"}` ## 5) Rendern (knit) Lassen Sie die Datei mit `Strg` + `Shift` + `K` (Windows) oder `Cmd` + `Shift` + `K` (Mac) rendern. Sie sollten nun im "Viewer" unten rechts eine "schön aufpolierte" Version ihrer Datei sehen. Falls das klappt: Herzlichen Glückwunsch! Ihr Code kann vollständig ohne Fehlermeldung gerendert werden. Falls nicht: Nur mut, das wird schon noch! Gehen Sie auf Fehlersuche! Ansonsten schaffen wir es ja in der Übung vielleicht gemeinsam. ## 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 Version ## Links [Exercise sheet as PDF](https://pzezula.pages.gwdg.de/sheet_manova.pdf) `r if(!params$soln) {"\\begin{comment}"}` **Exercise sheet with solutions included** [Exercise sheet with solutions included as PDF](https://pzezula.pages.gwdg.de/sheet_manova_solutions.pdf) [The source code of this sheet as .Rmd](https://pzezula.pages.gwdg.de/sheet_manova.Rmd) (Right click and "save as" to download ...) `r if(!params$soln) {"\\end{comment}"}` ## Some hints 1. Please try to solve this sheet in an .Rmd file. You can create one from scratch using `File > New file > R Markdown...`. You can delete the text beneath *Setup Chunk* (starting from line 11). Alternatively, you can download our template file unter [this link](https://pzezula.pages.gwdg.de/students_template.Rmd) (right click > save as...). 2. You'll find a lot of the important information on the [website of this course](https://www.psych.uni-goettingen.de/de/it/team/zezula/courses/multivariate) 3. Please don't hesitate to search the web for help with this sheet. In fact, being able to effectively search the web for problem solutions is a very useful skill, even R pros work this way all the time! The best starting point for this is the [R section on the programming site Stackoverflow](https://stackoverflow.com/questions/tagged/r) 4. On the R Studio website, you'll find highly helpful [cheat sheets](https://www.rstudio.com/resources/cheatsheets/) for many of R topics. The [base R cheat sheet](http://github.com/rstudio/cheatsheets/raw/master/base-r.pdf) might be a good starting point. ## Ressources Since this is a hands-on seminar, we won't be able to present each and every new command to you explicitly. Instead, you'll find here references to helpful ressources that you can use for completing this sheets. Ressource | Description -------------|-------------------------------- Field, chapter 16 | Book chapter explaining step by step the why and how of multivariate analyses of variance in R. **Highly recommended!** ## Hint of the week Note: This week's hint will be a bit longer than usual. ### 1) Using commands without loading packages Usually, you'll have to load a package using `library()` in order to use functions from that package. However, there's a trick for when you'll only need a function sporadically: Adding the package name in front of the function using a double colon like this: `package::function()`. This also helps in using a function with a name that occurs in multiple packages. E.g., you might have noticed the warning `dplyr::filter() masks stats::filter()` when loading the tidyverse. Here, `filter()` alone will use the dplyr-command. If you want to use the filter-function from the stats-commands, you can do so by typing `stats::filter()`. The relevant packages don't need to be loaded for this to work. However, they have to be installed on your computer. ### 2) How exactly does the pipe work? `%>%` Originally, the pipe stems from the package `magrittr` and is used frequently in the `tidyverse`. You mostly know piping from `dplyr`, but as long as you have loaded either `magrittr` or `dplyr` (using the `library()`-function), you'll be able to use the pipe for almost all R commands - if that is something you want to do. It works like this: R inserts the code on the *left hand side* of the pipe "under the hood" into the *right hand side* of the pipe, the default being using it as the first argument of a function. This way you can add additional arguments with a comma behind it. If you want the left code to be inserted at a different location, you can use the dot `.` to tell R where to put it (see example no. 4). #### Example no. 1 Say you only want to include subjects older than 18. ```{r, eval = FALSE} # without pipe filter(example_data, age > 18) # with pipe example_data %>% filter(age > 18) ``` #### Example no. 2 This is most useful when you have nested commands. Here, we only want to include subjects older than 18, and only look at the group vairable and our dependend variable. ```{r, eval = FALSE} # without pipe select(filter(example_data, age > 18), group, dependent_variable) # with pipe example_data %>% filter(age > 18) %>% select(group, dependent_variable) ``` #### Example no. 3 This works for almost all functions. Note: `na.rm = TRUE` excludes missing values when calculating the mean. ```{r, eval = FALSE} # without pipe mean(age, na.rm = TRUE) # with pipe age %>% mean(na.rm = TRUE) ``` #### Example no. 4 ```{r, eval = FALSE} # without pipe lm(dependent_variable ~ group, data = example_data) # with pipe example_data %>% lm(dependent_variable ~ group, data = .) ``` ## Procedure Note: Screenshot from Field, Miles & Field (2012). There are summarys like this for almost all statistical procedures covered in this seminar. They're invaluable for preparing for our exam and for life beyond this seminar. ![](data/screenshot_manova.png) ## 1) Load data 1. Load the necessary packages and set a reasonable working directory. 2. Load the dataframe `ocd_data.dat` from the link [https://pzezula.pages.gwdg.de/data/ocd_data.dat](https://pzezula.pages.gwdg.de/data/ocd_data.dat) into R 3. Code the `group` variable as a factor with a reasonable baseline. Give the "No Treatment Control" group the lavel "NT". `r if(!params$soln) {"\\begin{comment}"}` ### Solution #### Subtask 1 ```{r} library(tidyverse) ``` ```{r, eval = FALSE} setwd("~/ownCloud/_Arbeit/Hiwi Peter/gitlab_sheets") ``` #### Subtask 2 ```{r} ocd_data <- read_delim("https://pzezula.pages.gwdg.de/data/ocd_data.dat", delim = "\t") ``` #### Subtask 3 ```{r} ocd_data$group <- ocd_data$group %>% factor(levels = c("No Treatment Control", "BT", "CBT"), labels = c("NT", "BT", "CBT")) ``` `r if(!params$soln) {"\\end{comment}"}` ## 2) Data overview ### Meaning of the variables Our data example contains hypothetical data of an evaluation study on different therapies for obsessive compulsive disorder Variable | Meaning -------------|-------------------------------- group | Factor specifying the therapy method: NT = No Treatment, BT = behavioral therapy, CBT = cognitive behavioral therapy actions | Frequency of obsessive actions after the therapy thoughts | Frequency of obsessive thoughts after the therapy 1. First, create a simple scatterplot showing the relationship between obsessive actions on the x-axis and obsessive thoughts on the y-axis. 2. Add a regression line to the plot. 3. Use the command `facet_wrap()` to show the plot from 2.2 separatly for each group. 4. Use the command `ocd_data %>% select(actions, thoughts) %>% by(ocd_data$group, cov)` to get separate variance-covariance-matrices for each group. 5. Use the command `ocd_data %>% by(ocd_data$group, psych::describe)` to get descriptive statistics for each group. Are you able to understand this command? **If you get an error message, you might have to install the package `psych`.** 6. Use the code below to investigate whether the data satisfy the assumption of being normally distributed. Can you understand the code? Is the data multivariately normally distributed for each group? **If you get an error message, you might have to install the package `mvnormtest`.** 7. The Box's M-test checks for equal variance-covariance-matrices. Use the command `boxM()` of `library(heplots)` to check that. Do we violate this assumtion? **If you get an error message, you might have to install the package `heplots`.** ```{r, eval = FALSE} # prepare data nt <- ocd_data %>% filter(group == "NT") %>% select(2:3) %>% t() bt <- ocd_data %>% filter(group == "BT") %>% select(2:3) %>% t() cbt <- ocd_data %>% filter(group == "CBT") %>% select(2:3) %>% t() # run the tests mvnormtest::mshapiro.test(nt) mvnormtest::mshapiro.test(bt) mvnormtest::mshapiro.test(cbt) ``` `r if(!params$soln) {"\\begin{comment}"}` ### Solution #### Subtask 1 ```{r} # Create plot object baseplot <- ggplot(ocd_data, aes(x = actions, y = thoughts)) # Add point geom scatterplot <- baseplot + geom_point() # Display plot scatterplot ``` #### Subtask 2 ```{r} # Add regression line lineplot <- scatterplot + geom_smooth(method = "lm") # Display plot lineplot ``` #### Subtask 3 ```{r} # Split plot groupplot <- lineplot + facet_wrap(~ group) # Display plot groupplot ``` #### Subtask 4 ```{r} ocd_data %>% select(actions, thoughts) %>% by(ocd_data$group, cov) ``` #### Subtask 5 ```{r} ocd_data %>% by(ocd_data$group, psych::describe) ``` #### Subtask 6 The `t()` command transposes the matrices containing the group data. Before, the information was coded from top to bottom (i.e., in columns). Now, the information's saved left to right (i.e., in rows). This is pretty unusal, but neccessary for this test. ```{r} # Prepare data nt <- ocd_data %>% filter(group == "NT") %>% select(2:3) %>% t() bt <- ocd_data %>% filter(group == "BT") %>% select(2:3) %>% t() cbt <- ocd_data %>% filter(group == "CBT") %>% select(2:3) %>% t() # Run tests mvnormtest::mshapiro.test(nt) mvnormtest::mshapiro.test(bt) mvnormtest::mshapiro.test(cbt) ``` We get a significant result for the first group, indicating that the data in that group is not multivariately normally distributed. Inspite of this, we continue the analysis for the sake of this sheet. `r if(!params$soln) {"\\end{comment}"}` #### Subtask 7 ```{r} res.boxm <- heplots::boxM(ocd_data[,c('actions', 'thoughts')], group=ocd_data$group) res.boxm # summary(res.boxm) # for details ``` The p-value of our BoxM-test is not below 0.05, therefore we stay with H0 and state, that the variance-covariance-matrices don't differ significantly. ## 3) Run the MANOVA ### First, some explanation You can run a MANOVA in R using the command `manova()`. This works exactly like the `lm()` and `aov()` commands you already know, the pattern being `manova(outcome ~ predictor, data = data)`. The difference, however, is that you have to bind together all relevant outcome variables beforehand, using the `cbind()` command. 1. Set the contrasts for this analysis. This works the same way as for ANOVAs and regressions. If you're having trouble, consult your notes on earlier sheets covering these topics. a) First contrast: comparison between BT and NT b) Second contrast: comparison between CBT and NT *Note: These contrasts are non-orthogonal. Here, that's ok because we only have one predictor variable (cf. Field, ch. 16.6.6: Setting Contrasts)* *Note: We don't have to specify contrasts manually. Default contrasts would take `BT` as reference group and the difference to group `CBT` would be one of our effects.* 2. Create the nevessary `outcome` object by using `cbind()` to join `ocd_data$thoughts` and `ocd_data$actions` together. 3. Use the `manova()` command to run the analysis. Save the result as an object. 4. Use the `summary()`command on the MANOVA object, including the additional argument `intercept = TRUE`. 5. What conclusions can you draw from the output? `r if(!params$soln) {"\\begin{comment}"}` ### Solution #### Subtask 1 ```{r} .bt.vs.nt <- c(-1,1,0) .cbt.vs.nt <- c(-1,0,1) # Attach contrasts to factor contrasts(ocd_data$group) <- cbind(.bt.vs.nt, .cbt.vs.nt) mean(ocd_data$actions) ocd_data %>% group_by(group) %>% summarize(mean(actions)) ``` #### Subtask 2 ```{r} outcome <- cbind(ocd_data$actions, ocd_data$thoughts) ``` #### Subtask 3 ```{r} model1 <- manova(outcome ~ group, data = ocd_data) ``` #### Subtask 4 ```{r} summary(model1, intercept = TRUE) ``` #### Subtask 5 The F test for the groups is significant (F(4,54) = 2.56, p = .049). This means that the tyoe of therapy received had an influence on the frequency of obsessive symptoms, measured through both thoughts and actions. From these results alone, this is the most detailed conclusons we're able to draw right now. `r if(!params$soln) {"\\end{comment}"}` ## 4) Interpreting the MANOVA 1. Use the command `summary.aov()` on the MANOVA model. 2. Briefly interpret the results. a) Do these results justify analysing the planned contrasts? 3. Create a separate ANOVA model for each of the outcome variables and investigate the output in relation to your contrasts. 4. What do you conclude from these results? `r if(!params$soln) {"\\begin{comment}"}` ### Solution #### Subtask 1 ```{r} summary.aov(model1) ``` #### Subtask 2 Neither of the separate ANOVAs yield significant results. This is an indication for neither obsessive thoughts nor obsessive actions being especially influenced by the kind of therapy received. That's highly interesting because we do have found an effect on the combination of these two components. Based on this, we don't really have a good reason for investigating our planned contrasts. For the sake of this sheet, we'll do it anyway. #### Subtask 3 ```{r} out1 <- lm(actions ~ group, data = ocd_data) out2 <- lm(thoughts ~ group, data = ocd_data) summary(out1) summary(out2) ``` #### Subtask 4 For obsessive actions, the behavioral therapy works singificantly better than the no treatment control. However, for obsessive thoughts, it's the cognitive behavioral therapy beating the control group. `r if(!params$soln) {"\\end{comment}"}` ## 5) Rendering (knit) Render this file using `Ctrl` + `Shift` + `K` (Windows) or `Cmd` + `Shift` + `K` (Mac). In the viewer you should see a pretty versionn of your file. If this works: Congratulations! Your code can be rendered completely and without error codes! If it doesn't: No worries, you'll get there! Go hunting for errors in your code! Otherwise, we'll get it to render in the next seminar session! ## Literature *Note*: These sheets are based partially on exercises from the book *Dicovering Statistics Using R* (Field, Miles & Field, 2012). They've been modified for the porpuses of this seminar, and the R code was updated. Field, A., Miles, J., & Field, Z. (2012). *Discovering Statistics Using R*. London: SAGE Publications Ltd.