---
title: 'Übungszettel Logistische Regression, Exercise Sheet Logistic Regression'
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: 2
    toc_float: yes
  pdf_document: default
  #  latex_engine: xelatex
header-includes:
 - \usepackage{comment}
params:
 soln: TRUE
mainfont: Helvetica
linkcolor: blue
date: "`r format(Sys.time(), '%d %B, %Y %H:%M')`"
---

# Deutsche Version

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

## Links

[Übungszettel als PDF-Datei zum Drucken](https://pzezula.pages.gwdg.de/sheet_log_reg.pdf)

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

[Lösungszettel als PDF-Datei zum Drucken](https://pzezula.pages.gwdg.de/sheet_log_reg_solutions.pdf)

[Der gesamte Übugszettel als .Rmd-Datei](https://pzezula.pages.gwdg.de/sheet_log_reg.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 8 | Buchkapitel, das Schritt für Schritt erklärt, worum es geht, und wie man logistische Regressionen in R durchführt. **Große Empfehlung!**


## Tipp der Woche
Wenn Sie Zeilenumbrüche in Ihrem Code verwenden, ist er für Sie und Andere deutlich besser lesbar. Sie können z.B. nach jeder Pipe (` %>% `) oder nach dem `+` in ggplot2-Befehlen sinnvolle Zeilenumbrüche einbauen. Zögern Sie auch nicht, längere einzelne Befehle in mehrere Zeilen umzubrechen: Nach jedem Komma kann eine neue Zeile beginnen, wenn Sie es möchten.

Es kann sein, dass die Einrückung bei mehrzeiligen Befehlen irgendwann unordentlich wird. Wenn das passiert, markieren Sie den Code und drücken `strg` + `i` (windows), bzw. `cmd` + `i` (mac). Ihr Code wird dadurch autamtisch eingerückt. 

## 1 Daten einlesen

1. Laden Sie das `tidyverse` und fügen Sie die entsprechende Code-Zeile an den Anfang Ihrer .Rmd-Datei ein.
2. Setzen Sie ein sinnvolles Arbeitsverzeichnis und fügen Sie die entsprechende Code-Zeile an den Anfang Ihrer .Rmd-Datei ein.
3. Laden Sie den Datensatz `stats_lectures` unter [diesem Link](http://md.psych.bio.uni-goettingen.de/mv/data/div/stats_lectures.csv) herunter und speichern ihn in Ihrem Arbeitsverzeichnis, dort z.B. in einem Unterordner `data`.
4. Lesen Sie den Datensatz in R ein.

### Erklärung der Variablen

Jede Zeile enthält Daten von einem/r StudentIn. Beobachtet wurde, ob die Studierenden an Statistik-Lehrveranstaltungen teilnehmen. Erhoben wurden dazu die Temperatur, wie sehr die Studis Sonne genießen und wie sehr sie Statistik genießen. Letzteres wurde mit einem Fragebogen gemessen, der aus 5 Items besteht.

**Hinweis: Die Daten sind vollständig fiktiv und spiegeln natürlich auch nicht die Meinung der Lehrkräfte wieder. Dem Tutor ist nur gerade kein besseres Beispiel eingefallen.**

Name | Bedeutung
-------------|--------------------------------------
attend_or_not | Hat zwei Ausprägungen: "course not attended" und "course attended". Die Ausprägungen geben an, ob der/die Studi zur Lehrveranstaltung gegangen ist.
temperature | Die Temperatur. Höher ist heißer.
sun_joy | Wie sehr genießen die Studis Sonne. Höher ist stärkerer Genuss
stats_joy1 bis stats_joy5 | Wie sehr genießen die Studis Statistik. Höher ist stärkerer Genuss.

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

#### Unteraufgabe 1
```{r}
library(tidyverse)
```

#### Unteraufgabe 2
Beachten Sie hier, dass in R alle Ordner mit `/` getrennt werden. Deshalb können Sie Windows-Dateipfade nicht einfach kopieren, diese sind mit `\` getrennt. Ersetzen Sie dort `\` durch `/`.

Der Code ist:
```{r, eval = FALSE}
setwd()
```

Hier ein Beispiel mit meinem Dateipfad.
```{r, eval = FALSE}
# Beispiel
setwd("~/ownCloud/_Arbeit/Hiwi Peter/gitlab_sheets/public")
```

#### Unteraufgabe 3

Bitte folgen Sie den Anweisungen in der Aufgabenstellung.

#### Unteraufgabe 4
```{r}
#stats_lectures <- read_csv("public/mv/data/stats_lectures.csv")
stats_lectures <- read_csv("http://md.psych.bio.uni-goettingen.de/mv/data/div/stats_lectures.csv")
```
`r if(!params$soln) {"\\end{comment}"}`

## 2 Datenaufbereitung

1. Bilden Sie für jede Person im Datensatz den Mittelwert aus den 5 Items stats_joy1 bis stats_joy5 und fügen Sie diesen als neue Variable dem Datensatz hinzu. *Hinweis: Eine Google-Suche kann hier sehr hilfreich sein. Achten Sie auch auf Kleinigkeiten in den Befehlen, die Sie finden.* Falls Sie diese Aufgabe nicht lösen können, aber dennoch weiter machen möchten, finden Sie hier funktionierende Syntax, die Sie in Ihr Script kopieren können. [Link](https://pzezula.pages.gwdg.de/data/06_sheet_task2-1.R)
2. Nutzen Sie den Befehl `factor()`, um die Variable `attend_or_not` in einen Faktor umzuformen. Verwenden Sie das Argument `levels = c()`, um die Faktorstufen anzugeben. So stellen Sie sicher, dass Sie wissen, welche Kategorie die Baseline darstellt. Falls Sie diese Aufgabe nicht lösen können, aber dennoch weiter machen möchten, finden Sie hier funktionierende Syntax, die Sie in Ihr Script kopieren können. [Link](https://pzezula.pages.gwdg.de/data/06_sheet_task2-2.R)

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

Der Befehl `mean()` muss hier mit `mean(c())` verwendet werden, damit die Zeilenweisen Mittelwerte richtig gebildet werden.
Wir verwenden mutate(), um eine neue Varibale zu erstellen.
```{r}
stats_lectures <- stats_lectures %>% 
  rowwise() %>% # wir wollen zeilenweise vorgehen
  mutate(stats_joy = mean(c(stats_joy1, stats_joy2, # mittelwert bilden
                            stats_joy3, stats_joy4, 
                            stats_joy5)))
```

#### Unteraufgabe 2

```{r}
stats_lectures <- stats_lectures %>% mutate(attend_or_not = factor(attend_or_not, levels = c("course not attended", "course attended")))
```
`r if(!params$soln) {"\\end{comment}"}`

## 3 Erste logistische Regression

1. Nutzen Sie den Befehl `glm()` mit dem Argument `family = binomial()`, um die Kursbesuche der Studierenden mit deren Genuss von Statistik vorherzusagen. Benutzen Sie dafür den Mittelwert, den Sie in der vorherigen Aufgabe gebildet haben, als Prädiktor.
2. Lassen Sie sich den Output mit `summary()` anzeigen. "Null deviance" gibt die Deviance-Statistik für das Null-Modell an, "Residual Deviance" für das Alternativmodell (user spezifiziertes). Was können Sie aus der Deviance ableiten? (Siehe Field, Kapitel 8.3.2: Assessing the model: the deviance statistic)
3. Bilden Sie mit `exp()` das Exponential zur Basis e für die Regressionskoeffizienten. Dies ist das Odds-Ratio. *Tipp: Mit `$coefficients` können Sie direkt auf die Koeffizienten zugreifen, wenn Sie den Namen, den Sie dem Modell gegeben haben, vor das $-Zeichen schreiben.*
4. Interpretieren Sie das Odds Ratio für den vorliegenden Fall. (Siehe Field, Kapitel 8.3.6: The odds ratio)
5. Wenden Sie `confint()` auf das Modell an, um die Konfidenzintervalle für die Prädiktoren zu erhalten. Bilden Sie auch hiervon das Exponential mit `exp()`, um die Konfidenzinternalle für die Odds-Ratios zu erhalten.
6. Erstellen Sie einen Plot, um Ihr Modell zu visualisieren. 
    a) Wenden Sie den Befehl `fitted()` auf das Modell an, um die durch das Modell vorhergesagten Wahrscheinlichkeiten für Ihre Studierenden zu erhalten.
    b) Fügen Sie diese Werte Ihrem Datensatz hinzu.
    c) Nutzen Sie `gpplot`, um einen Plot zu erstellen. Verwenden Sie auf der X-Achse stats_joy und auf der y-Achse die in a) und b) erstellten vorhergesagten Wahrscheinlichkeiten. 

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

```{r}
model1 <- glm(attend_or_not ~ stats_joy, 
              data = stats_lectures, family = binomial())
```

#### Unteraufgabe 2

```{r}
summary(model1)
```

Die Deviance ist ein Indikator dafür, wie viel Information nicht durch das Modell aufgeklärt werden kann. Eine größere Deviance bedeutet, dass weniger Information erklärt werden kann. Mehr dazu können Sie hier nachlesen: Field, Kapitel 8.3.2: Assessing the model: the deviance statistic.

#### Unteraufgabe 3

```{r}
exp(model1$coefficients)
```

#### Unteraufgabe 4

Wenn hier bei Studi A die Freude an Statistik um einen Punkt höher ist als bei Studi B, dann heißt das, dass die Chance, dass Studi A die Statistik-Vorlesung besucht, ca. 2.86-mal größer ist, als dass Studi B es tut.

D.h. eine Steigerung der Freude an Statistik um einen Punkt führt zu einer 2.86-mal höheren Chance, die Statistik-Vorlesung zu besuchen.

#### Unteraufgabe 5

```{r}
exp(confint(model1))
```

#### Unteraufgabe 6

a)
Ich nutze hier head(), um nur die ersten 6 Werte anzeigen zu lassen. Andernfalls würde das Dokument sehr voll werden.
```{r}
head(fitted(model1))
```

b)
```{r}
stats_lectures$pred_prob <- fitted(model1)
```

c)
```{r}
# Objekt erstellen
plot1 <- ggplot(stats_lectures, aes(x = stats_joy, y = pred_prob))

# Plot befüllen
plot1 <- plot1 + geom_point() +
  labs(x = "Freude an Statistik",
       y = "Wahrscheinlichkeit Vorlesungsbesuch",
       title = "Vorlesungsbesuche in Abhängigkeit von Freude an Statistik") +
  theme_classic()

# plot anzeigen
plot1
```

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


## 4 Modellvergleich logistische Regression

1. Erstellen Sie drei weitere Modell mit dem `glm()`-Befehl:
    a) Fügen Sie den ersten Modell den Prädiktor temperature hinzu
    b) Fügen Sie dem Modell aus a) den Prädiktor sun_joy hinzu
    c) Fügen Sie dem Modell aus b) die Interaktion aus temperature und sun_joy als Prädiktor hinzu.
2. Testen Sie mit dem Befehl `anova()`, welches der Modelle am besten auf die Daten passt. Nutzen Sie die Option `test = "Chisq"`, um einen p-Wert zu erhalten.
    a) Sehen Sie sich mit `summary()` den Output für das beste Modell an.
    b) Welche Hypothesen würden Sie aus dieser Analyse ableiten?
3. Berechnen Sie für das beste Modell die Odds Ratios der Prädiktoren und deren Konfidenzintervalle.
4. Berechnen Sie für das beste Modell die drei Pseudo-R^2 Indices (Siehe Field, Kapitel 8.3.3: Assessing the model: R and R^2). Sie können dafür die unten stehende Funktion nutzen. *Mehr dazu in Field, R's Souls' Tip 8.2.* 
    a) Kopieren Sie den Code und führen Sie in aus. 
    b) Sie können nun die Funktion `logisticPseudoR2s()` auf ihr Modell anweden. 

```{r}
logisticPseudoR2s = function(LogModel) { 
  dev = LogModel$deviance
  nullDev = LogModel$null.deviance
  modelN = length(LogModel$fitted.values) 
  R.l = 1 - dev / nullDev
  R.cs = 1 - exp(-(nullDev - dev) / modelN)
  R.n = R.cs / (1 - (exp(-(nullDev / modelN)))) 
  
  cat("Pseudo R^2 for logistic regression\n") 
  cat("Hosmer and Lemeshow R^2 ", round(R.l, 3), "\n") 
  cat("Cox and Snell R^2 ", round(R.cs, 3),"\n") 
  cat("Nagelkerke R^2 ", round(R.n, 3), "\n")
}
```

`r if(!params$soln) {"\\begin{comment}"}`
### Lösung
#### Unteraufgabe 1
```{r}
# das erste Modell
model1 <- glm(attend_or_not ~ stats_joy, 
              data = stats_lectures, family = binomial())

# a)
model2 <- glm(attend_or_not ~ stats_joy + temperature, 
              data = stats_lectures, family = binomial())

# b)
model3 <- glm(attend_or_not ~ stats_joy + temperature + sun_joy, 
              data = stats_lectures, family = binomial())

# c) (kann auf zwei Arten geschrieben werden)
model4 <- glm(attend_or_not ~ stats_joy + temperature*sun_joy, 
              data = stats_lectures, family = binomial())

model4 <- glm(attend_or_not ~ stats_joy + temperature + 
                sun_joy + temperature:sun_joy, 
              data = stats_lectures, family = binomial())
```

#### Unteraufgabe 2

```{r}
anova(model1, model2, model3, model4, test = "Chisq")
```

Modell 4 passt mit deutlichem Abstand am besten auf die Daten.

a)
```{r}
summary(model4)
```

b)
Man kann folgende Hypothesen ableiten:
1. Je höher die Freude an Statistik ist, desto höher ist die Wahrscheinlichkeit, dass Studierende Statistik-Vorlesungen besuchen.
2. Es gibt eine Interaktion zwischen Temperatur und Freude an der Sonne. Möglicherweise ist es so, dass eine hohe Temperatur bei SonnenliebhaberInnen dazu führt, dass die Statistik-Veranstaltungen weniger häufig besucht werden.

#### Unteraufgabe 3

```{r}
# Odds Ratios
exp(model4$coefficients)

# Konfidenzintervalle. exp() nicht vergessen!
exp(confint(model4))
```

#### Unteraufgabe 4

a) Siehe Aufgabenstellung.

b)
```{r}
logisticPseudoR2s(model4)
```
  
`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. 


```{r eval=FALSE, include=FALSE}
#-----------------------#
#--- DATENSIMULATION ---#
#-----------------------#

# Hier werden die Daten für diesen Übungszettel simuliert.

#-----------------------#
#---      PAKETE     ---#
#-----------------------#
library(tidyverse)
library(VGAM) # für Funktion logit() zum Umwandeln von logit in Wahrscheinlichkeiten.

#---------------------------#
#--- VARIABLEN ERSTELLEN ---#
#---------------------------#
set.seed(12) # seed für Reproduzierbarkeit

temperature <- rnorm(m = 5, sd = 3, n = 400) %>% # Temperatur
  scale() %>% # Z-Standardisiert für leichteres Erstellen des Outcomes
  as.double() # scale() macht eine Matrix aus dem Objekt, das machen wir rückgängig.

sun_joy <- rnorm(m = 5, sd = 3, n = 400) %>% # Freude an der Sonne
  scale() %>% # Z-Standardisiert für leichteres Erstellen des Outcomes
  as.double() # scale() macht eine Matrix aus dem Objekt, das machen wir rückgängig.

# 5 Items eines Fragebogens
stats_joy1 <- rnorm(m = 5, sd = 3, n = 400) %>% round(0)
stats_joy2 <- rnorm(m = 5, sd = 3, n = 400) %>% round(0)
stats_joy3 <- rnorm(m = 5, sd = 3, n = 400) %>% round(0)
stats_joy4 <- rnorm(m = 5, sd = 3, n = 400) %>% round(0)
stats_joy5 <- rnorm(m = 5, sd = 3, n = 400) %>% round(0)

# Zusammenführen der 5 Items in einen Datensatz
stats_data <- tibble(stats_joy1, stats_joy2, 
                     stats_joy3, stats_joy4, stats_joy5)

# Mittelwert der 5 Items pro Person
stats_joy <- rowMeans(stats_data) 


#---------------------------#
#---  OUTCOME ERSTELLEN  ---#
#---------------------------#
attending_logit = 1 + 2*stats_joy - 0.1*temperature - 0.1*sun_joy -
  1*temperature*sun_joy + rnorm(0, sd = 2, n = 1)

#---------------------------#
#---     AV ERSTELLEN    ---#
#---------------------------#

# Logit in Wahrscheinlichkeit umwandeln
attending_prob = logit(attending_logit, inverse = TRUE)

# Wahrscheinlichkeit nutzen, um 0 oder 1 zu erstellen
attend_or_not <- rbinom(400, size = 1, prob = attending_prob)

# Faktor bilden
attend_or_not <- factor(attend_or_not, levels = c(0, 1), labels = c("course not attended", "course attended"))

#---------------------------#
#---     SPEICHERN       ---#
#---------------------------#

# Daten zu tibble zusammenbinden
stats_lectures <- tibble(attend_or_not, temperature, sun_joy, stats_joy1, stats_joy2, stats_joy3, stats_joy4, stats_joy5)

# abspeichern
write_csv(stats_lectures, "data/stats_lectures.csv")
```

## 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_log_reg.pdf)

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

[Exercise sheet with solutions included as PD](https://pzezula.pages.gwdg.de/sheet_log_reg_solutions.pdf)

[The source code of this sheet as .Rmdi](https://pzezula.pages.gwdg.de/sheet_log_reg.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 8 | Book chapter explaining step by step the why and how of logistic regression in R. **Highly recommended!** 


## Hint of the week

You should use new lines in your code to make it better readable for yourself and for others.
You could f. e. continue in a new line after every pipe (` %>% `) 
or after a `+` in a sequence of ggplot2 commands.
Don't hesitate to divide a complex and long command into several lines.
After every comma (`,`) you can continue in a new line, if you like.

If your code looks scattered because of unregular indentation, 
you may mark a section of code and press `strg` + `i` (Windows), or `cmd` + `i`(Mac) 
to automatically indent your source code.


## 1 Read data

1. Load package `tidyverse` and insert the apropriate code at the beginning of your .Rmd document.
2. Change to an adequate working directory. 
   Take care: The render process assumes, that your working directory is set to the location of your .Rmd document.
3. Load the data `stats_lectures` from [this link](http://md.psych.bio.uni-goettingen.de/mv/data/div/stats_lectures.csv) 
   and store it in your working directory or in a subfolder there. You might want to call it `data`.
4. Read the data in R. Alternatively read the data directly from the URL above.


### Explanation of the variables
Each line has the data of one student.
You find, whether the student took part in a statistic lession, the current temperature, 
to which extend the person enjoys sun and, on the other hand, statistics.
This was measured using a 5 items questionnaire.

**Hint: These are completely fictive data. They don't reflect the docents opinion. 
This is just an spontaneous idea of a student helper.**

name                     | meaning
------------------------ | ---------------------------------------------
attend_or_not            | dichotomous variable:  "course not attended" and "course attended". The values code, whether the student attendet the course or not.
temperature              | high values code high temperatures
sun_joy                  | the higher the value the higher the joy
stats_joy1 to stats_joy5 | the higher the values the more the student enjoys statistics

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

#### Subtask 1
```{r}
library(tidyverse)
# or
require(tidyverse)
```


#### Subtask 2

Take care: R uses `/` to specify pathes. 
So you cannot simply copy and paste Windows pathes into R. 
You have to replace all `\` with `/`.
Again: Remember the special behavior of the rendering with regard to working directories.

The code would be: 
```{r, eval = FALSE}
setwd()
```

Here an example with a local path. 
Adapt that to your needs.

```{r, eval = FALSE}
# example
setwd("P:/R/mv")
```

#### Subtask 3

Please follow the instructions in subtask 3. 

#### Subtask 4

```{r}
# stats_lectures <- read_csv("data/stats_lectures.csv") # relative path, working directory has to be set adequately
# alternatively read it from the URL
stats_lectures <- read_csv("http://md.psych.bio.uni-goettingen.de/mv/data/div/stats_lectures.csv")
```
`r if(!params$soln) {"\\end{comment}"}`


## 2 Data preparation

1. Calculate the mean of items stats_joy1 to stats_joy5 and store it in the dataframe resp. tibble.
   *Hint: A google search might help. Take also care of the details of the commands, you find.*
   If you aren't able to solve this part of the task and want to continue, 
   you find a working syntax for this part under 
   [https://pzezula.pages.gwdg.de/data/06_sheet_task2-1.R](https://pzezula.pages.gwdg.de/data/06_sheet_task2-1.R)
   that you might integrate in your script.
2. Use the command `factor()` to transform variable `attend_or_not` into a factor.
   Use the argument `levels = c()` to specify factor levels.
   This is how you set the correct reference group or baseline.
   If you want to keep on and don't know how to do tis, you find a working syntax under
   [https://pzezula.pages.gwdg.de/data/06_sheet_task2-2.R](https://pzezula.pages.gwdg.de/data/06_sheet_task2-2.R)

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

The command `mean()` has to be done like  `mean(c())` to have the means calculated rowwise.
We use `mutate()` to create a new variable.

```{r}
stats_lectures <- stats_lectures %>% 
  rowwise() %>% # wir want rowwise calculation
  mutate(stats_joy = mean(c(stats_joy1, stats_joy2, # calculate means
                            stats_joy3, stats_joy4, 
                            stats_joy5)))
# we could alternatively use base systems rowMeans()
# {}
```

#### Subtask 2

```{r}
stats_lectures <- stats_lectures %>% mutate(attend_or_not = factor(attend_or_not, levels = c("course not attended", "course attended")))
```
`r if(!params$soln) {"\\end{comment}"}`


## 3 First logistic regression

1. Use the command `glm()` with argument `family = binomial()` to predict, 
   which students attend the statistics course and whether they enjoy it.
   The predictor is the mean, that you calculated recently.
2. Check the output using `summary()`. "Null deviance" is the deviance statistic for the null model, 
   "residual deviance" is the same for the user defined alternative model. 
   What can we learn from the deviance?
   C. f. Field (2012), chapter 8.3.2: Assessing the model: the deviance statistic.
3. Use `exp()` to calculate the exponet to base e for the regression coefficient.
   The result is the odds ratio.
   *Tip: `$coefficients` gives you direct access to the coefficients. 
   As always you have to put the name of your result model and a $ before.*
4. Interprete the odds ratio for your example.
   (See Field (2012), chapter 8.3.6: The odds ratio)
5. Apply `confint()` to your result model to get confidence intervals for your predictors.
   Also use `exp()` to get the confidence intervals for the odds ratios.
6. Make a plot to visualize your model.
   a) Apply `fitted()` to your model to get the predicted probailities for your student subjects.
   b) Add these values to your data set.
   c) Use `ggplot()` to make a plot. 
      Put `stats_joy` to the x-axes and on the y-axes the predicted probilities, you calculated in a) and b).

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

```{r}
model1 <- glm(attend_or_not ~ stats_joy, 
              data = stats_lectures, family = binomial())
```

#### Subtask 2

```{r}
summary(model1)
```
The deviance is an indicator of how much information can *not* be explained by our model.
The bigger the deviance, the less information can be explained.
Read more about this in Field (2012), chapter 8.3.2 Assessing the model: The deviance statistic.

#### Subtask 3

```{r}
exp(model1$coefficients)
```

#### Subtask 4

If student A has a stats_joy of 1 point higher than student B, 
this means, that the chance, that student A visits the statistics class is 2.86 times higher than for student B.

In other words: An increase of joy in statistics of 1 point results in a 2.86 times higher probability to visit the statistics class. 

#### Subtask 5

```{r}
exp(confint(model1))
```

#### Subtask 6

a) We use `head()`, to see the first 6 values only. 
   Thus we avoid an extremely large document.

```{r}
head(fitted(model1))
```

b)
```{r}
stats_lectures$pred_prob <- fitted(model1)
```

c)
```{r}
# Create a result object
plot1 <- ggplot(stats_lectures, aes(x = stats_joy, y = pred_prob))

# Fill the plot
plot1 <- plot1 + geom_point() +
         labs(x = "Fun with statistics",
              y = "Probability of visiting the class",
         title = "Relation between fun with statistics and visit of statistics classes") +
         theme_classic()

# Show plot
plot1
```

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


## 4 Model comparison in logistic regression

1. Adapt three more models using the command `glm()`:
   a) Add predictor `temperature`
   b) Add predictor `sun_joy`
   c) Add to the model in b) the interaction of `temperature` and `sun_joy`.
2. Compare the models adapted using the command `anova()`. 
   Which of the models fit the data best, considering parsimonia.
   Use the parameter `test = "Chisq"` to get a p-value.
   a) Get the `summary()` of the best model.
   b) Which hypothesis would you deduct from this analysis.
3. Calculate the odds ratios and the confidence intervales for the parameters of the best model.
4. Calculate the pseudo R^2 for the best model (c. f. Field (2012), chapter 8.3.3: Assessing the model: R and R^2).
   You may use the function below for this. 
   *Read more on that in Field (2012), R's Souls' Tip 8.2.* 
   a) Copy the code and execute it
   b) You can apply `logisticPseudoR2s()` to your model.


```{r}
logisticPseudoR2s = function(LogModel) { 
  dev = LogModel$deviance
  nullDev = LogModel$null.deviance
  modelN = length(LogModel$fitted.values) 
  R.l = 1 - dev / nullDev
  R.cs = 1 - exp(-(nullDev - dev) / modelN)
  R.n = R.cs / (1 - (exp(-(nullDev / modelN)))) 
  
  cat("Pseudo R^2 for logistic regression\n") 
  cat("Hosmer and Lemeshow R^2 ", round(R.l, 3), "\n") 
  cat("Cox and Snell R^2 ", round(R.cs, 3),"\n") 
  cat("Nagelkerke R^2 ", round(R.n, 3), "\n")
}
```

`r if(!params$soln) {"\\begin{comment}"}`
### Solution
#### Subtask 1
```{r}
# the first model
model1 <- glm(attend_or_not ~ stats_joy, 
              data = stats_lectures, family = binomial())

# a)
model2 <- glm(attend_or_not ~ stats_joy + temperature, 
              data = stats_lectures, family = binomial())

# b)
model3 <- glm(attend_or_not ~ stats_joy + temperature + sun_joy, 
              data = stats_lectures, family = binomial())

# c) (we could write that in two manners)
model4 <- glm(attend_or_not ~ stats_joy + temperature*sun_joy, 
              data = stats_lectures, family = binomial())

model4 <- glm(attend_or_not ~ stats_joy + temperature + 
                sun_joy + temperature:sun_joy, 
              data = stats_lectures, family = binomial())
```

#### Subtask 2

```{r}
anova(model1, model2, model3, model4, test = "Chisq")
```

It is model 4 that fits the data best.

a)
```{r}
summary(model4)
```

b) We could deduct the hypothesis: 
   1. The more students enjoy statistics, the higher is the probability, that they visit statistics classes.
   2. Is there an interaction between temperature and enjoyment of sun.
      Might it be, that higher temperatures lead to less probability to visit classes 
      for students, that enjoy sun?

#### Subtask 3

```{r}
# Odds Ratios
exp(model4$coefficients)

# Confidence intervals, don't forget exp()
exp(confint(model4))
```

#### Subtask 4

a) See the text above

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

## 5 Render (knit)

Render your document using `Strg` + `Shift` + `K` (Windows) oder `Cmd` + `Shift` + `K` (Mac).
After that, you shoud see a nice looking version of your document.
I that works out, congratulations!
Your code could be rendered without any errors.
If not, don't be frustrated. 
Fix your errors. We will help you in our exercise class.


```{r eval=FALSE, include=FALSE}
#-----------------------#
#--- DATENSIMULATION ---#
#-----------------------#

# Hier werden die Daten für diesen Übungszettel simuliert.

#-----------------------#
#---      PAKETE     ---#
#-----------------------#
library(tidyverse)
library(VGAM) # für Funktion logit() zum Umwandeln von logit in Wahrscheinlichkeiten.

#---------------------------#
#--- VARIABLEN ERSTELLEN ---#
#---------------------------#
set.seed(12) # seed für Reproduzierbarkeit

temperature <- rnorm(m = 5, sd = 3, n = 400) %>% # Temperature
  scale() %>% # Z standardized for better calculation of our outcome
  as.double() # scale() creates a matrix, we convert that back to numbers

sun_joy <- rnorm(m = 5, sd = 3, n = 400) %>% # Joy with sun
  scale() %>% # Z standardized for better calculation of our outcome
  as.double() # scale() creates a matrix, we convert that back to numbers

# 5 Items of a questionnaire
stats_joy1 <- rnorm(m = 5, sd = 3, n = 400) %>% round(0)
stats_joy2 <- rnorm(m = 5, sd = 3, n = 400) %>% round(0)
stats_joy3 <- rnorm(m = 5, sd = 3, n = 400) %>% round(0)
stats_joy4 <- rnorm(m = 5, sd = 3, n = 400) %>% round(0)
stats_joy5 <- rnorm(m = 5, sd = 3, n = 400) %>% round(0)

# Join all in a dataframe or tibble
stats_data <- tibble(stats_joy1, stats_joy2, 
                     stats_joy3, stats_joy4, stats_joy5)

# We calculate the mean of the 5 items for all observations
stats_joy <- rowMeans(stats_data) 


#---------------------------#
#---  Creation of the OUTCOME  ---#
#---------------------------#
attending_logit = 1 + 2*stats_joy - 0.1*temperature - 0.1*sun_joy -
  1*temperature*sun_joy + rnorm(0, sd = 2, n = 1)

#---------------------------#
#--- Creation of dv      ---#
#---------------------------#

# Conversion of logits to probabilities
attending_prob = logit(attending_logit, inverse = TRUE)

# We use the probabilities to calculate the 0 or 1 values
attend_or_not <- rbinom(400, size = 1, prob = attending_prob)

# We create a factor
attend_or_not <- factor(attend_or_not, levels = c(0, 1), labels = c("course not attended", "course attended"))

#---------------------------#
#---     Store date      ---#
#---------------------------#

# Create a tibble of all data
stats_lectures <- tibble(attend_or_not, temperature, sun_joy, stats_joy1, stats_joy2, stats_joy3, stats_joy4, stats_joy5)

# ... and write that to a file
write_csv(stats_lectures, "data/stats_lectures.csv")
```


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


Version: `r format(Sys.time(), '%d %B, %Y %H:%M')`