---
lesson: Statistik4
knitr:
  opts_chunk: 
    collapse: false 
---

# Statistik 4: Demo

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

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

-   Datensatz *loyn.csv*

```{r}
library(readr)
library(ggplot2)
```


## Multiple lineare Regression

```{r}
loyn <- read_delim("datasets/stat/loyn.csv", delim = ";")

str(loyn)
summary(loyn)
```

### Korrelation zwischen den Prädiktoren

```{r}
# Wir setzen die Schwelle bei |0.7|
# Korrelationen rechnen details siehe: "?cor"
cor(loyn[, 2:7]) 

# oder mit Namen der columns resp. variablen
cor1 <- 
  loyn |>
  subset(select = AREA:ALT) |>
  cor()

# Korrelationen Visualisieren
library("corrplot")
corrplot.mixed(cor1, lower = "ellipse", upper = "number", order = "AOE")

# Nur Darstellung von variablen über Schwelenwert von 0.7
cor1[abs(cor1)<0.7] <- 0
cor1
```

-> Keine Korrelation ist \>\|0.7\|, so können wir alle Prädiktoren "behalten". 
Aber es gilt zu beachten, dass GRAZE ziemlich stark \|\>0.6\| mit AGE korreliert ist

```{r}
# Volles Modell definieren
names(loyn)
lm_1 <- lm(ABUND ~ AGE + AREA + DIST + LDIST + GRAZE + ALT, data = loyn)
```

```{r}
# Andere Variante, um korrelierte Prädiktoren zu finden 
# (üblicher Schwellenwert VIF = 5)
library("car")
vif(lm_1)
```

```{r}
# Modeldiagnostik
par(mfrow = c(2, 2))
plot(lm_1)
```

-> Plot sieht zwar OK aus, aber mit 6 Prädiktoren für \|\<60\| Beobachtungen ist das Modell wohl "overfitted"


### Modellselektion

```{r}
summary(lm_1)

drop1(lm_1, test = "F")
# Prädiktor mit grösstem p-Wert entfernen
lm_2 <- lm(ABUND ~ AGE + AREA + DIST  + GRAZE + ALT, data = loyn)
# oder
lm_2 <- update(lm_1, ~ . - LDIST) 

# Oben beschriebener Schritt wiederholen bis nur noch signifikante Prädiktoren 
# im Modell sind
drop1(lm_2, test = "F") 
lm_3 <- update(lm_2, ~ . - DIST)

drop1(lm_3, test = "F") 
lm_4 <- update(lm_3, ~ . - ALT)

drop1(lm_4, test = "F") 
lm_5 <- update(lm_4, ~ . - AGE)

drop1(lm_5, test = "F")

summary(lm_5) 

par(mfrow = c(2, 2))
plot(lm_5)
```

-> das minimal adäquate Modell enthält noch zwei Prädiktoren (AREA; GRAZE) und dessen Residualplots sehen ok aus.

### Hierarchical partitioning 

Wir können auch schauen wie bedeutsam die einzelnen Variablen sind:

```{r}
library("relaimpo")

# Berechnen
metrics <- calc.relimp(lm_1, type = c("lmg", "first"))
cbind(I = metrics$lmg, 
      J = metrics$first - metrics$lmg, 
      Total = metrics$first)
```

-> auch hier haben AREA und GRAZE die höchsten Werte (und an dritter Stelle AGE, 
der mit GRAZE am stärksten korreliert ist)

### Plot partielle regressionen

```{r}
# Beispiel GRAZE
lm_abund <- lm(ABUND ~ AREA, data = loyn)
lm_graze <- lm(GRAZE ~ AREA, data = loyn)

abundance_resid <- resid(lm_abund)
graze_resid <- resid(lm_graze)

ggplot(data = NULL, aes(x = graze_resid, y = abundance_resid)) +
  geom_point() +
  geom_smooth(method = "lm", se = FALSE, color = "blue") +
  labs(x = "Graze | others", y = "Abund | others") +
  theme_minimal()

# Einfacher geht es mit der function avPlots (package "car"). Nachteil ist, dass mit der funktion anders als mit der Methode oben, keine quadratische prädiktoren dargestellt werden können 
par(mfrow = c(1, 1))
avPlots(lm_5, ~GRAZE, ask = F)

#Für alle Prädiktoren im Modell
avPlots(lm_5, layout = c(1, 2) )
```

## Multimodel inference

```{r}
library("MuMIn")

global_model <- lm(ABUND ~ AGE + AREA + DIST + LDIST + GRAZE + ALT, data = loyn)

options(na.action = "na.fail")
allmodels <- dredge(global_model)
allmodels
# Wir haben mehre Modelle mit einem delta AICc <2, das heisst wir haben nicht 
# ein eindeutig bestes Modell (welches wir mit der funktion "get.models" 
# selektieren könnten)

# Variable importance
sw(allmodels)
```
-> Auch mit dieser Sichtweise sind AREA und GRAZE die wichtigste Prädiktoren

```{r}
# Model averaging
avgmodel <- model.avg(allmodels)
summary(avgmodel)

# Nur Estimates
summary(avgmodel)$coefficients

# Confindence intervals
confint(avgmodel)
```
