--- title: 'Übungszettel Anova' author: "Joshua Driesen (joshua.driesen@stud.uni-goettingen.de)" subtitle: "M.Psy.205 Multivariate Statistik. Dozent: Peter Zezula" output: html_document: theme: paper toc: yes toc_depth: 3 toc_float: yes pdf_document: default #latex_engine: xelatex header-includes: - \usepackage{comment} params: soln: TRUE mainfont: Arial linkcolor: blue --- # Deutsch ## Links [Übungszettel als PDF-Datei zum Drucken](https://pzezula.pages.gwdg.de/sheet_anova.pdf) `r if(!params$soln) {"\\begin{comment}"}` **Übungszettel mit Lösungen** [Lösungszettel als PDF-Datei zum Drucken](https://pzezula.pages.gwdg.de/sheet_anova_solutions.pdf) [Der gesamte Übugszettel als .Rmd-Datei](https://pzezula.pages.gwdg.de/sheet_anova.Rmd) (Zum Downloaden: Rechtsklick > Speichern unter...) `r if(!params$soln) {"\\end{comment}"}` ## Hinweise zur Bearbeitung 1. Bitte 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, Kap. 10|Eine ausführliche, sehr gute Einführung in die einfaktorielle ANOVA (und Kontraste) Field, Kap. 12|Für mehrfaktorielle Designs und Varianzanaylsen ## Tipp der Woche ### Befehle zu Paketen Die in unseren Skripten benutzten Befehle require() und library() akzeptieren den Namen der Pakete, auf die sie sich beziehen, einfach so, also z.B. require(dplyr). Wenn Sie auf Ihrem eigenen Rechner Pakete zum ersten Mal installieren wollen, erwartet Sie jedoch ein Stolperstein: Der install.packages()-Befehl verlangt den Paketnamen in Anführungszeichen (also als sogenannter String), z.B. install.packages("dplyr"). ## 1) Einfaktorielle ANOVA 1. Setzen Sie, wie Sie es bereits gewohnt sind, den von Ihnen gewählten Ordner als Arbeitsverzeichnis, und laden Sie die tidyverse-Pakete. 2. Lesen Sie den Datensatz `Superhero.dat` mit dem Befehl `read_delim()` direkt aus der URL `https://pzezula.pages.gwdg.de/data/Superhero.dat` in R ein. Dieser enthält Daten von Kindern, die in einem Superheldenkostüm in einer Notaufnahme vorstellig wurden. Dabei ist sowohl die Schwere der Verletzung vermerkt ("injury"), als auch, welchen Superheld die Kostüme nachbildeten ("hero"). Die Kodierung der Helden erfolgte nach folgendem Schema: 1 - Spiderman 2 - Superman 3 - Hulk 4 - Teenage Mutant Ninja Turtles Wir wollen nun mittels ANOVA herausfinden, ob verschiedene Superheldenkostüme mit unterschiedlich schweren Verletzungen assoziiert sind. 3. Ein FALSCHER Ansatz: Nutzen Sie den lm()-Befehl, um den Verletzungsgrad aus der Heldenkategorie vorherzusagen, und betrachten Sie das Ergebnis. Warum ist diese Art der Analyse NICHT valide? 4. Nutzen Sie den factor()-Befehl, um die Heldenvariable zu faktorisieren. Nutzen Sie auch labels(), um mit der Tabelle oben die Helden ihren jeweiligen Codes zuzuordnen. 5. Erzeugen Sie einen Bar-Plot mit Fehlerbalken. 6. Die unterschiedlich langen Fehlerbalken könnten auf eine Verletzung der Varianzhomogenität hinweisen. Prüfen Sie diese Vermutung mit dem Levene-Test. (Nutzen Sie hierzu mit "??levene" R's Suchfunktion, um den entsprechenden Befehl zu finden.) 7. Führen Sie mit dem aov()-Befehl eine ANOVA durch, und speichern Sie das Ergebnis als ANOVA1. Unterscheiden sich die verschiedenen Superhelden hinsichtlich der Verletzungsschwere? 8. In der Vorlesung wurde die Äquivalenz von ANOVA und Regression besprochen. Führen Sie dementsprechend die gleiche ANOVA nocheinmal durch, diesmal jedoch mit dem lm()-Befehl. Speichern Sie das Ergebnis als ANOVA2. Vergleichen Sie die beiden Ergebnisobjekte miteinander. Welche Werte entsprechen sich? Welche Analyseform beeinhaltet mehr Information? Bonus: Was genau bedeutet heroHulk = -6.250? 9. Führen Sie eine weitere Analyse durch, in der Sie untersuchen, ob sich die beiden bekannten Comic-Universen Marvel und DC hinsichtlich der assoziierten Verletzungsschwere unterscheiden. Definieren Sie hierzu mit contrasts(Superhero$hero) <- IhrKontrast einen geplanten Kontrast. Spiderman und Hulk gehören zu Marvel, Superman entstammt dem DC-Universum, und die Ninja Turtles gehen nicht in diese Analyse ein. `r if(!params$soln) {"\\begin{comment}"}` ### Lösung #### Unteraufgabe 1 ```{r} library(tidyverse) ``` #### Unteraufgabe 2 ```{r} Superhero <- read_delim("https://pzezula.pages.gwdg.de/data/Superhero.dat", delim = "\t") ``` #### Unteraufgabe 3 ```{r} Falsch <- lm(injury ~ hero, data = Superhero) summary(Falsch) ## die Heldenvariable wird als numerische Variable ausgewertet. ``` Diese Auswertung interpretiert die Zahlencodes der Helden als intervallskalierten Prädiktor. Solange Superman aber nicht um genausoviel heldiger als Spiderman ist, wie die Ninja Turtles als Hulk, ist dies aber Unsinn. #### Unteraufgabe 4 ```{r} # we define a new variable hero.f with factorized categorial levels and label them also Superhero$hero.f <- factor(Superhero$hero, levels = c(1:4), labels = c("Spiderman", "Superman", "Hulk", "Ninja Turtle") ) Superhero$hero.f ``` #### Unteraufgabe 5 ```{r} # library(Hmisc) bar <- ggplot(Superhero, aes(hero.f, injury)) bar + stat_summary(fun = mean, geom = "bar") + stat_summary(fun.data = mean_cl_normal, geom = "errorbar", width = 0.2, size = 1) + labs(x = "Superheld", y = "Verletzungsschwere") ``` #### Unteraufgabe 6 ```{r} ??levene car::leveneTest(Superhero$injury, Superhero$hero.f) ``` Nicht signifikant, auf Varianzhomogenität wird geschlossen #### Unteraufgabe 7 ```{r} ANOVA1 <- aov(injury ~ hero.f, Superhero) summary(ANOVA1) ``` Signifikant, die Superhelden unterscheiden sich hinsichtlich der Verletzungsschwere #### Unteraufgabe 8 ```{r} ANOVA2 <- lm(injury ~ hero.f, Superhero) summary(ANOVA2) ``` Die F-statistic in ANOVA2 entspricht ANOVA1. Zusätzlich enthält ANOVA2 Informationen über die verwendeten Dummy-Variablen, sodass gerichtete Aussagen über einzelne Mittelwerte möglich werden. Bonus: heroHulk = -6.250 bedeutet, dass Kinder im Hulkkostüm im Durchschnitt um 6.25 Einheiten weniger schlimme Verletzungen hatten als Kinder in der Referenzgruppe unserer Dummykodierung, also Kinder im Spidermankostüm. #### Unteraufgabe 9 ```{r} contrasts(Superhero$hero.f) <- c(1, -2, 1, 0) MarvelVsDC <- lm(injury ~ hero.f, Superhero) summary.lm(MarvelVsDC) ``` Der erste Kontrast hero1 wird signifikant, das negative Vorzeichen zeigt, dass DC mit schwereren Verletzungen assoziiert ist. Da wir nur einen der drei möglichen orthogonalen Kontraste definiert haben, füllt R die beiden anderen mit zu unserem Kontrast und zueinander orthogonalen, anderen Kontrasten auf. Da unsere Hypothese aber schon vom ersten Kontrast erschöpfend geprüft wird, können wir die Kontraste hero2 und hero3 ignorieren. `r if(!params$soln) {"\\end{comment}"}` ## 2) Zweifaktorielle ANOVA 1. Lesen Sie den Datensatz `ChickFlick.dat` mit dem Befehl `read_delim()` direkt aus der URL `https://pzezula.pages.gwdg.de/data/ChickFlick.dat` in R ein. In diesem finden Sie die physiologischen Erregungsmessungen von Männern und Frauen, die im Labor entweder den klassischen "Chick-Flick" Bridget Jones' Diary, oder den Thriller Memento zu sehen bekamen. 2. Erzeugen Sie einen Bar-Plot mit Fehlerbalken. 3. Prüfen Sie die Validität des Konzeptes "Chick-Flick", indem Sie mittels ANOVA beantworten, ob die beiden Geschlechter unterschiedlich auf verschiedene Filme reagieren. 4. Wieviel Vertrauen haben Sie in die Generalisierbarkeit Ihrer Schlussfolgerung aus der vorherigen Aufgabe? Welche Alternativerklärungen bleiben offen? Wie liesse sich die Aussagekraft erhöhen? `r if(!params$soln) {"\\begin{comment}"}` ### Lösung #### Unteraufgabe 1 ```{r} ChickFlick <- read_delim("https://pzezula.pages.gwdg.de/data/ChickFlick.dat", delim = "\t") # we set factor ChickFlick$gender <- factor(ChickFlick$gender) ChickFlick$film <- factor(ChickFlick$film) ``` #### Unteraufgabe 2 ```{r} # library(Hmisc) bar <- ggplot(ChickFlick, aes(gender, arousal, group = film, fill = film)) bar + stat_summary(fun = mean, geom = "bar", position = "dodge") + stat_summary(fun.data = mean_cl_normal, geom = "errorbar", position = position_dodge(width = 0.90), width = 0.2) + labs(x = "Gender", y = "Mean Arousal", fill = "Film") ``` #### Unteraufgabe 3 ```{r} ChickFlickANOVA <- aov(arousal ~ gender * film, ChickFlick) summary(ChickFlickANOVA) ``` Der Interaktionsterm wird nicht signifikant, unser Experiment kann das Konzept, dass Männer und Frauen unterschiedlich auf verschiedene Filme reagieren, nicht untermauern. #### Unteraufgabe 4 Für ein so schwammiges Konzept wie "Chick-Flick" benötigen wir valideres, und vor allem mehr Stimulusmaterial. Die beiden gewählten Filme scheinen zwar recht gute Prototypen ihrer Kategorien zu sein, jedoch ist nicht auszuschliessen, dass etwa Bridget Jones' Diary einfach insofern ein schlechter Chick-Flick ist, dass er es nicht schafft, die frauenspezifische Erregung dieses Genres zu erzeugen. `r if(!params$soln) {"\\end{comment}"}` ## 3) Rendern Rendern, bzw. knitten Sie nun das Dokument über die Tastenkombination `strg` + `shift` + `k` (Windows) oder `cmd` + `shift` + `k`. Wenn das funktioniert: Top gemacht! Wenn nicht: Schauen Sie sich die Fehlermeldung an, und betrachten Sie insbesondere die Zeilen Ihrer Syntax, die in der Fehlermeldung auftauchen. Suchen Sie nach dem Fehler und probieren Sie es erneut! ## 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://pzezula.pages.gwdg.de/sheet_anova.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_anova_solutions.pdf) [The source code of this sheet as .Rmd](https://pzezula.pages.gwdg.de/sheet_anova.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://pzezula.pages.gwdg.de/students_template.Rmd) by donloading it. 2. You may find the informations useful that you can find on the [start 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 to 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 This is a hands on course. We cannot present you all the useful commands in detail. Instead we give you links to useful ressources, where you might find hints to help you with the exercises. Ressource | Description ----------|-------------- Field, Chapter 10 | A detailed introduction to single factor ANOVA and contrasts Field, Chapter 12 | This chapter deals with multi factorial designs and variance analyses ## Tip of the week ### Dealing with package names We use two commands to load packages in our scripts, `require()` and `library()`, both of which understand quoted names of packages, f. e. `require("psych")`. To make life a bit easier, we can also use unquoted package names, f. e. `require(psych)` or `library(dplyr)`. But take care: We still need to use quoted package names when we want to install using `install.packages()` f. e. `install.packages("psych")`. ## 1) Single factor ANOVA (oneway) 1. Set your working directory and load tidyverse as usual. 2. Read the data file 'Superhero.dat' using either `read_delim()` or `read_tsv()` from the URL `https://pzezula.pages.gwdg.de/data/Superhero.dat` and store it in a data object in R. We find data of injured children that were dressed up like a superhero when they were presented in an emergency room of a hospital. 'injury' gives us severity of the injury and 'hero' specifies, which superhero was beeing copied: 1 - Spiderman 2 - Superman 3 - Hulk 4 - Teenage Mutant Ninja Turtles Using a ANOVA we want to find out, whether there is an association between type of superhero and severiy of injuries. 3. A wrog approach: Use the command `lm()` and take a look at the results. Why is this analysis **not** valid? 4. Use the command `factor()` to factorize variable `superhero`. Also define `labels()` to connect type of hero and the referred code. 5. Generate a barplot with errorbars. 6. The varying errorbars might indicate inhomogeneous variances. Test this using the levene test. Try `??levene` to let R help you to learn about this command. 7. Use the command `aov()` to get an ANOVA and store the result under the name 'ANOVA1'. Are there significant differences between the superheroes severity of injuries? 8. You heared of equivalence of ANOVA and regression in the lecture. Adapt the same ANOVA to the data using the command `lm()` and store the result under 'ANOVA2'. Compare the two result objects. Can you identify corresponding results? Which analysis provides more information? Bonus: what exactly is meant by "heroHulk = -6.250"? 9. Do further analyses where you find out, whether the two comic universes "Marvel" and "DC" differ in severity of injuries. Do this by applying planned contrasts and use `contrasts(Superhero$hero) <- ` for this. Spiderman and Hulk are members of Marvel, Superman is in the DC universe and Ninja Turtles are not analysed here. `r if(!params$soln) {"\\begin{comment}"}` ### Solutions #### Subtask 1 ```{r} library(tidyverse) ``` #### Unteraufgabe 2 ```{r} Superhero <- read_delim("https://pzezula.pages.gwdg.de/data/Superhero.dat", delim = "\t") # or Superhero <- read_tsv("https://pzezula.pages.gwdg.de/data/Superhero.dat") ``` #### Subtask 3 ```{r} Wrong <- lm(injury ~ hero, data = Superhero) summary(Wrong) ``` In this calculation we wold interprete group codes as continuous variable (interval-scaled). This is nonsense. #### Subtask 4 ```{r} Superhero$hero.f <- factor(Superhero$hero, levels = c(1:4), labels = c("Spiderman", "Superman", "Hulk", "Ninja Turtle") ) Superhero$hero.f ``` #### Subtask 5 ```{r} # library(Hmisc) bar <- ggplot(Superhero, aes(hero.f, injury)) bar + stat_summary(fun.y = mean, geom = "bar") + stat_summary(fun.data = mean_cl_normal, geom = "errorbar", width = 0.2, size = 1) + labs(x = "Superhero", y = "Severity of Injuries") ``` #### Subtask 6 ```{r} # we wold enter # ??levene # but this doesnt make sense in a Rmd that will be rendered # so we comment it out # leveneTest() is part of package car # we call it without loading it car::leveneTest(Superhero$injury, Superhero$hero.f) ``` No significant differences of variances. So we assume homogenity of variances. #### Subtask 7 ```{r} ANOVA1 <- aov(injury ~ hero.f, Superhero) summary(ANOVA1) ``` We find significant differences of the severity of injuries between the types of superheroes. #### Subtask 8 ```{r} ANOVA2 <- lm(injury ~ hero.f, Superhero) summary(ANOVA2) ``` The F-statistics correspond in ANOVA1 and ANOVA2. ANOVA2 gives us additional informations about the used dummy variables. So there are some specific tests or comparisons between pairs of means. Bonus: `heroHulk = -6.250` tells us, that children dressed like Hulk had 6.25 units less severe injuries than the reference group of our dummy coding system. Here the reference group would be children dressed up like Spidermen. #### Subtask 9 ```{r} contrasts(Superhero$hero.f) <- c(1, -2, 1, 0) MarvelVsDC <- lm(injury ~ hero.f, Superhero) summary.lm(MarvelVsDC) ``` The first contrast `hero1` is highly significant. The '-' shows, that DC is associated with more severe injuries. As we defined only one of the three possible orthogonal contrasts, R completes the two missing orthogonal contrasts. But our hypothesis is already tested by our first contrast, so we can ignore the results of the other ones (hero2 and hero3). `r if(!params$soln) {"\\end{comment}"}` ## 2) Two factorial ANOVA 1. Read the data `ChickFlick.dat` using either `read_delim()` or `read_tsv()` from the URL `https://pzezula.pages.gwdg.de/data/ChickFlick.dat` and store it in a data object in R. You find data of physiological activation of men and women, that saw the classical "Chick-Flick" movie "Bridget Jones' Diary" or the Thriller "Memento". 2. Generate a bar plot with error bars included. 3. Check the validity of the concept "chick-flick" by finding out whether we have different reactions to the films due to gender. 4. Do you trust in the generalizibility of the results of the above computation? Do you see alternative explanations? How could we make our results more trustworthy? `r if(!params$soln) {"\\begin{comment}"}` ### Solutions #### Subtask 1 ```{r} ChickFlick <- read_delim("https://pzezula.pages.gwdg.de/data/ChickFlick.dat", delim = "\t") # or ChickFlick <- read_tsv("https://pzezula.pages.gwdg.de/data/ChickFlick.dat") # we set factor ChickFlick$gender <- factor(ChickFlick$gender) ChickFlick$film <- factor(ChickFlick$film) ``` #### Subtask 2 ```{r} # library(Hmisc) bar <- ggplot(ChickFlick, aes(gender, arousal, group = film, fill = film)) bar + stat_summary(fun.y = mean, geom = "bar", position = "dodge") + stat_summary(fun.data = mean_cl_normal, geom = "errorbar", position = position_dodge(width = 0.90), width = 0.2) + labs(x = "Gender", y = "Mean Arousal", fill = "Film") ``` #### Subtask 3 ```{r} ChickFlickANOVA <- aov(arousal ~ gender * film, ChickFlick) summary(ChickFlickANOVA) ``` We cannot prove our interaction term be significant. So we don't have evidence in our data that men and women react different to the movies in question. #### Subtask 4 We really need more and especially more valid stimulus material to get clearer about a concept as nebulous as "Chick-Flick". Although both movies seem to be good prototypes of teir category, "Bridget Jones' Diary" might not cause the specific female activation pattern that this genre is supposed to produce. `r if(!params$soln) {"\\end{comment}"}` ## 3) Render Render (or knit) your Rmd file using the shortkey `strg` + `shift` + `k` (Windows) or `cmd` + `shift` + `k` (Mac). If it works, well done! If not, check your error messages, inspect the lines of your code, where the error is supposed to occur. Correct the error and over again. ## 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. Version: `r format(Sys.time(), '%d %B, %Y %H:%M')`