---
lesson: Statistik 4
knitr:
  opts_chunk: 
    collapse: false
musterloesung: false
---

# Statistik 4: Übung

- Lesen Sie den Datensatz *steprasen_ukraine.csv* in R  ein. Dieser enthält Pflanzenartenzahlen (Species.richness) von 199
10 m² grossen Plots (Vegetationsaufnahmen) von Steppenrasen in der Ukraine sowie zahlreiche Umweltvariablen, deren Bedeutung und Einheiten im Kopf der ExcelTabelle angegeben sind.
- **Ermitteln Sie ein minimal adäquates Modell, das den Artenreichtum in den Plots durch die Umweltvariablen erklärt.**
- Bitte erklären und begründen Sie die einzelnen Schritte, die Sie unternehmen, um zu diesem Ergebnis zu kommen. Dazu erstellen Sie bitte ein Word-Dokument, in das Sie Schritt für Schritt den verwendeten R-Code, die dazu gehörigen Ausgaben von R, Ihre Interpretation derselben und die sich ergebenden Schlussfolgerungen für das weitere Vorgehen dokumentieren.
- Dieser Ablauf sollte insbesondere beinhalten:
  - Überprüfen der Datenstruktur nach dem Einlesen: welches sind die abhängige(n)
    und welches die unabängige(n) Variablen, sind alle Variablen für die Analyse
    geeignet?
  - Explorative Datenanalyse, um zu sehen, ob die abhängige Variable in der
    vorliegenden Form für die Analyse geeignet ist
  - Definition eines globalen Modelles und dessen Reduktion zu einem minimal
    adäquaten Modell
  - Durchführen der Modelldiagnostik für dieses
  - Generieren aller Zahlen, Statistiken und Tabellen, die für eine wiss.
    Ergebnisdarstellung benötigt werden
  - Formulieren Sie abschliessend einen Methoden- und Ergebnisteil (ggf. incl.
    adäquaten Abbildungen) zu dieser Untersuchung in der Form einer
    wissenschaftlichen Arbeit (ausformulierte schriftliche Zusammenfassung, mit je
    einem Absatz von ca. 60-100 Worten, resp. 3-8 Sätzen für den Methoden- und
    Ergebnisteil). D. h. alle wichtigen Informationen sollten enthalten sein, unnötige
    Redundanz dagegen vermieden werden.
  - **Zu erstellen sind (a) Ein lauffähiges R-Skript; (b) begründeter Lösungsweg
    (Kombination aus R-Code, R Output und dessen Interpretation) und (c)
    ausformulierter Methoden- und Ergebnisteil (für eine wiss. Arbeit).**


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

:::{.callout-note}

## Lösung Übung 4

[Download dieses Demoscript ](../purl/Statistik4_Uebung.R) (rechtsklick → speichern unter)


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

```{r}
library(pacman)
p_load("tidyverse")

ukr <- read_delim("datasets/stat/steprasen_ukraine.csv", delim = ";")
str(ukr)
```
Man erkennt, dass alle Spalten bis auf die erste mit der Plot ID numerisch (num) sind. Wir können nun die Variable "PlotID" als rownames setzten

```{r}
# Setze erste Spalte (colum) als rownames
ukr <- ukr |>
  column_to_rownames(var = "PlotID")
```
Nun steht die abhängige Variable in Spalte 1 und die Prediktorvariablen in den Spalten  2 bis 19 stehen.


```{r}
summary(ukr)
```

Da es sich bei Artenzahlen um Zähldaten handelt, müsste man theoretisch ein glm mit Poisson-Verteilung rechnen; bei einem Mittelwert, der hinreichend von Null verschieden ist (hier: ca. 40), ist eine Poisson-Verteilung aber praktisch nicht von einer Normalverteilung zu unterscheiden und wir können uns den Aufwand auch sparen).


```{r}
#| eval: false

cor <- cor(ukr[, 2:19])
cor[abs(cor)<0.7] <- 0
cor
```
Die Korrelationsanalyse dient dazu, zu entscheiden, ob die Prädiktorvariablen hinreichend voneinander unabhängig sind, um alle in das globale Modell hinein zu nehmen. Bei Pearson’s Korrelationskoeffizienten r, die betragsmässig grösser als 0.7 sind, würde es problematisch. Alternativ hätten wir auch den VIF (Variance Inflation Factor) als Kriterium für den möglichen Ausschluss von Variablen aus dem globalen Modell nehmen können. 

Wenn man auf cor nun doppel-klickt und es in einem separaten Fenster öffnet, sieht man, wo es problematische Korrelationen zwischen Variablenpaaren gibt. Es sind dies Altitude vs. Temperature und N.total vs. C.org. Wir müssen aus jedem dieser Paare jetzt eine Variable rauswerfen, am besten jene, die weniger gut interpretierbar ist. Ich entscheide mich dafür Temperature statt Altitude (weil das der direktere ökologische Wirkfaktor ist) und C.org statt N.total zu behalten (weil es in der Literatur mehr Daten zum Humusgehalt als zum N-Gehalt gibt, damit eine bessere Vergleichbarkeit erzielt wird). Die Aussagen, die wir für die beibehaltene Variable erzielen, stehen aber +/- auch für die entfernte.

