---
lesson: Statistik6
knitr:
  opts_chunk: 
    collapse: true
---

# Statistik 6: Demo

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


-   Datensatz *rm_plants.csv*
-   Datensatz *DeerEcervi.csv*

```{r}
# Packete laden
library("pacman")
p_load("tidyverse", "sjPlot", "car", "lme4", "glmmTMB", "performance", "DHARMa")
```

## Repeated measurement

Beispiel mit Daten von Lepš & Šmilauer (2020)


```{r}
# Daten laden
plantf <- read_delim("datasets/stat/rm_plants.csv", delim = ";") |>
  mutate(across(where(is.character), as.factor))

str(plantf)
summary(plantf)
table(plantf$Treatment, plantf$PlantID)

# Mit aov
pf_eaov <- aov(PlantHeight ~ Treatment * Time  + Error(PlantID), data = plantf)
summary(pf_eaov)

# Als lmm 
pf_lmm <- glmmTMB(PlantHeight ~ Treatment * Time + (1 | PlantID), 
                 family = gaussian, 
                 REML = TRUE,
                 data = plantf)

Anova(pf_lmm)
summary(pf_lmm)

# Darstellung (Interaktions Plot)
plot_model(pf_lmm, 
           type = "pred", pred.type = "re",
           terms = c("Time", "Treatment") ) +
  theme_classic()
```


## Split-plot Design

```{r}
# Import Data
glex <- read_delim("datasets/stat/Riesch_et_al_2020_grassland.csv", delim = ";") |>
  mutate(across(where(is.character), as.factor))

str(glex)
glex$year <- as.factor(glex$year)
summary(glex)

table(glex$site_code, glex$treatment, glex$plot_type, glex$year)
# balanced Design
```

## 2. Fit LMM with random interecept

```{r}
# REML = TRUE : (Restricted maximum likelihood) v.s 
# REML = FALSE: Maximum likelihood (ML) 
# Bei REML sind die estimates genauer, aber REML sollte nicht für likelihood 
# ratio test (drop1) benutzt werden
# Dies ist nur relevant für Gaussian mixed models (LMM) nicht für GLMMs

# 2.1 Model fitten
lmm_1 <- glmmTMB(SR ~ year * treatment * plot_type + 
               (1| site_code/treatment/plot_type), 
               family = gaussian, 
               REML = FALSE,
               data = glex)

Anova(lmm_1)

# Model optimierung
drop1(lmm_1, test = "Chi")

lmm_2 <- update(lmm_1,~. -year:treatment:plot_type)
drop1(lmm_2, test = "Chi")

lmm_3 <- update(lmm_2,~. -treatment:plot_type)
drop1(lmm_3, test = "Chi")

# Refit with REML
lmm_4 <- update(lmm_3, REML = TRUE)

# Resultat
Anova(lmm_4)
summary(lmm_4)
r2(lmm_4)

# Model validierung mit dem package "DHARMa"

# Detailierte Model validation mit dem package 
set.seed(123)
simulationOutput <- simulateResiduals(fittedModel = lmm_4, plot = TRUE, n = 1000)

# Plot residuals vs, covariates des models
plotResiduals(simulationOutput, form = glex$year)
plotResiduals(simulationOutput, form = glex$treatment)
plotResiduals(simulationOutput, form = glex$plot_type)

# Darstellung
plot_model(lmm_4, 
           type = "pred", pred.type = "re", 
           terms = c("year", "treatment", "plot_type") ) +
  theme_classic()
```

## Random slope & random intercept

```{r}
# Daten laden
data(sleepstudy)
?sleepstudy

# Daten ohne Trainingsphase
sleepstudy_2 <- subset(sleepstudy, Days>=2)
str(sleepstudy_2)
summary(sleepstudy_2)

table(sleepstudy_2$Subject)

# Visualisierung
ggplot(sleepstudy_2, aes(y = Reaction, x = Days)) +
  geom_point() +
  xlab("Average reaction time (ms)") +
  ylab("Number of days of sleep deprivation") +
  geom_smooth(method = "lm", formula = 'y ~ x', se = F, fullrange = T) +
  theme_classic() +
  facet_wrap(~Subject)
```

Wie man in der Visualisierung sehen kann, unterscheiden sich nicht nur die Intercepts (Reaktionszeit ohne Schlafmangel), sondern der Schlafmangel scheint sich auch unterschiedlich stark auf die Reaktionszeit der Probanden auszuwirken. In diesem Fall ist es daher sinnvoll, nicht nur einen Random Intercept, sondern auch einen Random Slope zu fitten.

```{r}
# Fit model
lmm_1 <- glmmTMB(Reaction ~ Days + (Days | Subject),
                 family = gaussian,
                 REML = TRUE,
                 data = sleepstudy_2)

summary(lmm_1)

# Modelvalidierung

# Model validation mit DHARMa
set.seed(123)
simulationOutput <- simulateResiduals(fittedModel = lmm_1, plot = TRUE, n = 1000)

plot_model(lmm_1, 
           type = "pred", pred.type = "re",
           show.data = TRUE) +
  theme_classic()
```

## GLMM

```{r}
# Daten laden und für GLMM aufbereiten
DeerEcervi <- read_delim("datasets/stat/DeerEcervi.csv", delim = ";") |>
  mutate(across(where(is.character), as.factor))

# Daten anschauen
str(DeerEcervi)
head(DeerEcervi)
summary(DeerEcervi)
table(DeerEcervi$Farm)
```


```{r}
# Kontinuierliche variable Hischlänge standardisieren
Std <- function(x) { (x - mean(x)) / sd(x)}
DeerEcervi$Lengt_std <-  Std(DeerEcervi$Length)
# oder
DeerEcervi$Length_std <- as.vector( scale(DeerEcervi$Length, center = TRUE))

# Model fitten
glmm_1 <- glmmTMB(Ecervi ~ Length_std * Sex + (1 | Farm), 
                  family = binomial, 
                  data = DeerEcervi)

car::Anova(glmm_1)
drop1(glmm_1, test = "Chi")

summary(glmm_1)
r2(glmm_1)
```

```{r}
# Model validation

# Test overdisepersion
check_overdispersion(glmm_1)

# Model validation in detail using DHARMa
set.seed(123)
simulationOutput <- simulateResiduals(fittedModel = glmm_1, plot = TRUE, n = 1000)

# Plot the scaled quantile residuals versus fitted values.
plotResiduals(simulationOutput, form = DeerEcervi$Sex)
plotResiduals(simulationOutput, form = DeerEcervi$Length_std) 
```

```{r}
#Visualize result
plot_model(glmm_1, 
           type = "pred", pred.type = "re", 
           terms = c("Length_std[all]", "Sex"),  
           show.data = TRUE) +
  theme_classic()
```
