---
lesson: Statistik 5
knitr:
  opts_chunk: 
    collapse: false
musterloesung: true
---


# Statistik 5: Übung

**Datensatz *polis.csv***

Der Datensatz polis.csv beschreibt für 19 Inseln im Golf von Kalifornien, ob 
Eidechsen der Gattung Uta vorkommen (presence/absence: PA) in Abhängigkeit von 
der Form der Inseln (Verhältnis Umfang zu Fläche: RATIO)

Bitte prüft mit einer logistischen Regression, ob und ggf. wie die Inselform die Präsenz der Eidechsen beinflusst


::::{.content-hidden unless-meta="musterloesung"}
:::{.callout-note}

[Demoscript herunterladen (.R)](../purl/Statistik5_Uebung.R){.dld}

[Demoscript herunterladen (.qmd)](../purl/Statistik5_Uebung.qmd){.dld}


- [Lösungstext als Download](Statistik5_Loesung.pdf)



## Musterlösung Übung 5

```{r}
#| message: false
library("readr")
library("ggplot2")
library("dplyr")
```

```{r}
#| message: false
polis <- read_delim("datasets/stat/polis.csv", delim = ";") |>
  mutate(across(where(is.character), as.factor))
                                     
str(polis)
summary(polis)
```

Man erkennt, dass polis 19 Beobachtungen von drei Parametern enthält, wobei 
ISLAND ein Faktor mit den Inselnamen ist, während RATIO metrisch ist und PA 
nur 0 oder 1 enthält. Prädiktorvariable ist RATIO, abhängige Variable PA, 
mithin ist das korrekte statistische Verfahren eine logistische Regression (GLM).


```{r}
# Definition des logistischen Modells
glm_1 <- glm(PA ~ RATIO, family = "binomial", data = polis)
summary(glm_1)
```

Modell ist signifikant (p-Wert in Zeile RATIO ist < 0.05). Jetzt müssen wir noch prüfen, ob es auch valide ist:

```{r}
# Modelldiagnostik für das gewählte Modell (wenn nicht signifikant, dann OK)
1 - pchisq(glm_1$deviance, glm_1$df.resid)

library(DHARMa)
set.seed(123)
simulateResiduals(fittedModel = glm_1, plot = TRUE, n = 1000)
```

Ist ok
Jetzt brauchen wir noch die Modellgüte (Pseudo-R2):

```{r}
# Modellgüte (pseudo-R²)
library("performance")
r2(glm_1)
```

Um unser Modell zu interpretieren müssen wir noch in Betracht ziehen, dass 
wir nicht die Vorkommenswahrscheinlichkeit selbst, sondern logit (Vorkommenswahrscheinlichkeit) modelliert haben. Unser Ergebnis (die beiden Parameterschätzungen von oben, also 3.6061 und -0.2196) muss also zunächst in 
etwas Interpretierbares übersetzt werden:

```{r}
# Steilheit der Beziehung in Modellen mit nur einem Parameter
exp(glm_1$coef[2])
```

< 1, d. h. Vorkommenswahrscheinlichkeit sinkt mit zunehmender Isolation.

```{r}
# LD50 für 1-Parameter-Modelle (hier also x-Werte, bei der 50% der Inseln 
# besiedelt sind)
library(MASS)
dose.p(glm_1, p = 0.5)
```

Am besten stellen wir auch unsere Funktionsgleichung dar. Dazu müssen wir das „Rohergebnis“ (mit P = Vorkommenswahrscheinlichkeit)

ln (P/ (1- P)) = 3.606 – 0.220 RATIO

so umformen, dass wir links nur P stehen haben:

P = exp(3.606 – 0.220 RATIO) / (1 + exp(3.606 – 0.220 RATIO))

Das ist also unsere vorhergesagte Regressionsfunktion, die wir in einem letzten Schritt auch noch visualisieren können (und sollten):

```{r}
# Ergebnisplots
ggplot(data = polis, aes(x = RATIO, y = PA)) +
  geom_point() +
  labs(x = "Umfang-Flächen-Verhältnis", y = "Vorkommenswahrscheinlichkeit") +
  geom_smooth(method = "glm", method.args = list(family = "binomial")) +
  theme_classic()
```

:::
:::: 