Das Problem ist aber, dass wir immer noch 16 Variablen haben, was einen sehr leistungsfähigen Rechner oder sehr lange Rechenzeit erfordern würde. Wir sollten also unter 15 Variablen kommen. Wir könnten uns jetzt überlegen, welche uns ökologisch am wichtigsten sind, oder ein noch strengeres Kriterium bei r verwenden, etwa 0.6


```{r}
cor <- cor(ukr[, c(2:19)])
cor[abs(cor) < 0.6] <- 0
cor
```

```{r}
#| include: false

# write.table(cor, file = "stat/Correlation.csv", sep = ";", dec = ".", col.names = NA)
```

Entsprechend „werfen“ wir auch noch die folgenden Variablen „raus“: Temperature.range (positiv mit Temperature), Precipitation (negativ mit Temperature) sowie Conductivity (positiv mit pH).

Nun können wir das globale Modell definieren, indem wir alle verbleibenden Variablen aufnehmen, das sind 13. (Wenn das nicht eh schon so viele wären, dass es uns an die Grenze der Rechenleistung bringt, hätten wir auch noch darüber nachdenken können, einzelne quadratische Terme oder Interaktionsterme zu berücksichtigen).

```{r}
global_model <- lm(Species_richness ~ Inclination + Heat_index + Microrelief + Grazing_intensity +
    Litter + Stones_and_rocks + Gravel + Fine_soil + pH + CaCO3 + C_org + CN_ratio + Temperature, data = ukr)
```

Nun gibt es im Prinzip zwei Möglichkeiten, vom globalen (vollen) Modell zu einem minimal adäquaten Modell zu kommen. (1) Der Ansatz der „frequentist statistic“, in dem man aus dem vollen Modell so lange schrittweise Variablen entfernt, bis nur noch signifikante Variablen verbleiben. (2) Den informationstheoretischen Ansatz, bei dem alle denkbaren Modelle berechnet und verglichen werden (also alle möglichen Kombinationen von 13,12,…, 1, 0 Parametern). Diese Lösung stelle ich im Folgenden vor:

```{r}
# Multimodel inference
p_load("MuMIn")

options(na.action = "na.fail")
allmodels <- dredge(global_model)
```

```{r}
#| eval: false

allmodels
```

Jetzt bekommen wir die besten der insgesamt 8192 möglichen Modelle gelistet mit ihren Parameterschätzungen und ihrem AICc.

Das beste Modell umfasst 5 Parameter (CaCO3, CN.ratio, Grazing.intensity. Heat.index, Litter). Allerdings ist das nächstbeste Modell (mit 6 Parametern) nur wenig schlechter (delta AICc = 0.71), was sich in fast gleichen (und zudem sehr niedrigen) Akaike weights bemerkbar macht. Nach dem Verständnis des Information theoretician approach, sollte man in einer solchen Situation nicht das eine „beste“ Modell benennen, sondern eine Aussage über die Gruppe der insgesamt brauchbaren Modelle treffen. Hierzu kann man (a) Importance der Parameter über alle Modelle hinweg berechnen (= Summe der Akaike weights aller Modelle, die den betreffenden Parameter enthalten) und/oder (b) ein nach Akaike weights gemitteltes Modell berechnen.

```{r}
# Importance values der Variablen
sw(allmodels)
```

Demnach ist Heat.index die wichtigste Variable (in 100% aller relevanten Modelle), während ferner Litter, CaCO3, CN_ratio und Grazing_intensity in mehr als 50% der relevanten Modelle enthalten sind.

```{r}
#| eval: false

# Modelaveraging (Achtung: dauert mit 13 Variablen einige Minuten)
avgmodel <- model.avg(allmodels)
summary(avgmodel)
```

Aus dem gemittelten Modell können wir die Richtung der Beziehung (positiv oder negativ) und ggf. die Effektgrössen (wie verändert sich die Artenzahl, wenn die Prädiktorvariable um eine Einheit zunimmt?) ermitteln.

```{r}
# Modelldiagnostik nicht vergessen
par(mfrow = c(2, 2))
plot(global_model)
plot(lm(Species_richness ~ Heat_index + Litter + CaCO3 + CN_ratio + Grazing_intensity, data = ukr))

```

Wie immer kommt am Ende die Modelldiagnostik. Wir können uns entweder das globale Modell oder das Modell mit den 5 Variablen mit importance > 50% anschauen. Das Bild sieht fast identisch aus und zeigt keinerlei problematische Abweichungen, d. h. links oben weder ein Keil, noch eine Banane, rechts oben eine nahezu perfekte Gerade. 

:::
:::: 