Data_modelling.Rmd

---
title: "Data_modelling"
output: html_document
date: '2023-03-17'
---

```{r setup, include=FALSE}
renv::restore()
library(tidyverse)
set.seed(42)
library(mgcv)
library(lubridate)
library(timetk)
library(RcppRoll)
library(dLagM)
library(dlnm)
library(zoo)
theme_set(theme_gray())
library(knitr)
library(kableExtra)
library(sf)
library(tmap)
library(gridExtra)
library(lubridate)
library(timetk)
library(RcppRoll)
library(dLagM)
library(dlnm)
library(zoo)
theme_set(theme_gray())
library(js)
library(sjPlot)
library(verification)
library(ROCR)
library(pROC)
# Read data
source("data_read.R")
```

# GLM Approaches 
## ZV: LowLevel (Binary)
## temp
``` {r}
# summer

y <- hydro_summer_10304$lowlevel
x <- hydro_summer_10304$airtmp
cor(x, y, method = "pearson")
glm(y~x, data = data.frame(x, y), family = binomial())
event_summer_10304_tmp <- ggplot(data = data.frame(x, y), aes(x = x , y = y)) +
  geom_point() +
  geom_smooth(method = "glm", method.args = list(family = "binomial"), se = TURE) +
  labs(title = "Summer 10304", x = "Temp", y = "LowLevel")
event_summer_10304_tmp0 <- ggplot(data.frame(y, x) %>% filter(y == 0), aes(x = x, y = y)) +
  geom_boxplot() +
  labs(title = "LowLevel = 0", x = "tmp", y = "Low Level")
event_summer_10304_tmp1 <- ggplot(data.frame(y, x) %>% filter(y == 1), aes(x = x, y = y)) +
  geom_boxplot() +
  labs(title = "LowLevel = 1", x = "tmp", y = "Low Level")
grid_event_summer_10304_tmp <- grid.arrange(event_summer_10304_tmp, event_summer_10304_tmp0, event_summer_10304_tmp1, ncol = 1)


y <- hydro_summer_11502$lowlevel
x <- hydro_summer_11502$airtmp
cor(x, y, method = "pearson")
glm(y~x, data = data.frame(x, y), family = binomial())
event_summer_11502_tmp <- ggplot(data = data.frame(x, y), aes(x = x , y = y)) +
  geom_point() +
  geom_smooth(method = "glm", method.args = list(family = "binomial"), se = TRUE) +
  labs(title = "Summer 11502", x = "Temp", y = "LowLevel")
event_summer_11502_tmp0 <- ggplot(data.frame(y, x) %>% filter(y == 0), aes(x = x, y = y)) +
  geom_boxplot() +
  labs(title = "LowLevel = 0", x = "tmp", y = "Low Level")
event_summer_11502_tmp1 <- ggplot(data.frame(y, x) %>% filter(y == 1), aes(x = x, y = y)) +
  geom_boxplot() +
  labs(title = "LowLevel = 1", x = "tmp", y = "Low Level")
grid_event_summer_11502_tmp <- grid.arrange(event_summer_11502_tmp, event_summer_11502_tmp0, event_summer_11502_tmp1, ncol = 1)


y <- hydro_summer_20203$lowlevel
x <- hydro_summer_20203$airtmp
cor(x, y, method = "pearson")
glm(y~x, data = data.frame(x, y), family = binomial())
event_summer_20203_tmp <- ggplot(data = data.frame(x, y), aes(x = x , y = y)) +
  geom_point() +
  geom_smooth(method = "glm", method.args = list(family = "binomial"), se = TRUE) +
  labs(title = "Summer 20203", x = "Temp", y = "LowLevel")
event_summer_20203_tmp0 <- ggplot(data.frame(y, x) %>% filter(y == 0), aes(x = x, y = y)) +
  geom_boxplot() +
  labs(title = "LowLevel = 0", x = "tmp", y = "Low Level")
event_summer_20203_tmp1 <- ggplot(data.frame(y, x) %>% filter(y == 1), aes(x = x, y = y)) +
  geom_boxplot() +
  labs(title = "LowLevel = 1", x = "tmp", y = "Low Level")
grid_event_summer_20203_tmp <- grid.arrange(event_summer_20203_tmp, event_summer_20203_tmp0, event_summer_20203_tmp1, ncol = 1)

grid.arrange(grid_event_summer_10304_tmp, grid_event_summer_11502_tmp, grid_event_summer_20203_tmp, ncol = 2)

# winter
y <- hydro_winter_10304$lowlevel
x <- hydro_winter_10304$airtmp
cor(x, y, method = "pearson")
glm(y~x, data = data.frame(x, y), family = binomial())
event_winter_10304_tmp <- ggplot(data = data.frame(x, y), aes(x = x , y = y)) +
  geom_point() +
  geom_smooth(method = "glm", method.args = list(family = "binomial"), se = TRUE) +
  labs(title = "Winter 10304", x = "Temp", y = "LowLevel")
event_winter_10304_tmp0 <- ggplot(data.frame(y, x) %>% filter(y == 0), aes(x = x, y = y)) +
  geom_boxplot() +
  labs(title = "LowLevel = 0", x = "tmp", y = "Low Level")
event_winter_10304_tmp1 <- ggplot(data.frame(y, x) %>% filter(y == 1), aes(x = x, y = y)) +
  geom_boxplot() +
  labs(title = "LowLevel = 1", x = "tmp", y = "Low Level")
grid_event_winter_10304_tmp <- grid.arrange(event_winter_10304_tmp, event_winter_10304_tmp0, event_winter_10304_tmp1, ncol = 1)


y <- hydro_winter_11502$lowlevel
x <- hydro_winter_11502$airtmp
cor(x, y, method = "pearson")
glm(y~x, data = data.frame(x, y), family = binomial())
event_winter_11502_tmp <- ggplot(data = data.frame(x, y), aes(x = x , y = y)) +
  geom_point() +
  geom_smooth(method = "glm", method.args = list(family = "binomial"), se = TRUE) +
  labs(title = "Winter 11502", x = "Temp", y = "LowLevel")
event_winter_11502_tmp0 <- ggplot(data.frame(y, x) %>% filter(y == 0), aes(x = x, y = y)) +
  geom_boxplot() +
  labs(title = "LowLevel = 0", x = "tmp", y = "Low Level")
event_winter_11502_tmp1 <- ggplot(data.frame(y, x) %>% filter(y == 1), aes(x = x, y = y)) +
  geom_boxplot() +
  labs(title = "LowLevel = 1", x = "tmp", y = "Low Level")
grid_event_winter_11502_tmp <- grid.arrange(event_winter_11502_tmp, event_winter_11502_tmp0, event_winter_11502_tmp1, ncol = 1)


y <- hydro_winter_20203$lowlevel
x <- hydro_winter_20203$airtmp
cor(x, y, method = "pearson")
glm(y~x, data = data.frame(x, y), family = binomial())
event_winter_20203_tmp <- ggplot(data = data.frame(x, y), aes(x = x , y = y)) +
  geom_point() +
  geom_smooth(method = "glm", method.args = list(family = "binomial"), se = TRUE) +
  labs(title = "Winter 20203", x = "Temp", y = "LowLevel")
event_winter_20203_tmp0 <- ggplot(data.frame(y, x) %>% filter(y == 0), aes(x = x, y = y)) +
  geom_boxplot() +
  labs(title = "LowLevel = 0", x = "tmp", y = "Low Level")
event_winter_20203_tmp1 <- ggplot(data.frame(y, x) %>% filter(y == 1), aes(x = x, y = y)) +
  geom_boxplot() +
  labs(title = "LowLevel = 1", x = "tmp", y = "Low Level")
grid_event_winter_20203_tmp <- grid.arrange(event_winter_20203_tmp, event_winter_20203_tmp0, event_winter_20203_tmp1, ncol = 1)

grid.arrange(grid_event_summer_10304_tmp, grid_event_summer_11502_tmp, grid_event_summer_20203_tmp, ncol = 2)
grid.arrange(grid_event_summer_10304_tmp, grid_event_summer_11502_tmp, grid_event_summer_20203_tmp, grid_event_summer_10304_tmp, grid_event_summer_11502_tmp, grid_event_summer_20203_tmp, ncol = 3)

```

## precip
``` {r}
# summer

y <- hydro_summer_10304$lowlevel
x <- hydro_summer_10304$precip
cor(x, y, method = "pearson")
glm(y~x, data = data.frame(x, y), family = binomial())
event_summer_10304_precip <- ggplot(data = data.frame(x, y), aes(x = x , y = y)) +
  geom_point() +
  geom_smooth(method = "glm", method.args = list(family = "binomial"), se = TRUE) +
  labs(title = "Summer 10304", x = "precip", y = "LowLevel")
event_summer_10304_precip0 <- ggplot(data.frame(y, x) %>% filter(y == 0), aes(x = x, y = y)) +
  geom_boxplot() +
  labs(title = "LowLevel = 0", x = "tmp", y = "Low Level")
event_summer_10304_precip1 <- ggplot(data.frame(y, x) %>% filter(y == 1), aes(x = x, y = y)) +
  geom_boxplot() +
  labs(title = "LowLevel = 1", x = "tmp", y = "Low Level")
grid_event_summer_10304_precip <- grid.arrange(event_summer_10304_precip, event_summer_10304_precip0, event_summer_10304_precip1, ncol = 1)



y <- hydro_summer_11502$lowlevel
x <- hydro_summer_11502$precip
cor(x, y, method = "pearson")
glm(y~x, data = data.frame(x, y), family = binomial())
event_summer_11502_precip <- ggplot(data = data.frame(x, y), aes(x = x , y = y)) +
  geom_point() +
  geom_smooth(method = "glm", method.args = list(family = "binomial"), se = TRUE) +
  labs(title = "Summer 11502", x = "precip", y = "LowLevel")
event_summer_11502_precip0 <- ggplot(data.frame(y, x) %>% filter(y == 0), aes(x = x, y = y)) +
  geom_boxplot() +
  labs(title = "LowLevel = 0", x = "tmp", y = "Low Level")
event_summer_11502_precip1 <- ggplot(data.frame(y, x) %>% filter(y == 1), aes(x = x, y = y)) +
  geom_boxplot() +
  labs(title = "LowLevel = 1", x = "tmp", y = "Low Level")
grid_event_summer_11502_precip <- grid.arrange(event_summer_11502_precip, event_summer_11502_precip0, event_summer_11502_precip1, ncol = 1)



y <- hydro_summer_20203$lowlevel
x <- hydro_summer_20203$precip
cor(x, y, method = "pearson")
glm(y~x, data = data.frame(x, y), family = binomial())
event_summer_20203_precip <- ggplot(data = data.frame(x, y), aes(x = x , y = y)) +
  geom_point() +
  geom_smooth(method = "glm", method.args = list(family = "binomial"), se = TRUE) +
  labs(title = "Summer 20203", x = "precip", y = "LowLevel")
event_summer_20203_precip0 <- ggplot(data.frame(y, x) %>% filter(y == 0), aes(x = x, y = y)) +
  geom_boxplot() +
  labs(title = "LowLevel = 0", x = "tmp", y = "Low Level")
event_summer_20203_precip1 <- ggplot(data.frame(y, x) %>% filter(y == 1), aes(x = x, y = y)) +
  geom_boxplot() +
  labs(title = "LowLevel = 1", x = "tmp", y = "Low Level")
grid_event_summer_20203_precip <- grid.arrange(event_summer_20203_precip, event_summer_20203_precip0, event_summer_20203_precip1, ncol = 1)
grid.arrange(grid_event_summer_10304_precip, grid_event_summer_11502_precip, grid_event_summer_20203_precip, ncol = 2)

# winter
y <- hydro_winter_10304$lowlevel
x <- hydro_winter_10304$precip
cor(x, y, method = "pearson")
glm(y~x, data = data.frame(x, y), family = binomial())
event_winter_10304_precip <- ggplot(data = data.frame(x, y), aes(x = x , y = y)) +
  geom_point() +
  geom_smooth(method = "glm", method.args = list(family = "binomial"), se = TRUE) +
  labs(title = "Winter 10304", x = "precip", y = "LowLevel")
event_winter_10304_precip0 <- ggplot(data.frame(y, x) %>% filter(y == 0), aes(x = x, y = y)) +
  geom_boxplot() +
  labs(title = "LowLevel = 0", x = "tmp", y = "Low Level")
event_winter_10304_precip1 <- ggplot(data.frame(y, x) %>% filter(y == 1), aes(x = x, y = y)) +
  geom_boxplot() +
  labs(title = "LowLevel = 1", x = "tmp", y = "Low Level")
grid_event_winter_10304_precip <- grid.arrange(event_winter_10304_precip, event_winter_10304_precip0, event_winter_10304_precip1, ncol = 1)


y <- hydro_winter_11502$lowlevel
x <- hydro_winter_11502$precip
cor(x, y, method = "pearson")
glm(y~x, data = data.frame(x, y), family = binomial())
event_winter_11502_precip <- ggplot(data = data.frame(x, y), aes(x = x , y = y)) +
  geom_point() +
  geom_smooth(method = "glm", method.args = list(family = "binomial"), se = TRUE) +
  labs(title = "Winter 11502", x = "precip", y = "LowLevel")
event_winter_11502_precip0 <- ggplot(data.frame(y, x) %>% filter(y == 0), aes(x = x, y = y)) +
  geom_boxplot() +
  labs(title = "LowLevel = 0", x = "tmp", y = "Low Level")
event_winter_11502_precip1 <- ggplot(data.frame(y, x) %>% filter(y == 1), aes(x = x, y = y)) +
  geom_boxplot() +
  labs(title = "LowLevel = 1", x = "tmp", y = "Low Level")
grid_event_winter_11502_precip <- grid.arrange(event_winter_11502_precip, event_winter_11502_precip0, event_winter_11502_precip1, ncol = 1)


y <- hydro_winter_20203$lowlevel
x <- hydro_winter_20203$precip
cor(x, y, method = "pearson")
glm(y~x, data = data.frame(x, y), family = binomial())
event_winter_20203_precip <- ggplot(data = data.frame(x, y), aes(x = x , y = y)) +
  geom_point() +
  geom_smooth(method = "glm", method.args = list(family = "binomial"), se = TRUE) +
  labs(title = "Winter 20203", x = "precip", y = "LowLevel")
event_winter_20203_precip0 <- ggplot(data.frame(y, x) %>% filter(y == 0), aes(x = x, y = y)) +
  geom_boxplot() +
  labs(title = "LowLevel = 0", x = "tmp", y = "Low Level")
event_winter_20203_precip1 <- ggplot(data.frame(y, x) %>% filter(y == 1), aes(x = x, y = y)) +
  geom_boxplot() +
  labs(title = "LowLevel = 1", x = "tmp", y = "Low Level")
grid_event_winter_20203_precip <- grid.arrange(event_winter_20203_precip, event_winter_20203_precip0, event_winter_20203_precip1, ncol = 1)
grid.arrange(grid_event_winter_10304_precip, grid_event_winter_11502_precip, grid_event_winter_20203_precip, ncol = 2)
grid.arrange(grid_event_summer_10304_precip, grid_event_summer_11502_precip, grid_event_summer_20203_precip, grid_event_winter_10304_precip, grid_event_winter_11502_precip, grid_event_winter_20203_precip, ncol = 3)

```

## radiation
``` {r}
# summer

y <- hydro_summer_10304$lowlevel
x <- hydro_summer_10304$glorad
cor(x, y, method = "pearson")
glm(y~x, data = data.frame(x, y), family = binomial())
event_summer_10304_glorad<- ggplot(data = data.frame(x, y), aes(x = x , y = y)) +
  geom_point() +
  geom_smooth(method = "glm", method.args = list(family = "binomial"), se = TRUE) +
  labs(title = "Summer 10304", x = "glorad", y = "LowLevel")
event_summer_10304_glorad0 <- ggplot(data.frame(y, x) %>% filter(y == 0), aes(x = x, y = y)) +
  geom_boxplot() +
  labs(title = "LowLevel = 0", x = "glorad", y = "Low Level")
event_summer_10304_glorad1 <- ggplot(data.frame(y, x) %>% filter(y == 1), aes(x = x, y = y)) +
  geom_boxplot() +
  labs(title = "LowLevel = 1", x = "glorad", y = "Low Level")
grid_event_summer_10304_glorad <- grid.arrange(event_summer_10304_glorad, event_summer_10304_glorad0, event_summer_10304_glorad1, ncol = 1)


y <- hydro_summer_11502$lowlevel
x <- hydro_summer_11502$glorad
cor(x, y, method = "pearson")
glm(y~x, data = data.frame(x, y), family = binomial())
event_summer_11502_glorad <- ggplot(data = data.frame(x, y), aes(x = x , y = y)) +
  geom_point() +
  geom_smooth(method = "glm", method.args = list(family = "binomial"), se = TRUE) +
  labs(title = "Summer 11502", x = "glorad", y = "LowLevel")
event_summer_11502_glorad0 <- ggplot(data.frame(y, x) %>% filter(y == 0), aes(x = x, y = y)) +
  geom_boxplot() +
  labs(title = "LowLevel = 0", x = "glorad", y = "Low Level")
event_summer_11502_glorad1 <- ggplot(data.frame(y, x) %>% filter(y == 1), aes(x = x, y = y)) +
  geom_boxplot() +
  labs(title = "LowLevel = 1", x = "glorad", y = "Low Level")
grid_event_summer_11502_glorad <- grid.arrange(event_summer_11502_glorad, event_summer_11502_glorad0, event_summer_11502_glorad1, ncol = 1)


y <- hydro_summer_20203$lowlevel
x <- hydro_summer_20203$glorad
cor(x, y, method = "pearson")
glm(y~x, data = data.frame(x, y), family = binomial())
event_summer_20203_glorad <- ggplot(data = data.frame(x, y), aes(x = x , y = y)) +
  geom_point() +
  geom_smooth(method = "glm", method.args = list(family = "binomial"), se = TRUE) +
  labs(title = "Summer 20203", x = "glorad", y = "LowLevel")
event_summer_20203_glorad0 <- ggplot(data.frame(y, x) %>% filter(y == 0), aes(x = x, y = y)) +
  geom_boxplot() +
  labs(title = "LowLevel = 0", x = "glorad", y = "Low Level")
event_summer_20203_glorad1 <- ggplot(data.frame(y, x) %>% filter(y == 1), aes(x = x, y = y)) +
  geom_boxplot() +
  labs(title = "LowLevel = 1", x = "glorad", y = "Low Level")
grid_event_summer_20203_glorad <- grid.arrange(event_summer_20203_glorad, event_summer_20203_glorad0, event_summer_20203_glorad1, ncol = 1)

grid.arrange(grid_event_summer_10304_glorad, grid_event_summer_11502_glorad, grid_event_summer_20203_glorad, ncol = 2)

# winter
y <- hydro_winter_10304$lowlevel
x <- hydro_winter_10304$glorad
cor(x, y, method = "pearson")
glm(y~x, data = data.frame(x, y), family = binomial())
event_winter_10304_glorad <- ggplot(data = data.frame(x, y), aes(x = x , y = y)) +
  geom_point() +
  geom_smooth(method = "glm", method.args = list(family = "binomial"), se = TRUE) +
  labs(title = "Winter 10304", x = "glorad", y = "LowLevel")
event_winter_10304_glorad0 <- ggplot(data.frame(y, x) %>% filter(y == 0), aes(x = x, y = y)) +
  geom_boxplot() +
  labs(title = "LowLevel = 0", x = "glorad", y = "Low Level")
event_winter_10304_glorad1 <- ggplot(data.frame(y, x) %>% filter(y == 1), aes(x = x, y = y)) +
  geom_boxplot() +
  labs(title = "LowLevel = 1", x = "glorad", y = "Low Level")
grid_event_winter_10304_glorad <- grid.arrange(event_winter_10304_glorad, event_winter_10304_glorad0, event_winter_10304_glorad1, ncol = 1)                           


y <- hydro_winter_11502$lowlevel
x <- hydro_winter_11502$glorad
cor(x, y, method = "pearson")
glm(y~x, data = data.frame(x, y), family = binomial())
event_winter_11502_glorad <- ggplot(data = data.frame(x, y), aes(x = x , y = y)) +
  geom_point() +
  geom_smooth(method = "glm", method.args = list(family = "binomial"), se = TRUE) +
  labs(title = "Winter 11502", x = "glorad", y = "LowLevel")
event_winter_11502_glorad0 <- ggplot(data.frame(y, x) %>% filter(y == 0), aes(x = x, y = y)) +
  geom_boxplot() +
  labs(title = "LowLevel = 0", x = "glorad", y = "Low Level")
event_winter_11502_glorad1 <- ggplot(data.frame(y, x) %>% filter(y == 1), aes(x = x, y = y)) +
  geom_boxplot() +
  labs(title = "LowLevel = 1", x = "glorad", y = "Low Level")
grid_event_winter_11502_glorad <- grid.arrange(event_winter_11502_glorad, event_winter_11502_glorad0, event_winter_11502_glorad1, ncol = 1)


y <- hydro_winter_20203$lowlevel
x <- hydro_winter_20203$glorad
cor(x, y, method = "pearson")
glm(y~x, data = data.frame(x, y), family = binomial())
event_winter_20203_glorad <- ggplot(data = data.frame(x, y), aes(x = x , y = y)) +
  geom_point() +
  geom_smooth(method = "glm", method.args = list(family = "binomial"), se = TRUE) +
  labs(title = "Winter 20203", x = "glorad", y = "LowLevel")
event_winter_20203_glorad0 <- ggplot(data.frame(y, x) %>% filter(y == 0), aes(x = x, y = y)) +
  geom_boxplot() +
  labs(title = "LowLevel = 0", x = "glorad", y = "Low Level")
event_winter_20203_glorad1 <- ggplot(data.frame(y, x) %>% filter(y == 1), aes(x = x, y = y)) +
  geom_boxplot() +
  labs(title = "LowLevel = 1", x = "glorad", y = "Low Level")
grid_event_winter_20203_glorad <- grid.arrange(event_winter_20203_glorad, event_winter_20203_glorad0, event_winter_20203_glorad1, ncol = 1)
grid.arrange(grid_event_winter_10304_glorad, grid_event_winter_11502_glorad, grid_event_winter_20203_glorad, ncol = 2)
grid.arrange(grid_event_summer_10304_glorad, grid_event_summer_11502_glorad, grid_event_summer_20203_glorad, grid_event_winter_10304_glorad, grid_event_winter_11502_glorad, grid_event_winter_20203_glorad, ncol = 3)

```

## rel.humidity
``` {r}
# summer

y <- hydro_summer_10304$lowlevel
x <- hydro_summer_10304$relhum
cor(x, y, method = "pearson")
glm(y~x, data = data.frame(x, y), family = binomial())
event_summer_10304_relhum <- ggplot(data = data.frame(x, y), aes(x = x , y = y)) +
  geom_point() +
  geom_smooth(method = "glm", method.args = list(family = "binomial"), se = TRUE) +
  labs(title = "Summer 10304", x = "relhum", y = "LowLevel")
event_summer_10304_relhum0 <- ggplot(data.frame(y, x) %>% filter(y == 0), aes(x = x, y = y)) +
  geom_boxplot() +
  labs(title = "LowLevel = 0", x = "relhum", y = "Low Level")
event_summer_10304_relhum1 <- ggplot(data.frame(y, x) %>% filter(y == 1), aes(x = x, y = y)) +
  geom_boxplot() +
  labs(title = "LowLevel = 1", x = "relhum", y = "Low Level")
grid_event_summer_10304_relhum <- grid.arrange(event_summer_10304_relhum, event_summer_10304_relhum0, event_summer_10304_relhum1, ncol = 1)

y <- hydro_summer_11502$lowlevel
x <- hydro_summer_11502$relhum
cor(x, y, method = "pearson")
glm(y~x, data = data.frame(x, y), family = binomial())
event_summer_11502_relhum <- ggplot(data = data.frame(x, y), aes(x = x , y = y)) +
  geom_point() +
  geom_smooth(method = "glm", method.args = list(family = "binomial"), se = TRUE) +
  labs(title = "Summer 11502", x = "relhum", y = "LowLevel")
event_summer_11502_relhum0 <- ggplot(data.frame(y, x) %>% filter(y == 0), aes(x = x, y = y)) +
  geom_boxplot() +
  labs(title = "LowLevel = 0", x = "relhum", y = "Low Level")
event_summer_11502_relhum1 <- ggplot(data.frame(y, x) %>% filter(y == 1), aes(x = x, y = y)) +
  geom_boxplot() +
  labs(title = "LowLevel = 1", x = "relhum", y = "Low Level")
grid_event_summer_11502_relhum <- grid.arrange(event_summer_11502_relhum, event_summer_11502_relhum0, event_summer_11502_relhum1, ncol = 1)


y <- hydro_summer_20203$lowlevel
x <- hydro_summer_20203$relhum
cor(x, y, method = "pearson")
glm(y~x, data = data.frame(x, y), family = binomial())
event_summer_20203_relhum <- ggplot(data = data.frame(x, y), aes(x = x , y = y)) +
  geom_point() +
  geom_smooth(method = "glm", method.args = list(family = "binomial"), se = TRUE) +
  labs(title = "Summer 20203", x = "relhum", y = "LowLevel")
event_summer_20203_relhum0 <- ggplot(data.frame(y, x) %>% filter(y == 0), aes(x = x, y = y)) +
  geom_boxplot() +
  labs(title = "LowLevel = 0", x = "relhum", y = "Low Level")
event_summer_20203_relhum1 <- ggplot(data.frame(y, x) %>% filter(y == 1), aes(x = x, y = y)) +
  geom_boxplot() +
  labs(title = "LowLevel = 1", x = "relhum", y = "Low Level")
grid_event_summer_20203_relhum <- grid.arrange(event_summer_20203_relhum, event_summer_20203_relhum0, event_summer_20203_relhum1, ncol = 1)

grid.arrange(grid_event_summer_10304_relhum, grid_event_summer_11502_relhum, grid_event_summer_20203_relhum, ncol = 2)

# winter
y <- hydro_winter_10304$lowlevel
x <- hydro_winter_10304$relhum
cor(x, y, method = "pearson")
glm(y~x, data = data.frame(x, y), family = binomial())
event_winter_10304_relhum <- ggplot(data = data.frame(x, y), aes(x = x , y = y)) +
  geom_point() +
  geom_smooth(method = "glm", method.args = list(family = "binomial"), se = TRUE) +
  labs(title = "Winter 10304", x = "relhum", y = "LowLevel")
event_winter_10304_relhum0 <- ggplot(data.frame(y, x) %>% filter(y == 0), aes(x = x, y = y)) +
  geom_boxplot() +
  labs(title = "LowLevel = 0", x = "relhum", y = "Low Level")
event_winter_10304_relhum1 <- ggplot(data.frame(y, x) %>% filter(y == 1), aes(x = x, y = y)) +
  geom_boxplot() +
  labs(title = "LowLevel = 1", x = "relhum", y = "Low Level")
grid_event_winter_10304_relhum <- grid.arrange(event_winter_10304_relhum, event_winter_10304_relhum0, event_winter_10304_relhum1, ncol = 1)       

y <- hydro_winter_11502$lowlevel
x <- hydro_winter_11502$relhum
cor(x, y, method = "pearson")
glm(y~x, data = data.frame(x, y), family = binomial())
event_winter_11502_relhum <- ggplot(data = data.frame(x, y), aes(x = x , y = y)) +
  geom_point() +
  geom_smooth(method = "glm", method.args = list(family = "binomial"), se = TRUE) +
  labs(title = "Winter 11502", x = "relhum", y = "LowLevel")
event_winter_11502_relhum0 <- ggplot(data.frame(y, x) %>% filter(y == 0), aes(x = x, y = y)) +
  geom_boxplot() +
  labs(title = "LowLevel = 0", x = "relhum", y = "Low Level")
event_winter_11502_relhum1 <- ggplot(data.frame(y, x) %>% filter(y == 1), aes(x = x, y = y)) +
  geom_boxplot() +
  labs(title = "LowLevel = 1", x = "relhum", y = "Low Level")
grid_event_winter_11502_relhum <- grid.arrange(event_winter_11502_relhum, event_winter_11502_relhum0, event_winter_11502_relhum1, ncol = 1)


y <- hydro_winter_20203$lowlevel
x <- hydro_winter_20203$relhum
cor(x, y, method = "pearson")
glm(y~x, data = data.frame(x, y), family = binomial())
event_winter_20203_relhum <- ggplot(data = data.frame(x, y), aes(x = x , y = y)) +
  geom_point() +
  geom_smooth(method = "glm", method.args = list(family = "binomial"), se = TRUE) +
  labs(title = "Winter 20203", x = "relhum", y = "LowLevel")
event_winter_20203_relhum0 <- ggplot(data.frame(y, x) %>% filter(y == 0), aes(x = x, y = y)) +
  geom_boxplot() +
  labs(title = "LowLevel = 0", x = "relhum", y = "Low Level")
event_winter_20203_relhum1 <- ggplot(data.frame(y, x) %>% filter(y == 1), aes(x = x, y = y)) +
  geom_boxplot() +
  labs(title = "LowLevel = 1", x = "relhum", y = "Low Level")
grid_event_winter_20203_relhum <- grid.arrange(event_winter_20203_relhum, event_winter_20203_relhum0, event_winter_20203_relhum1, ncol = 1)

grid.arrange(grid_event_winter_10304_relhum, grid_event_winter_11502_relhum, grid_event_winter_20203_relhum, ncol = 2)
grid.arrange(grid_event_summer_10304_relhum, grid_event_summer_11502_relhum, grid_event_summer_20203_relhum, grid_event_winter_10304_relhum, grid_event_winter_11502_relhum, grid_event_winter_20203_relhum, ncol = 3)


```

## groundwater depth
``` {r}
# summer

y <- hydro_summer_10304$lowlevel
x <- hydro_summer_10304$groundwaterdepth
cor(x, y, method = "pearson")
glm(y~x, data = data.frame(x, y), family = binomial())
event_summer_10304_groundwaterdepth <- ggplot(data = data.frame(x, y), aes(x = x , y = y)) +
  geom_point() +
  geom_smooth(method = "glm", method.args = list(family = "binomial"), se = TRUE) +
  labs(title = "Summer 10304", x = "groundwaterdepth", y = "LowLevel")
event_summer_10304_groundwaterdepth0 <- ggplot(data.frame(y, x) %>% filter(y == 0), aes(x = x, y = y)) +
  geom_boxplot() +
  labs(title = "LowLevel = 0", x = "groundwaterdepth", y = "Low Level")
event_summer_10304_groundwaterdepth1 <- ggplot(data.frame(y, x) %>% filter(y == 1), aes(x = x, y = y)) +
  geom_boxplot() +
  labs(title = "LowLevel = 1", x = "groundwaterdepth", y = "Low Level")
grid_event_summer_10304_groundwaterdepth <- grid.arrange(event_summer_10304_groundwaterdepth, event_summer_10304_groundwaterdepth0, event_summer_10304_groundwaterdepth1, ncol = 1)



y <- hydro_summer_11502$lowlevel
x <- hydro_summer_11502$groundwaterdepth
cor(x, y, method = "pearson")
glm(y~x, data = data.frame(x, y), family = binomial())
event_summer_11502_groundwaterdepth <- ggplot(data = data.frame(x, y), aes(x = x , y = y)) +
  geom_point() +
  geom_smooth(method = "glm", method.args = list(family = "binomial"), se = TRUE) +
  labs(title = "Summer 11502", x = "groundwaterdepth", y = "LowLevel")
event_summer_11502_groundwaterdepth0 <- ggplot(data.frame(y, x) %>% filter(y == 0), aes(x = x, y = y)) +
  geom_boxplot() +
  labs(title = "LowLevel = 0", x = "groundwaterdepth", y = "Low Level")
event_summer_11502_groundwaterdepth1 <- ggplot(data.frame(y, x) %>% filter(y == 1), aes(x = x, y = y)) +
  geom_boxplot() +
  labs(title = "LowLevel = 1", x = "groundwaterdepth", y = "Low Level")
grid_event_summer_11502_groundwaterdepth<- grid.arrange(event_summer_11502_groundwaterdepth, event_summer_11502_groundwaterdepth0, event_summer_11502_groundwaterdepth1, ncol = 1)



y <- hydro_summer_20203$lowlevel
x <- hydro_summer_20203$groundwaterdepth
cor(x, y, method = "pearson")
glm(y~x, data = data.frame(x, y), family = binomial())
event_summer_20203_groundwaterdepth <- ggplot(data = data.frame(x, y), aes(x = x , y = y)) +
  geom_point() +
  geom_smooth(method = "glm", method.args = list(family = "binomial"), se = TRUE) +
  labs(title = "Summer 20203", x = "groundwaterdepth", y = "LowLevel")
event_summer_20203_groundwaterdepth0 <- ggplot(data.frame(y, x) %>% filter(y == 0), aes(x = x, y = y)) +
  geom_boxplot() +
  labs(title = "LowLevel = 0", x = "groundwaterdepth", y = "Low Level")
event_summer_20203_groundwaterdepth1 <- ggplot(data.frame(y, x) %>% filter(y == 1), aes(x = x, y = y)) +
  geom_boxplot() +
  labs(title = "LowLevel = 1", x = "groundwaterdepth", y = "Low Level")
grid_event_summer_20203_groundwaterdepth <- grid.arrange(event_summer_20203_groundwaterdepth, event_summer_20203_groundwaterdepth0, event_summer_20203_groundwaterdepth1, ncol = 1)

grid.arrange(grid_event_summer_10304_groundwaterdepth, grid_event_summer_11502_groundwaterdepth, grid_event_summer_20203_groundwaterdepth, ncol = 2)

# winter
y <- hydro_winter_10304$lowlevel
x <- hydro_winter_10304$groundwaterdepth
cor(x, y, method = "pearson")
glm(y~x, data = data.frame(x, y), family = binomial())
event_winter_10304_groundwaterdepth <- ggplot(data = data.frame(x, y), aes(x = x , y = y)) +
  geom_point() +
  geom_smooth(method = "glm", method.args = list(family = "binomial"), se = TRUE) +
  labs(title = "Winter 10304", x = "groundwaterdepth", y = "LowLevel")
event_winter_10304_groundwaterdepth0 <- ggplot(data.frame(y, x) %>% filter(y == 0), aes(x = x, y = y)) +
  geom_boxplot() +
  labs(title = "LowLevel = 0", x = "groundwaterdepth", y = "Low Level")
event_winter_10304_groundwaterdepth1 <- ggplot(data.frame(y, x) %>% filter(y == 1), aes(x = x, y = y)) +
  geom_boxplot() +
  labs(title = "LowLevel = 1", x = "groundwaterdepth", y = "Low Level")
grid_event_winter_10304_groundwaterdepth <- grid.arrange(event_winter_10304_groundwaterdepth, event_winter_10304_groundwaterdepth0, event_winter_10304_groundwaterdepth1, ncol = 1)  


y <- hydro_winter_11502$lowlevel
x <- hydro_winter_11502$groundwaterdepth
cor(x, y, method = "pearson")
glm(y~x, data = data.frame(x, y), family = binomial())
event_winter_11502_groundwaterdepth <- ggplot(data = data.frame(x, y), aes(x = x , y = y)) +
  geom_point() +
  geom_smooth(method = "glm", method.args = list(family = "binomial"), se = TRUE) +
  labs(title = "Winter 11502", x = "groundwaterdepth", y = "LowLevel")
event_winter_11502_groundwaterdepth0 <- ggplot(data.frame(y, x) %>% filter(y == 0), aes(x = x, y = y)) +
  geom_boxplot() +
  labs(title = "LowLevel = 0", x = "groundwaterdepth", y = "Low Level")
event_winter_11502_groundwaterdepth1 <- ggplot(data.frame(y, x) %>% filter(y == 1), aes(x = x, y = y)) +
  geom_boxplot() +
  labs(title = "LowLevel = 1", x = "groundwaterdepth", y = "Low Level")
grid_event_winter_11502_groundwaterdepth <- grid.arrange(event_winter_11502_groundwaterdepth, event_winter_11502_groundwaterdepth0, event_winter_11502_groundwaterdepth1, ncol = 1)



y <- hydro_winter_20203$lowlevel
x <- hydro_winter_20203$groundwaterdepth
cor(x, y, method = "pearson")
glm(y~x, data = data.frame(x, y), family = binomial())
event_winter_20203_groundwaterdepth <- ggplot(data = data.frame(x, y), aes(x = x , y = y)) +
  geom_point() +
  geom_smooth(method = "glm", method.args = list(family = "binomial"), se = TRUE) +
  labs(title = "Winter 20203", x = "groundwaterdepth", y = "LowLevel")
event_winter_20203_groundwaterdepth0 <- ggplot(data.frame(y, x) %>% filter(y == 0), aes(x = x, y = y)) +
  geom_boxplot() +
  labs(title = "LowLevel = 0", x = "groundwaterdepth", y = "Low Level")
event_winter_20203_groundwaterdepth1 <- ggplot(data.frame(y, x) %>% filter(y == 1), aes(x = x, y = y)) +
  geom_boxplot() +
  labs(title = "LowLevel = 1", x = "groundwaterdepth", y = "Low Level")
grid_event_winter_20203_groundwaterdepth <- grid.arrange(event_winter_20203_groundwaterdepth, event_winter_20203_groundwaterdepth0, event_winter_20203_groundwaterdepth1, ncol = 1)

grid.arrange(grid_event_winter_10304_groundwaterdepth, grid_event_winter_11502_groundwaterdepth, grid_event_winter_20203_groundwaterdepth, ncol = 2)
grid.arrange(grid_event_summer_10304_groundwaterdepth, grid_event_summer_11502_groundwaterdepth, grid_event_summer_20203_groundwaterdepth, grid_event_winter_10304_groundwaterdepth, grid_event_winter_11502_groundwaterdepth, grid_event_winter_20203_groundwaterdepth, ncol = 3)


```

## soilwater
``` {r}
# summer

y <- hydro_summer_10304$lowlevel
x <- hydro_summer_10304$soilwater
cor(x, y, method = "pearson")
glm(y~x, data = data.frame(x, y), family = binomial())
event_summer_10304_soilwater <- ggplot(data = data.frame(x, y), aes(x = x , y = y)) +
  geom_point() +
  geom_smooth(method = "glm", method.args = list(family = "binomial"), se = TRUE) +
  labs(title = "Summer 10304", x = "soilwater", y = "LowLevel")
event_summer_10304_soilwater0 <- ggplot(data.frame(y, x) %>% filter(y == 0), aes(x = x, y = y)) +
  geom_boxplot() +
  labs(title = "LowLevel = 0", x = "soilwater", y = "Low Level")
event_summer_10304_soilwater1 <- ggplot(data.frame(y, x) %>% filter(y == 1), aes(x = x, y = y)) +
  geom_boxplot() +
  labs(title = "LowLevel = 1", x = "soilwater", y = "Low Level")
grid_event_summer_10304_soilwater <- grid.arrange(event_summer_10304_soilwater, event_summer_10304_soilwater0, event_summer_10304_soilwater1, ncol = 1)


y <- hydro_summer_11502$lowlevel
x <- hydro_summer_11502$soilwater
cor(x, y, method = "pearson")
glm(y~x, data = data.frame(x, y), family = binomial())
event_summer_11502_soilwater <- ggplot(data = data.frame(x, y), aes(x = x , y = y)) +
  geom_point() +
  geom_smooth(method = "glm", method.args = list(family = "binomial"), se = TRUE) +
  labs(title = "Summer 11502", x = "soilwater", y = "LowLevel")
event_summer_11502_soilwater0 <- ggplot(data.frame(y, x) %>% filter(y == 0), aes(x = x, y = y)) +
  geom_boxplot() +
  labs(title = "LowLevel = 0", x = "soilwater", y = "Low Level")
event_summer_11502_soilwater1 <- ggplot(data.frame(y, x) %>% filter(y == 1), aes(x = x, y = y)) +
  geom_boxplot() +
  labs(title = "LowLevel = 1", x = "soilwater", y = "Low Level")
grid_event_summer_11502_soilwater <- grid.arrange(event_summer_11502_soilwater, event_summer_11502_soilwater0, event_summer_11502_soilwater1, ncol = 1)


y <- hydro_summer_20203$lowlevel
x <- hydro_summer_20203$soilwater
cor(x, y, method = "pearson")
glm(y~x, data = data.frame(x, y), family = binomial())
event_summer_20203_soilwater <- ggplot(data = data.frame(x, y), aes(x = x , y = y)) +
  geom_point() +
  geom_smooth(method = "glm", method.args = list(family = "binomial"), se = TRUE) +
  labs(title = "Summer 20203", x = "soilwater", y = "LowLevel")
event_summer_20203_soilwater0 <- ggplot(data.frame(y, x) %>% filter(y == 0), aes(x = x, y = y)) +
  geom_boxplot() +
  labs(title = "LowLevel = 0", x = "soilwater", y = "Low Level")
event_summer_20203_soilwater1 <- ggplot(data.frame(y, x) %>% filter(y == 1), aes(x = x, y = y)) +
  geom_boxplot() +
  labs(title = "LowLevel = 1", x = "soilwater", y = "Low Level")
grid_event_summer_20203_soilwater <- grid.arrange(event_summer_20203_soilwater, event_summer_20203_soilwater0, event_summer_20203_soilwater1, ncol = 1)

grid.arrange(grid_event_summer_10304_soilwater, grid_event_summer_11502_soilwater, grid_event_summer_20203_soilwater, ncol = 2)

# winter
y <- hydro_winter_10304$lowlevel
x <- hydro_winter_10304$soilwater
cor(x, y, method = "pearson")
glm(y~x, data = data.frame(x, y), family = binomial())
event_winter_10304_soilwater <- ggplot(data = data.frame(x, y), aes(x = x , y = y)) +
  geom_point() +
  geom_smooth(method = "glm", method.args = list(family = "binomial"), se = TRUE) +
  labs(title = "Winter 10304", x = "soilwater", y = "LowLevel")
event_winter_10304_soilwater0 <- ggplot(data.frame(y, x) %>% filter(y == 0), aes(x = x, y = y)) +
  geom_boxplot() +
  labs(title = "LowLevel = 0", x = "soilwater", y = "Low Level")
event_winter_10304_soilwater1 <- ggplot(data.frame(y, x) %>% filter(y == 1), aes(x = x, y = y)) +
  geom_boxplot() +
  labs(title = "LowLevel = 1", x = "soilwater", y = "Low Level")
grid_event_winter_10304_soilwater <- grid.arrange(event_winter_10304_soilwater, event_winter_10304_soilwater0, event_winter_10304_soilwater1, ncol = 1)       

y <- hydro_winter_11502$lowlevel
x <- hydro_winter_11502$soilwater
cor(x, y, method = "pearson")
glm(y~x, data = data.frame(x, y), family = binomial())
event_winter_11502_soilwater <- ggplot(data = data.frame(x, y), aes(x = x , y = y)) +
  geom_point() +
  geom_smooth(method = "glm", method.args = list(family = "binomial"), se = TRUE) +
  labs(title = "Winter 11502", x = "soilwater", y = "LowLevel")
event_winter_11502_soilwater0 <- ggplot(data.frame(y, x) %>% filter(y == 0), aes(x = x, y = y)) +
  geom_boxplot() +
  labs(title = "LowLevel = 0", x = "soilwater", y = "Low Level")
event_winter_11502_soilwater1 <- ggplot(data.frame(y, x) %>% filter(y == 1), aes(x = x, y = y)) +
  geom_boxplot() +
  labs(title = "LowLevel = 1", x = "soilwater", y = "Low Level")
grid_event_winter_11502_soilwater <- grid.arrange(event_winter_11502_soilwater, event_winter_11502_soilwater0, event_winter_11502_soilwater1, ncol = 1)


y <- hydro_winter_20203$lowlevel
x <- hydro_winter_20203$soilwater
cor(x, y, method = "pearson")
glm(y~x, data = data.frame(x, y), family = binomial())
event_winter_20203_soilwater <- ggplot(data = data.frame(x, y), aes(x = x , y = y)) +
  geom_point() +
  geom_smooth(method = "glm", method.args = list(family = "binomial"), se = TRUE) +
  labs(title = "Winter 20203", x = "soilwater", y = "LowLevel")
event_winter_20203_soilwater0 <- ggplot(data.frame(y, x) %>% filter(y == 0), aes(x = x, y = y)) +
  geom_boxplot() +
  labs(title = "LowLevel = 0", x = "soilwater", y = "Low Level")
event_winter_20203_soilwater1 <- ggplot(data.frame(y, x) %>% filter(y == 1), aes(x = x, y = y)) +
  geom_boxplot() +
  labs(title = "LowLevel = 1", x = "soilwater", y = "Low Level")
grid_event_winter_20203_soilwater <- grid.arrange(event_winter_20203_soilwater, event_winter_20203_soilwater0, event_winter_20203_soilwater1, ncol = 1)

grid.arrange(grid_event_summer_10304_soilwater, grid_event_summer_11502_soilwater, grid_event_summer_20203_soilwater, ncol = 2)
grid.arrange(grid_event_summer_10304_soilwater, grid_event_summer_11502_soilwater, grid_event_summer_20203_soilwater, grid_event_summer_10304_soilwater, grid_event_summer_11502_soilwater, grid_event_summer_20203_soilwater, ncol = 3)

```

## Snow
``` {r}

# summer

y <- hydro_summer_10304$lowlevel
x <- hydro_summer_10304$snowstorage
cor(x, y, method = "pearson")
glm(y~x, data = data.frame(x, y), family = binomial())
event_summer_10304_snowstorage <- ggplot(data = data.frame(x, y), aes(x = x , y = y)) +
  geom_point() +
  geom_smooth(method = "glm", method.args = list(family = "binomial"), se = TRUE) +
  labs(title = "Summer 10304", x = "snowstorage", y = "LowLevel")
event_summer_10304_snowstorage0 <- ggplot(data.frame(y, x) %>% filter(y == 0), aes(x = x, y = y)) +
  geom_boxplot() +
  labs(title = "LowLevel = 0", x = "snowstorage", y = "Low Level")
event_summer_10304_snowstorage1 <- ggplot(data.frame(y, x) %>% filter(y == 1), aes(x = x, y = y)) +
  geom_boxplot() +
  labs(title = "LowLevel = 1", x = "snowstorage", y = "Low Level")
grid_event_summer_10304_snowstorage <- grid.arrange(event_summer_10304_snowstorage, event_summer_10304_snowstorage0, event_summer_10304_snowstorage1, ncol = 1)


y <- hydro_summer_11502$lowlevel
x <- hydro_summer_11502$snowstorage
cor(x, y, method = "pearson")
glm(y~x, data = data.frame(x, y), family = binomial())
event_summer_11502_snowstorage <- ggplot(data = data.frame(x, y), aes(x = x , y = y)) +
  geom_point() +
  geom_smooth(method = "glm", method.args = list(family = "binomial"), se = TRUE) +
  labs(title = "Summer 11502", x = "snowstorage", y = "LowLevel")
event_summer_11502_snowstorage0 <- ggplot(data.frame(y, x) %>% filter(y == 0), aes(x = x, y = y)) +
  geom_boxplot() +
  labs(title = "LowLevel = 0", x = "snowstorage", y = "Low Level")
event_summer_11502_snowstorage1 <- ggplot(data.frame(y, x) %>% filter(y == 1), aes(x = x, y = y)) +
  geom_boxplot() +
  labs(title = "LowLevel = 1", x = "snowstorage", y = "Low Level")
grid_event_summer_11502_snowstorage <- grid.arrange(event_summer_11502_snowstorage, event_summer_11502_snowstorage0, event_summer_11502_snowstorage1, ncol = 1)


y <- hydro_summer_20203$lowlevel
x <- hydro_summer_20203$snowstorage
cor(x, y, method = "pearson")
glm(y~x, data = data.frame(x, y), family = binomial())
event_summer_20203_snowstorage <- ggplot(data = data.frame(x, y), aes(x = x , y = y)) +
  geom_point() +
  geom_smooth(method = "glm", method.args = list(family = "binomial"), se = TRUE) +
  labs(title = "Summer 20203", x = "snowstorage", y = "LowLevel")
event_summer_20203_snowstorage0 <- ggplot(data.frame(y, x) %>% filter(y == 0), aes(x = x, y = y)) +
  geom_boxplot() +
  labs(title = "LowLevel = 0", x = "snowstorage", y = "Low Level")
event_summer_20203_snowstorage1 <- ggplot(data.frame(y, x) %>% filter(y == 1), aes(x = x, y = y)) +
  geom_boxplot() +
  labs(title = "LowLevel = 1", x = "snowstorage", y = "Low Level")
grid_event_summer_20203_snowstorage <- grid.arrange(event_summer_20203_snowstorage, event_summer_20203_snowstorage0, event_summer_20203_snowstorage1, ncol = 1)

grid.arrange(grid_event_summer_10304_snowstorage, grid_event_summer_11502_snowstorage, grid_event_summer_20203_snowstorage, ncol = 2)

# winter
y <- hydro_winter_10304$lowlevel
x <- hydro_winter_10304$snowstorage
cor(x, y, method = "pearson")
glm(y~x, data = data.frame(x, y), family = binomial())
event_winter_10304_snowstorage <- ggplot(data = data.frame(x, y), aes(x = x , y = y)) +
  geom_point() +
  geom_smooth(method = "glm", method.args = list(family = "binomial"), se = TRUE) +
  labs(title = "Winter 10304", x = "snowstorage", y = "LowLevel")
event_winter_10304_snowstorage0 <- ggplot(data.frame(y, x) %>% filter(y == 0), aes(x = x, y = y)) +
  geom_boxplot() +
  labs(title = "LowLevel = 0", x = "snowstorage", y = "Low Level")
event_winter_10304_snowstorage1 <- ggplot(data.frame(y, x) %>% filter(y == 1), aes(x = x, y = y)) +
  geom_boxplot() +
  labs(title = "LowLevel = 1", x = "snowstorage", y = "Low Level")
grid_event_winter_10304_snowstorage <- grid.arrange(event_winter_10304_snowstorage, event_winter_10304_snowstorage0, event_winter_10304_snowstorage1, ncol = 1)       

y <- hydro_winter_11502$lowlevel
x <- hydro_winter_11502$snowstorage
cor(x, y, method = "pearson")
glm(y~x, data = data.frame(x, y), family = binomial())
event_winter_11502_snowstorage <- ggplot(data = data.frame(x, y), aes(x = x , y = y)) +
  geom_point() +
  geom_smooth(method = "glm", method.args = list(family = "binomial"), se = TRUE) +
  labs(title = "Winter 11502", x = "snowstorage", y = "LowLevel")
event_winter_11502_snowstorage0 <- ggplot(data.frame(y, x) %>% filter(y == 0), aes(x = x, y = y)) +
  geom_boxplot() +
  labs(title = "LowLevel = 0", x = "snowstorage", y = "Low Level")
event_winter_11502_snowstorage1 <- ggplot(data.frame(y, x) %>% filter(y == 1), aes(x = x, y = y)) +
  geom_boxplot() +
  labs(title = "LowLevel = 1", x = "snowstorage", y = "Low Level")
grid_event_winter_11502_snowstorage <- grid.arrange(event_winter_11502_snowstorage, event_winter_11502_snowstorage0, event_winter_11502_snowstorage1, ncol = 1)


y <- hydro_winter_20203$lowlevel
x <- hydro_winter_20203$snowstorage
cor(x, y, method = "pearson")
glm(y~x, data = data.frame(x, y), family = binomial())
event_winter_20203_snowstorage <- ggplot(data = data.frame(x, y), aes(x = x , y = y)) +
  geom_point() +
  geom_smooth(method = "glm", method.args = list(family = "binomial"), se = TRUE) +
  labs(title = "Winter 20203", x = "snowstorage", y = "LowLevel")
event_winter_20203_snowstorage0 <- ggplot(data.frame(y, x) %>% filter(y == 0), aes(x = x, y = y)) +
  geom_boxplot() +
  labs(title = "LowLevel = 0", x = "snowstorage", y = "Low Level")
event_winter_20203_snowstorage1 <- ggplot(data.frame(y, x) %>% filter(y == 1), aes(x = x, y = y)) +
  geom_boxplot() +
  labs(title = "LowLevel = 1", x = "snowstorage", y = "Low Level")
grid_event_winter_20203_snowstorage <- grid.arrange(event_winter_20203_snowstorage, event_winter_20203_snowstorage0, event_winter_20203_snowstorage1, ncol = 1)

grid.arrange(grid_event_summer_10304_snowstorage, grid_event_summer_11502_snowstorage, grid_event_summer_20203_snowstorage, ncol = 2)
grid.arrange(grid_event_summer_10304_snowstorage, grid_event_summer_11502_snowstorage, grid_event_summer_20203_snowstorage, grid_event_summer_10304_snowstorage, grid_event_summer_11502_snowstorage, grid_event_summer_20203_snowstorage, ncol = 3)

```

## Forst plot
``` {r}
## Sommer

varnames <- colnames(hydro_summer_10304)
varnames[15] <- "Lufttemperatur"
varnames[18] <- "Niederschlag"
varnames[16] <- "KurzwelligeStrahlung"
varnames[19] <- "RelativeLuftfeuchte"
varnames[17] <- "Grundwasserstand"
varnames[37] <- "relOberflächennaheBodenfeuchte"
varnames[20] <- "Schneespeicher"
varnames[26] <- "Niedriwasserevent"
varnames[22] <- "Wasserleitfähigkeit"
## 10304
hydro_summer_10304_named <- hydro_summer_10304
colnames(hydro_summer_10304_named) <- varnames
sommer10304GLM <- glm(Niedriwasserevent ~ Lufttemperatur + Niederschlag + KurzwelligeStrahlung + RelativeLuftfeuchte + Grundwasserstand + relOberflächennaheBodenfeuchte + Schneespeicher + Wasserleitfähigkeit, family = binomial(link='logit'), data = hydro_summer_10304_named)
#11502
hydro_summer_11502_named <- hydro_summer_11502
colnames(hydro_summer_11502_named) <- varnames
sommer11502GLM <- glm(Niedriwasserevent ~ Lufttemperatur + Niederschlag + KurzwelligeStrahlung + RelativeLuftfeuchte + Grundwasserstand + relOberflächennaheBodenfeuchte + Schneespeicher + Wasserleitfähigkeit, family = binomial(link='logit'), data = hydro_summer_11502_named)
#20203
hydro_summer_20203_named <- hydro_summer_20203
colnames(hydro_summer_20203_named) <- varnames
sommer20203GLM <- glm(Niedriwasserevent ~ Lufttemperatur + Niederschlag + KurzwelligeStrahlung + RelativeLuftfeuchte + Grundwasserstand + relOberflächennaheBodenfeuchte + Schneespeicher + Wasserleitfähigkeit, family = binomial(link='logit'), data = hydro_summer_20203_named)

##Winter
##10304
hydro_summer_10304_named <- hydro_summer_10304
colnames(hydro_summer_10304_named) <- varnames
winter10304GLM <- glm(Niedriwasserevent ~ Lufttemperatur + Niederschlag + KurzwelligeStrahlung + RelativeLuftfeuchte + Grundwasserstand + relOberflächennaheBodenfeuchte + Schneespeicher + Wasserleitfähigkeit, family = binomial(link='logit'), data = hydro_summer_10304_named)
#11502
hydro_winter_11502_named <- hydro_winter_11502
colnames(hydro_winter_11502_named) <- varnames
winter11502GLM <- glm(Niedriwasserevent ~ Lufttemperatur + Niederschlag + KurzwelligeStrahlung + RelativeLuftfeuchte + Grundwasserstand + relOberflächennaheBodenfeuchte + Schneespeicher + Wasserleitfähigkeit, family = binomial(link='logit'), data = hydro_winter_11502_named)
#20203
hydro_summer_20203_named <- hydro_summer_20203
colnames(hydro_summer_20203_named) <- varnames
winter20203GLM <- glm(Niedriwasserevent ~ Lufttemperatur + Niederschlag + KurzwelligeStrahlung + RelativeLuftfeuchte + Grundwasserstand + relOberflächennaheBodenfeuchte + Schneespeicher + Wasserleitfähigkeit, family = binomial(link='logit'), data = hydro_summer_20203_named)

all.models <- list()
all.models[[1]] <- sommer10304GLM
all.models[[2]] <- sommer11502GLM
all.models[[3]] <- sommer20203GLM
all.models[[4]] <- winter10304GLM
all.models[[5]] <- winter11502GLM
all.models[[6]] <- winter20203GLM
plot_models(all.models)
p <- last_plot() +
  scale_color_manual(values = c("yellow", "red", "blue","pink", "green", "purple" ), 
                     name = "Gewässer/Halbjahr", 
                     labels = c("Winter Fraenkische Saale Salz",  "Winter Iller Kempten",  "Winter Isar Mittenwald","Sommer Fraenkische Saale Salz", "Sommer Iller Kempten", "Sommer Isar Mittenwald")) +
  theme(legend.title = element_text(size = 12, face = "bold"))+
  ylab("Odds Ratio")

print(p)
```
## Rescaled forest plot
``` {r}
## Sommer

varnames <- colnames(hydro_summer_10304)
varnames1 <- colnames(hydro_summer_11502)
#Variablennamen für 10304
varnames[47] <- "avgLufttemperatur"
varnames[48] <- "avgNiederschlag"
varnames[49] <- "avgKurzwelligeStrahlung"
varnames[50] <- "avgRelativeLuftfeuchte"
varnames[52] <- "Grundwasserstand"
varnames[45] <- "avgOberflächennaheBodenfeuchte"
varnames[46] <- "avgSchneespeicher"
varnames[26] <- "Niedrigwasserevent"
varnames[51] <- "avgWasserleitfähigkeit"
varnames[53] <- "avgSchneeAbtauen"
# Variablennamen für 11502 und 20203
varnames1[47] <- "avgLufttemperatur"
varnames1[48] <- "avgNiederschlag"
varnames1[49] <- "avgKurzwelligeStrahlung"
varnames1[50] <- "avgRelativeLuftfeuchte"
varnames1[52] <- "Grundwasserstand"
varnames1[45] <- "avgOberflächennaheBodenfeuchte"
varnames1[46] <- "avgSchneespeicher"
varnames1[26] <- "Niedrigwasserevent"
varnames1[51] <- "avgWasserleitfähigkeit"
varnames1[53] <- "avgSchneeAbtauen"
## 10304
hydro_summer_10304_named <- hydro_summer_10304
colnames(hydro_summer_10304_named) <- varnames

sommer10304GLM <- glm(Niedrigwasserevent ~ avgLufttemperatur + avgNiederschlag + avgKurzwelligeStrahlung + avgRelativeLuftfeuchte + Grundwasserstand + avgOberflächennaheBodenfeuchte + avgSchneespeicher + avgWasserleitfähigkeit + avgSchneeAbtauen, family = binomial(link='logit'), data = hydro_summer_10304)
#11502
hydro_summer_11502_named <- hydro_summer_11502
colnames(hydro_summer_11502_named) <- varnames1

sommer11502GLM <- glm(Niedrigwasserevent ~ avgLufttemperatur + avgNiederschlag + avgKurzwelligeStrahlung + avgRelativeLuftfeuchte + Grundwasserstand + avgOberflächennaheBodenfeuchte + avgSchneespeicher + avgWasserleitfähigkeit + avgSchneeAbtauen, family = binomial(link='logit'), data = hydro_summer_11502_named)
#20203
hydro_summer_20203_named <- hydro_summer_20203
colnames(hydro_summer_20203_named) <- varnames1

sommer20203GLM <- glm(Niedrigwasserevent ~ avgLufttemperatur + avgNiederschlag + avgKurzwelligeStrahlung + avgRelativeLuftfeuchte + Grundwasserstand + avgOberflächennaheBodenfeuchte + avgSchneespeicher + avgWasserleitfähigkeit + avgSchneeAbtauen, family = binomial(link='logit'), data = hydro_summer_20203_named)

##Winter
##10304
hydro_winter_10304_named <- hydro_winter_10304
colnames(hydro_winter_10304_named) <- varnames

winter10304GLM <- glm(Niedrigwasserevent ~ avgLufttemperatur + avgNiederschlag + avgKurzwelligeStrahlung + avgRelativeLuftfeuchte + Grundwasserstand + avgOberflächennaheBodenfeuchte + avgSchneespeicher + avgWasserleitfähigkeit + avgSchneeAbtauen, family = binomial(link='logit'), data = hydro_winter_10304_named)

#11502
hydro_winter_11502_named <- hydro_winter_11502
colnames(hydro_winter_11502_named) <- varnames1

winter11502GLM <- glm(Niedrigwasserevent ~ avgLufttemperatur + avgNiederschlag + avgKurzwelligeStrahlung + avgRelativeLuftfeuchte + Grundwasserstand + avgOberflächennaheBodenfeuchte + avgSchneespeicher + avgWasserleitfähigkeit + avgSchneeAbtauen, family = binomial(link='logit'), data = hydro_winter_11502_named)
#20203
hydro_winter_20203_named <- hydro_winter_20203
colnames(hydro_winter_20203_named) <- varnames1

winter20203GLM <- glm(Niedrigwasserevent ~ avgLufttemperatur + avgNiederschlag + avgKurzwelligeStrahlung + avgRelativeLuftfeuchte + Grundwasserstand + avgOberflächennaheBodenfeuchte + avgSchneespeicher + avgWasserleitfähigkeit + avgSchneeAbtauen, family = binomial(link='logit'), data = hydro_winter_20203_named)

all.models <- list()
all.models[[1]] <- sommer10304GLM
all.models[[2]] <- sommer11502GLM
all.models[[3]] <- sommer20203GLM
all.models[[4]] <- winter10304GLM
all.models[[5]] <- winter11502GLM
all.models[[6]] <- winter20203GLM
plot_models(all.models)
p <- last_plot() +
  scale_color_manual(values = c("yellow", "red", "blue","pink", "green", "purple" ), 
                     name = "Gewässer/Halbjahr", 
                     labels = c("Winter Fraenkische Saale Salz",  "Winter Iller Kempten",  "Winter Isar Mittenwald","Sommer Fraenkische Saale Salz", "Sommer Iller Kempten", "Sommer Isar Mittenwald")) +
  theme(legend.title = element_text(size = 12, face = "bold"))+
  ylab("Odds Ratio")+
  labs(title = "Forest-Plot zum GLM der Reskalierten Variablen")

print(p)

### Nochmal ohne Schneespeicher, da starker Ausreißer
## 10304
sommer10304GLM_noSnow <- glm(Niedrigwasserevent ~ avgLufttemperatur + avgNiederschlag + avgKurzwelligeStrahlung + avgRelativeLuftfeuchte + Grundwasserstand + avgOberflächennaheBodenfeuchte + avgWasserleitfähigkeit + avgSchneeAbtauen, family = binomial(link='logit'), data = hydro_summer_10304)

#11502
sommer11502GLM_noSnow <- glm(Niedrigwasserevent ~ avgLufttemperatur + avgNiederschlag + avgKurzwelligeStrahlung + avgRelativeLuftfeuchte + Grundwasserstand + avgOberflächennaheBodenfeuchte + avgWasserleitfähigkeit + avgSchneeAbtauen, family = binomial(link='logit'), data = hydro_summer_11502_named)

#20203
sommer20203GLM_noSnow <- glm(Niedrigwasserevent ~ avgLufttemperatur + avgNiederschlag + avgKurzwelligeStrahlung + avgRelativeLuftfeuchte + Grundwasserstand + avgOberflächennaheBodenfeuchte + avgWasserleitfähigkeit + avgSchneeAbtauen, family = binomial(link='logit'), data = hydro_summer_20203_named)

##Winter
##10304

winter10304GLM_noSnow <- glm(Niedrigwasserevent ~ avgLufttemperatur + avgNiederschlag + avgKurzwelligeStrahlung + avgRelativeLuftfeuchte + Grundwasserstand + avgOberflächennaheBodenfeuchte + avgWasserleitfähigkeit + avgSchneeAbtauen, family = binomial(link='logit'), data = hydro_winter_10304_named)

#11502
winter11502GLM_noSnow <- glm(Niedrigwasserevent ~ avgLufttemperatur + avgNiederschlag + avgKurzwelligeStrahlung + avgRelativeLuftfeuchte + Grundwasserstand + avgOberflächennaheBodenfeuchte + avgWasserleitfähigkeit + avgSchneeAbtauen, family = binomial(link='logit'), data = hydro_winter_11502_named)

#20203
winter20203GLM_noSnow <- glm(Niedrigwasserevent ~ avgLufttemperatur + avgNiederschlag + avgKurzwelligeStrahlung + avgRelativeLuftfeuchte + Grundwasserstand + avgOberflächennaheBodenfeuchte + avgWasserleitfähigkeit + avgSchneeAbtauen, family = binomial(link='logit'), data = hydro_winter_20203_named)

all.models <- list()
all.models[[1]] <- sommer10304GLM_noSnow
all.models[[2]] <- sommer11502GLM_noSnow
all.models[[3]] <- sommer20203GLM_noSnow
all.models[[4]] <- winter10304GLM_noSnow
all.models[[5]] <- winter11502GLM_noSnow
all.models[[6]] <- winter20203GLM_noSnow
plot_models(all.models)
p1 <- last_plot() +
  scale_color_manual(values = c("yellow", "red", "blue","pink", "green", "purple" ), 
                     name = "Gewässer/Halbjahr", 
                     labels = c("Winter Fraenkische Saale Salz",  "Winter Iller Kempten",  "Winter Isar Mittenwald","Sommer Fraenkische Saale Salz", "Sommer Iller Kempten", "Sommer Isar Mittenwald")) +
  theme(legend.title = element_text(size = 12, face = "bold"))+
  ylab("Odds Ratio")+
  labs(title = "Forest-Plot zum GLM der Reskalierten Variablen",
       caption = "Ohne Schneespeicher")

print(p1)
```
### Glm - Linear trend over time 
```{r}
model_1 <- glm(lowlevel ~ date  + avg_airtmp + avg_soilwater + avg_precip + avg_precip * avg_soilwater, data = hydro_summer_10304_kbe, family = binomial)
summary(model_1)

model_2 <- glm(lowlevel ~ date  + avg_airtmp + avg_soilwater + avg_precip + avg_precip * avg_soilwater + avg_snowstorage + avg_glorad  + avg_infiltration, data = hydro_summer_10304_kbe, family = binomial)
summary(model_2)

# generalized linear mixed-effects model with random intercept for member (no interactions)
library(lme4)
glmer_1 <- glmer(lowlevel ~ precip + airtmp + glorad + relhum + soilwater + snowstorage + groundwaterdepth + infiltration + (1|member), data = hydro_winter_10304, family = binomial())
summary(glmer_1)
```

# First GAM Approaches ----
## GAM for lowlevel with "normal" drivers
```{r}
library(mgcv)
gam_1 <- gam(lowlevel ~ date + s(avg_airtmp) + s(avg_soilwater) + s(avg_precip) + avg_precip * avg_soilwater + avg_snowstorage + s(avg_glorad) + s(avg_infiltration), data = hydro_summer_10304_kbe, family = binomial(), method = "REML")
summary(gam_1)

gam_1_s_interac <- gam(lowlevel ~ s(precip) + s(airtmp) + s(glorad) + s(relhum) + s(soilwater) + s(snowstorage) + s(groundwaterdepth) + s(infiltration) + s(precip, by = soilwater), data = hydro_summer_10304, family = binomial(), method = "REML")
summary(gam_1_s_interac)
plot(gam_1_s_interac)

## gam mit glatten Effekten (penalized B-splines) und member als random Effekt
# 10304, winter
gam_randinter_member_10304_winter <- gam(lowlevel ~ date + s(precip, bs = "ps") + s(airtmp, bs = "ps") + s(glorad, bs = "ps") + s(relhum, bs = "ps") + s(soilwater, bs = "ps") + s(snowstorage, bs = "ps") + s(groundwaterdepth, bs = "ps") + s(infiltration, bs = "ps") + s(member, bs = "re"), data = hydro_winter_10304, family = binomial(), method = "REML")
summary(gam_randinter_member_10304_winter)
plot(gam_randinter_member_10304_winter)

# 10304, summer
gam_randinter_member_10304_summer <- gam(lowlevel ~ date + s(precip, bs = "ps") + s(airtmp, bs = "ps") + s(glorad, bs = "ps") + s(relhum, bs = "ps") + s(soilwater, bs = "ps") + s(snowstorage, bs = "ps") + s(groundwaterdepth, bs = "ps") + s(infiltration, bs = "ps") + s(member, bs = "re"), data = hydro_summer_10304, family = binomial(), method = "REML")
summary(gam_randinter_member_10304_summer)
plot(gam_randinter_member_10304_summer)

# 11502, winter
gam_randinter_member_11502_winter <- gam(lowlevel ~ date + s(precip, bs = "ps") + s(airtmp, bs = "ps") + s(glorad, bs = "ps") + s(relhum, bs = "ps") + s(soilwater, bs = "ps") + s(snowstorage, bs = "ps") + s(groundwaterdepth, bs = "ps") + s(infiltration, bs = "ps") + s(member, bs = "re"), data = hydro_winter_11502, family = binomial(), method = "REML")
summary(gam_randinter_member_11502_winter)
plot(gam_randinter_member_11502_winter)

# 11502, summer
gam_randinter_member_11502_summer <- gam(lowlevel ~ date + s(precip, bs = "ps") + s(airtmp, bs = "ps") + s(glorad, bs = "ps") + s(relhum, bs = "ps") + s(soilwater, bs = "ps") + s(snowstorage, bs = "ps") + s(groundwaterdepth, bs = "ps") + s(infiltration, bs = "ps") + s(member, bs = "re"), data = hydro_summer_11502, family = binomial(), method = "REML")
summary(gam_randinter_member_11502_summer)
plot(gam_randinter_member_11502_summer)

# 20203, winter
gam_randinter_member_20203_winter <- gam(lowlevel ~ date + s(precip, bs = "ps") + s(airtmp, bs = "ps") + s(glorad, bs = "ps") + s(relhum, bs = "ps") + s(soilwater, bs = "ps") + s(snowstorage, bs = "ps") + s(groundwaterdepth, bs = "ps") + s(infiltration, bs = "ps") + s(member, bs = "re"), data = hydro_winter_20203, family = binomial(), method = "REML")
summary(gam_randinter_member_20203_winter)
plot(gam_randinter_member_20203_winter)

# 20203, summer
gam_randinter_member_20203_summer <- gam(lowlevel ~ date + s(precip, bs = "ps") + s(airtmp, bs = "ps") + s(glorad, bs = "ps") + s(relhum, bs = "ps") + s(soilwater, bs = "ps") + s(snowstorage, bs = "ps") + s(groundwaterdepth, bs = "ps") + s(infiltration, bs = "ps") + s(member, bs = "re"), data = hydro_summer_20203, family = binomial(), method = "REML")
summary(gam_randinter_member_20203_summer)
plot(gam_randinter_member_20203_summer)

gam_randinter_member_ps <- gam(lowlevel ~ s(precip, bs = "ps") + s(airtmp, bs = "ps") + s(glorad, bs = "ps") + s(relhum, bs = "ps") + s(soilwater, bs = "ps") + s(snowstorage, bs = "ps") + s(groundwaterdepth, bs = "ps") + s(infiltration, bs = "ps") + s(precip, by = soilwater, bs = "ps") + s(member, bs = "re"), data = hydro_winter_10304, family = binomial(), method = "REML")
summary(gam_randinter_member_ps)
plot(gam_randinter_member_ps)


# generalized additive mixed model with ARMA-correlation-structure for member kbe
gamm_randinter_member_ps_corAR1 <- gamm(lowlevel ~ s(precip, bs = "ps") + s(airtmp, bs = "ps") + s(glorad, bs = "ps") + s(relhum, bs = "ps") + s(soilwater, bs = "ps") + s(snowstorage, bs = "ps") + s(groundwaterdepth, bs = "ps") + s(infiltration, bs = "ps") + s(precip, by = soilwater, bs = "ps") + s(member, bs = "re"), data = hydro_winter_10304, family = binomial(), correlation = corCAR1(form = ~ date | member))
# Ausprobieren: corARMA(p = 4, q = 1)
summary(gamm_randinter_member_ps_corAR1)


gam_2 <- gam(lowlevel ~ date + s(airtmp) + s(soilwater) + s(precip) + s(soilwater, precip), data = hydro_summer_10304_kbe, family = binomial())
summary(gam_2)


```

## GAM for lowlevel with only one avg/max driver
```{r}
library(mgcv)
## which drivers have a high explanatory power per water level and per season 
# 10304, winter
gam_kbe_avg_10304_winter_test <- gam(lowlevel ~ s(avg_snowstorage_drain, bs = "ps"), data = hydro_winter_10304_kbe, family = binomial(), method = "REML")
summary(gam_kbe_avg_10304_winter_test)
plot(gam_kbe_avg_10304_winter_test)
# viel Erklärungskraft alleine (>= 10%): avg_airtmp, avg_glorad, avg_soilwater, avg_snowstorage, groundwaterdepth, avg_infiltration, max_precip, avg_snowstorage_drain
# wenig Erklärungskraft alleine (< 10%): avg_precip, avg_relhum

# 10304, summer
gam_kbe_avg_10304_summer_test <- gam(lowlevel ~ s(avg_snowstorage_drain, bs = "ps"), data = hydro_summer_10304_kbe, family = binomial(), method = "REML")
summary(gam_kbe_avg_10304_summer_test)
# viel Erklärungskraft alleine: avg_precip, avg_glorad, avg_soilwater, avg_snowstorage, avg_infiltration, avg_snowstorage_drain
# wenig Erklärungskraft alleine: avg_airtmp, avg_relhum, groundwaterdepth, max_precip

# 11502, winter
gam_kbe_avg_11502_winter_test <- gam(lowlevel ~ s(avg_snowstorage_drain, bs = "ps"), data = hydro_winter_11502_kbe, family = binomial(), method = "REML")
summary(gam_kbe_avg_11502_winter_test)
# viel Erklärungskraft alleine: avg_precip, avg_glorad, avg_soilwater, avg_snowstorage, avg_infiltration, max_precip, avg_snowstorage_drain
# wenig Erklärungskraft alleine: avg_airtmp, avg_relhum, groundwaterdepth

# 11502, summer
gam_kbe_avg_11502_summer_test <- gam(lowlevel ~ s(avg_snowstorage_drain, bs = "ps"), data = hydro_summer_11502_kbe, family = binomial(), method = "REML")
summary(gam_kbe_avg_11502_summer_test)
# viel Erklärungskraft alleine: avg_precip, avg_soilwater, avg_snowstorage, groundwaterdepth, avg_infiltration
# wenig Erklärungskraft alleine: avg_airtmp, avg_glorad, avg_relhum, max_precip, avg_snowstorage_drain

# 20203, winter
gam_kbe_avg_20203_winter_test <- gam(lowlevel ~ s(avg_snowstorage_drain, bs = "ps"), data = hydro_winter_20203_kbe, family = binomial(), method = "REML")
summary(gam_kbe_avg_20203_winter_test)
# viel Erklärungskraft alleine: avg_precip, avg_relhum, avg_soilwater, avg_snowstorage, avg_infiltration, avg_snowstorage_drain
# wenig Erklärungskraft alleine: avg_airtmp, avg_glorad, groundwaterdepth, max_precip

# 20203, summer
gam_kbe_avg_20203_summer_test <- gam(lowlevel ~ s(avg_snowstorage, bs = "ps"), data = hydro_summer_20203_kbe, family = binomial(), method = "REML")
summary(gam_kbe_avg_20203_summer_test)
# viel Erklärungskraft alleine: avg_soilwater, avg_snowstorage_drain
# wenig Erklärungskraft alleine: avg_precip, avg_airtmp, avg_glorad, avg_relhum, avg_snowstorage, groundwaterdepth, avg_infiltration, max_precip
```

## GAM for lowlevel with avg/max drivers: full models without interactions etc.
```{r}
## gam mit avg/max drivers glatten Effekten (penalized B-splines) für member kbe
# 10304, winter
gam_kbe_avg_10304_winter <- gam(lowlevel ~ YY + s(avg_precip, bs = "ps") + s(avg_airtmp, bs = "ps") + s(avg_glorad, bs = "ps") + s(avg_relhum, bs = "ps") + s(avg_soilwater, bs = "ps") + s(avg_snowstorage, bs = "ps") + s(groundwaterdepth, bs = "ps") + s(avg_infiltration, bs = "ps") + s(max_precip, bs = "ps") + s(avg_snowstorage_drain, bs = "ps"), data = hydro_winter_10304_kbe, family = binomial(), method = "REML")
summary(gam_kbe_avg_10304_winter)
plot(gam_kbe_avg_10304_winter)

# 10304, summer
gam_kbe_avg_10304_summer <- gam(lowlevel ~ YY + s(avg_precip, bs = "ps") + s(avg_airtmp, bs = "ps") + s(avg_glorad, bs = "ps") + s(avg_relhum, bs = "ps") + s(avg_soilwater, bs = "ps") + s(avg_snowstorage, bs = "ps") + s(groundwaterdepth, bs = "ps") + s(avg_infiltration, bs = "ps") + s(max_precip, bs = "ps") + s(avg_snowstorage_drain, bs = "ps"), data = hydro_summer_10304_kbe, family = binomial(), method = "REML")
summary(gam_kbe_avg_10304_summer)
plot(gam_kbe_avg_10304_summer)

# 11502, winter
gam_kbe_avg_11502_winter <- gam(lowlevel ~ YY + s(avg_precip, bs = "ps") + s(avg_airtmp, bs = "ps") + s(avg_glorad, bs = "ps") + s(avg_relhum, bs = "ps") + s(avg_soilwater, bs = "ps") + s(avg_snowstorage, bs = "ps") + s(groundwaterdepth, bs = "ps") + s(avg_infiltration, bs = "ps") + s(max_precip, bs = "ps") + s(avg_snowstorage_drain, bs = "ps"), data = hydro_winter_11502_kbe, family = binomial(), method = "REML")
summary(gam_kbe_avg_11502_winter)
plot(gam_kbe_avg_11502_winter)

# 11502, summer
gam_kbe_avg_11502_summer <- gam(lowlevel ~ YY + s(avg_precip, bs = "ps") + s(avg_airtmp, bs = "ps") + s(avg_glorad, bs = "ps") + s(avg_relhum, bs = "ps") + s(avg_soilwater, bs = "ps") + s(avg_snowstorage, bs = "ps") + s(groundwaterdepth, bs = "ps") + s(avg_infiltration, bs = "ps") + s(max_precip, bs = "ps") + s(avg_snowstorage_drain, bs = "ps"), data = hydro_summer_11502_kbe, family = binomial(), method = "REML")
summary(gam_kbe_avg_11502_summer)
plot(gam_kbe_avg_11502_summer)

# 20203, winter
gam_kbe_avg_20203_winter <- gam(lowlevel ~ YY + s(avg_precip, bs = "ps") + s(avg_airtmp, bs = "ps") + s(avg_glorad, bs = "ps") + s(avg_relhum, bs = "ps") + s(avg_soilwater, bs = "ps") + s(avg_snowstorage, bs = "ps") + s(groundwaterdepth, bs = "ps") + s(avg_infiltration, bs = "ps") + s(max_precip, bs = "ps") + s(avg_snowstorage_drain, bs = "ps"), data = hydro_winter_20203_kbe, family = binomial(), method = "REML")
summary(gam_kbe_avg_20203_winter)
plot(gam_kbe_avg_20203_winter)

# 20203, summer
gam_kbe_avg_20203_summer <- gam(lowlevel ~ YY + s(avg_precip, bs = "ps") + s(avg_airtmp, bs = "ps") + s(avg_glorad, bs = "ps") + s(avg_relhum, bs = "ps") + s(avg_soilwater, bs = "ps") + s(avg_snowstorage, bs = "ps") + s(groundwaterdepth, bs = "ps") + s(avg_infiltration, bs = "ps") + s(max_precip, bs = "ps") + s(avg_snowstorage_drain, bs = "ps"), data = hydro_summer_20203_kbe, family = binomial(), method = "REML")
summary(gam_kbe_avg_20203_summer)
plot(gam_kbe_avg_20203_summer)
```

```{r}
#stepwise AIC-variable-selection for winter 10304
AIC(gam_kbe_avg_10304_winter)

gam_kbe_avg_10304_winter_wo_avg_relhum <- gam(lowlevel ~ YY + s(avg_precip, bs = "ps") + s(avg_airtmp, bs = "ps") + s(avg_glorad, bs = "ps") + s(avg_soilwater, bs = "ps") + s(avg_snowstorage, bs = "ps") + s(groundwaterdepth, bs = "ps") + s(avg_infiltration, bs = "ps") + s(max_precip, bs = "ps"), data = hydro_winter_10304_kbe, family = binomial(), method = "REML")
AIC(gam_kbe_avg_10304_winter_wo_avg_relhum)
BIC(gam_kbe_avg_10304_winter_wo_avg_relhum)

gam_kbe_avg_10304_winter_wo_avg_relhum_wo_max_precip <- gam(lowlevel ~ YY + s(avg_precip, bs = "ps") + s(avg_airtmp, bs = "ps") + s(avg_soilwater, bs = "ps") + s(avg_snowstorage, bs = "ps") + s(avg_glorad, bs = "ps") + s(groundwaterdepth, bs = "ps") + s(avg_infiltration, bs = "ps"), data = hydro_winter_10304_kbe, family = binomial(), method = "REML")
summary(gam_kbe_avg_10304_winter_wo_avg_relhum_wo_max_precip)
AIC(gam_kbe_avg_10304_winter_wo_avg_relhum_wo_max_precip)

gam_kbe_avg_10304_winter_wo_avg_relhum_wo_max_precip_wo_avg_airtmp <- gam(lowlevel ~ YY + s(avg_precip, bs = "ps") + s(avg_glorad, bs = "ps") + s(avg_soilwater, bs = "ps") + s(avg_snowstorage, bs = "ps") + s(groundwaterdepth, bs = "ps") + s(avg_infiltration, bs = "ps"), data = hydro_winter_10304_kbe, family = binomial(), method = "REML")
summary(gam_kbe_avg_10304_winter_wo_avg_relhum_wo_max_precip_wo_avg_airtmp)
AIC(gam_kbe_avg_10304_winter_wo_avg_relhum_wo_max_precip_wo_avg_airtmp)

gam_kbe_avg_10304_winter_wo_avg_relhum_wo_max_precip_wo_avg_airtmp_wo_avg_precip <- gam(lowlevel ~ YY + s(avg_glorad, bs = "ps") + s(avg_soilwater, bs = "ps") + s(avg_snowstorage, bs = "ps") + s(groundwaterdepth, bs = "ps") + s(avg_infiltration, bs = "ps"), data = hydro_winter_10304_kbe, family = binomial(), method = "REML")
summary(gam_kbe_avg_10304_winter_wo_avg_relhum_wo_max_precip_wo_avg_airtmp_wo_avg_precip)
plot(gam_kbe_avg_10304_winter_wo_avg_relhum_wo_max_precip_wo_avg_airtmp_wo_avg_precip)
BIC(gam_kbe_avg_10304_winter_wo_avg_relhum_wo_max_precip_wo_avg_airtmp_wo_avg_precip)

#stepwise AIC-variable-selection for summer 10304
summary(gam_kbe_avg_10304_summer)

gam_kbe_avg_10304_summer_wo_avg_relhum <- gam(lowlevel ~ YY + s(avg_precip, bs = "ps") + s(avg_airtmp, bs = "ps") + s(avg_glorad, bs = "ps") + s(avg_soilwater, bs = "ps") + s(avg_snowstorage, bs = "ps") + s(groundwaterdepth, bs = "ps") + s(avg_infiltration, bs = "ps") + s(max_precip, bs = "ps"), data = hydro_summer_10304_kbe, family = binomial(), method = "REML")
summary(gam_kbe_avg_10304_summer_wo_avg_relhum)

gam_kbe_avg_10304_summer_wo_avg_relhum_wo_max_precip <- gam(lowlevel ~ YY + s(avg_precip, bs = "ps") + s(avg_airtmp, bs = "ps") + s(avg_glorad, bs = "ps") + s(avg_soilwater, bs = "ps") + s(avg_snowstorage, bs = "ps") + s(groundwaterdepth, bs = "ps") + s(avg_infiltration, bs = "ps"), data = hydro_summer_10304_kbe, family = binomial(), method = "REML")
summary(gam_kbe_avg_10304_summer_wo_avg_relhum_wo_max_precip)

gam_kbe_avg_10304_summer_wo_avg_relhum_wo_max_precip_wo_avg_precip <- gam(lowlevel ~ YY + s(avg_airtmp, bs = "ps") + s(avg_glorad, bs = "ps") + s(avg_soilwater, bs = "ps") + s(avg_snowstorage, bs = "ps") + s(groundwaterdepth, bs = "ps") + s(avg_infiltration, bs = "ps"), data = hydro_summer_10304_kbe, family = binomial(), method = "REML")
summary(gam_kbe_avg_10304_summer_wo_avg_relhum_wo_max_precip_wo_avg_precip)

gam_kbe_avg_10304_summer_wo_avg_relhum_wo_max_precip_wo_avg_precip_wo_snowstorage <- gam(lowlevel ~ YY + s(avg_airtmp, bs = "ps") + s(avg_glorad, bs = "ps") + s(avg_soilwater, bs = "ps") + s(groundwaterdepth, bs = "ps") + s(avg_infiltration, bs = "ps"), data = hydro_summer_10304_kbe, family = binomial(), method = "REML")
summary(gam_kbe_avg_10304_summer_wo_avg_relhum_wo_max_precip_wo_avg_precip_wo_snowstorage)
plot(gam_kbe_avg_10304_summer_wo_avg_relhum_wo_max_precip_wo_avg_precip_wo_snowstorage)

#stepwise AIC-variable-selection for winter 11502
summary(gam_kbe_avg_11502_winter)

gam_kbe_avg_11502_winter_wo_avg_relhum <- gam(lowlevel ~ YY + s(avg_precip, bs = "ps") + s(avg_airtmp, bs = "ps") + s(avg_glorad, bs = "ps") + s(avg_soilwater, bs = "ps") + s(avg_snowstorage, bs = "ps") + s(groundwaterdepth, bs = "ps") + s(avg_infiltration, bs = "ps") + s(max_precip, bs = "ps"), data = hydro_winter_11502_kbe, family = binomial(), method = "REML")
summary(gam_kbe_avg_11502_winter_wo_avg_relhum)

gam_kbe_avg_11502_winter_wo_avg_relhum_wo_avg_glorad <- gam(lowlevel ~ YY + s(avg_precip, bs = "ps") + s(avg_airtmp, bs = "ps") + s(avg_soilwater, bs = "ps") + s(avg_snowstorage, bs = "ps") + s(groundwaterdepth, bs = "ps") + s(avg_infiltration, bs = "ps") + s(max_precip, bs = "ps"), data = hydro_winter_11502_kbe, family = binomial(), method = "REML")
summary(gam_kbe_avg_11502_winter_wo_avg_relhum_wo_avg_glorad)

gam_kbe_avg_11502_winter_wo_avg_relhum_wo_avg_glorad_wo_avg_precip <- gam(lowlevel ~ YY + s(avg_airtmp, bs = "ps") + s(avg_soilwater, bs = "ps") + s(avg_snowstorage, bs = "ps") + s(groundwaterdepth, bs = "ps") + s(avg_infiltration, bs = "ps") + s(max_precip, bs = "ps"), data = hydro_winter_11502_kbe, family = binomial(), method = "REML")
summary(gam_kbe_avg_11502_winter_wo_avg_relhum_wo_avg_glorad_wo_avg_precip)

gam_kbe_avg_11502_winter_wo_avg_relhum_wo_avg_glorad_wo_avg_precip_wo_avg_airtmp <- gam(lowlevel ~ YY + s(avg_soilwater, bs = "ps") + s(avg_snowstorage, bs = "ps") + s(groundwaterdepth, bs = "ps") + s(avg_infiltration, bs = "ps") + s(max_precip, bs = "ps"), data = hydro_winter_11502_kbe, family = binomial(), method = "REML")
summary(gam_kbe_avg_11502_winter_wo_avg_relhum_wo_avg_glorad_wo_avg_precip_wo_avg_airtmp)

#stepwise AIC-variable-selection for summer 11502
summary(gam_kbe_avg_11502_summer)

gam_kbe_avg_11502_summer_wo_avg_precip <- gam(lowlevel ~ YY + s(avg_airtmp, bs = "ps") + s(avg_glorad, bs = "ps") + s(avg_relhum, bs = "ps") + s(avg_soilwater, bs = "ps") + s(avg_snowstorage, bs = "ps") + s(groundwaterdepth, bs = "ps") + s(avg_infiltration, bs = "ps") + s(max_precip, bs = "ps"), data = hydro_summer_11502_kbe, family = binomial(), method = "REML")
summary(gam_kbe_avg_11502_summer_wo_avg_precip)

gam_kbe_avg_11502_summer_wo_avg_precip_wo_max_precip <- gam(lowlevel ~ YY + s(avg_airtmp, bs = "ps") + s(avg_glorad, bs = "ps") + s(avg_relhum, bs = "ps") + s(avg_soilwater, bs = "ps") + s(avg_snowstorage, bs = "ps") + s(groundwaterdepth, bs = "ps") + s(avg_infiltration, bs = "ps"), data = hydro_summer_11502_kbe, family = binomial(), method = "REML")
summary(gam_kbe_avg_11502_summer_wo_avg_precip_wo_max_precip)

gam_kbe_avg_11502_summer_wo_avg_precip_wo_max_precip_wo_avg_relhum <- gam(lowlevel ~ YY + s(avg_airtmp, bs = "ps") + s(avg_glorad, bs = "ps") + s(avg_soilwater, bs = "ps") + s(avg_snowstorage, bs = "ps") + s(groundwaterdepth, bs = "ps") + s(avg_infiltration, bs = "ps"), data = hydro_summer_11502_kbe, family = binomial(), method = "REML")
summary(gam_kbe_avg_11502_summer_wo_avg_precip_wo_max_precip_wo_avg_relhum)

gam_kbe_avg_11502_summer_wo_avg_precip_wo_max_precip_wo_avg_relhum_wo_avg_airtmp <- gam(lowlevel ~ YY + s(avg_glorad, bs = "ps") + s(avg_soilwater, bs = "ps") + s(avg_snowstorage, bs = "ps") + s(groundwaterdepth, bs = "ps") + s(avg_infiltration, bs = "ps"), data = hydro_summer_11502_kbe, family = binomial(), method = "REML")
plot(gam_kbe_avg_11502_summer_wo_avg_precip_wo_max_precip_wo_avg_relhum_wo_avg_airtmp)
summary(gam_kbe_avg_11502_summer_wo_avg_precip_wo_max_precip_wo_avg_relhum_wo_avg_airtmp)

#stepwise AIC-variable-selection for winter 20203
summary(gam_kbe_avg_20203_winter)

gam_kbe_avg_20203_winter_wo_avg_precip <- gam(lowlevel ~ YY + s(avg_airtmp, bs = "ps") + s(avg_glorad, bs = "ps") + s(avg_relhum, bs = "ps") + s(avg_soilwater, bs = "ps") + s(avg_snowstorage, bs = "ps") + s(groundwaterdepth, bs = "ps") + s(avg_infiltration, bs = "ps") + s(max_precip, bs = "ps"), data = hydro_winter_20203_kbe, family = binomial(), method = "REML")
summary(gam_kbe_avg_20203_winter_wo_avg_precip)

gam_kbe_avg_20203_winter_wo_avg_precip_wo_avg_airtmp <- gam(lowlevel ~ YY + s(avg_glorad, bs = "ps") + s(avg_relhum, bs = "ps") + s(avg_soilwater, bs = "ps") + s(avg_snowstorage, bs = "ps") + s(groundwaterdepth, bs = "ps") + s(avg_infiltration, bs = "ps") + s(max_precip, bs = "ps"), data = hydro_winter_20203_kbe, family = binomial(), method = "REML")
summary(gam_kbe_avg_20203_winter_wo_avg_precip_wo_avg_airtmp)

gam_kbe_avg_20203_winter_wo_avg_precip_wo_avg_airtmp_wo_avg_glorad <- gam(lowlevel ~ YY + s(avg_relhum, bs = "ps") + s(avg_soilwater, bs = "ps") + s(avg_snowstorage, bs = "ps") + s(groundwaterdepth, bs = "ps") + s(avg_infiltration, bs = "ps") + s(max_precip, bs = "ps"), data = hydro_winter_20203_kbe, family = binomial(), method = "REML")
summary(gam_kbe_avg_20203_winter_wo_avg_precip_wo_avg_airtmp_wo_avg_glorad)

gam_kbe_avg_20203_winter_wo_avg_precip_wo_avg_airtmp_wo_avg_glorad_wo_avg_snowstorage <- gam(lowlevel ~ YY + s(avg_relhum, bs = "ps") + s(avg_soilwater, bs = "ps") + s(groundwaterdepth, bs = "ps") + s(avg_infiltration, bs = "ps") + s(max_precip, bs = "ps"), data = hydro_winter_20203_kbe, family = binomial(), method = "REML")
summary(gam_kbe_avg_20203_winter_wo_avg_precip_wo_avg_airtmp_wo_avg_glorad_wo_avg_snowstorage)

gam_kbe_avg_20203_winter_wo_avg_precip_wo_avg_airtmp_wo_avg_glorad_wo_avg_snowstorage_wo_groundwaterdepth <- gam(lowlevel ~ YY + s(avg_relhum, bs = "ps") + s(avg_soilwater, bs = "ps") + s(avg_infiltration, bs = "ps") + s(max_precip, bs = "ps"), data = hydro_winter_20203_kbe, family = binomial(), method = "REML")
summary(gam_kbe_avg_20203_winter_wo_avg_precip_wo_avg_airtmp_wo_avg_glorad_wo_avg_snowstorage_wo_groundwaterdepth)

#stepwise AIC-variable-selection for winter 20203
BIC(gam_kbe_avg_20203_summer)

gam_kbe_avg_20203_summer_wo_avg_glorad <- gam(lowlevel ~ YY + s(avg_precip, bs = "ps") + s(avg_airtmp, bs = "ps") + s(avg_relhum, bs = "ps") + s(avg_soilwater, bs = "ps") + s(avg_snowstorage, bs = "ps") + s(groundwaterdepth, bs = "ps") + s(avg_infiltration, bs = "ps") + s(max_precip, bs = "ps"), data = hydro_summer_20203_kbe, family = binomial(), method = "REML")
summary(gam_kbe_avg_20203_summer_wo_avg_glorad)

gam_kbe_avg_20203_summer_wo_avg_glorad_wo_avg_airtmp <- gam(lowlevel ~ YY + s(avg_precip, bs = "ps") + s(avg_relhum, bs = "ps") + s(avg_soilwater, bs = "ps") + s(avg_snowstorage, bs = "ps") + s(groundwaterdepth, bs = "ps") + s(avg_infiltration, bs = "ps") + s(max_precip, bs = "ps"), data = hydro_summer_20203_kbe, family = binomial(), method = "REML")
summary(gam_kbe_avg_20203_summer_wo_avg_glorad_wo_avg_airtmp)

gam_kbe_avg_20203_summer_wo_avg_glorad_wo_avg_airtmp_wo_max_precip <- gam(lowlevel ~ YY + s(avg_precip, bs = "ps") + s(avg_relhum, bs = "ps") + s(avg_soilwater, bs = "ps") + s(avg_snowstorage, bs = "ps") + s(groundwaterdepth, bs = "ps") + s(avg_infiltration, bs = "ps"), data = hydro_summer_20203_kbe, family = binomial(), method = "REML")
BIC(gam_kbe_avg_20203_summer_wo_avg_glorad_wo_avg_airtmp_wo_max_precip)

gam_kbe_avg_20203_summer_wo_avg_glorad_wo_avg_airtmp_wo_max_precip_wo_avg_precip <- gam(lowlevel ~ YY + s(avg_relhum, bs = "ps") + s(avg_soilwater, bs = "ps") + s(avg_snowstorage, bs = "ps") + s(groundwaterdepth, bs = "ps") + s(avg_infiltration, bs = "ps"), data = hydro_summer_20203_kbe, family = binomial(), method = "REML")
summary(gam_kbe_avg_20203_summer_wo_avg_glorad_wo_avg_airtmp_wo_max_precip_wo_avg_precip)
```

## Variable selection of full models
```{r}
# 10304, winter
# NAs durch Mittelwert der Spalte imputieren
for (i in 1:length(names(hydro_winter_10304_kbe))) {
  if (is.numeric(hydro_winter_10304_kbe[, i])) {
    for (k in 1:nrow(hydro_winter_10304_kbe)) {
      if (is.na(hydro_winter_10304_kbe[k, i])) {
        hydro_winter_10304_kbe[k, i] <- mean(hydro_winter_10304_kbe[, i], na.rm = TRUE)
      }
    }
  }
}
hydro_winter_10304_kbe$max_precip <- zoo::rollmax(hydro_winter_10304_kbe$precip, 7, fill = mean(hydro_winter_10304_kbe$precip))

gam_kbe_avg_10304_winter <- gam(lowlevel ~ YY + s(avg_precip, bs = "ps") + s(avg_airtmp, bs = "ps") + s(avg_glorad, bs = "ps") + s(avg_relhum, bs = "ps") + s(avg_soilwater, bs = "ps") + s(avg_snowstorage, bs = "ps") + s(groundwaterdepth, bs = "ps") + s(avg_infiltration, bs = "ps") + s(max_precip, bs = "ps"), data = hydro_winter_10304_kbe, family = binomial(), method = "REML")
# Null Space
gam_kbe_avg_10304_winter_null_space <- gam(lowlevel ~ YY + s(avg_precip, bs = "ps") + s(avg_airtmp, bs = "ps") + s(avg_glorad, bs = "ps") + s(avg_relhum, bs = "ps") + s(avg_soilwater, bs = "ps") + s(avg_snowstorage, bs = "ps") + s(groundwaterdepth, bs = "ps") + s(avg_infiltration, bs = "ps") + s(max_precip, bs = "ps"), data = hydro_winter_10304_kbe, family = binomial(), method = "REML", select = T)
summary(gam_kbe_avg_10304_winter_null_space)

# Shrinkage Smoothers
gam_kbe_avg_10304_winter_shrinkage <- gam(lowlevel ~ YY + s(avg_precip, bs = "cs") + s(avg_airtmp, bs = "cs") + s(avg_glorad, bs = "cs") + s(avg_relhum, bs = "cs") + s(avg_soilwater, bs = "cs") + s(avg_snowstorage, bs = "cs") + s(groundwaterdepth, bs = "cs") + s(avg_infiltration, bs = "cs") + s(max_precip, bs = "cs"), data = hydro_winter_10304_kbe, family = binomial(), method = "REML")
summary(gam_kbe_avg_10304_winter_shrinkage)



gam_kbe_avg_10304_winter_select <- gam(lowlevel ~ YY + s(avg_precip, bs = "ps") + s(avg_airtmp, bs = "ps") + s(avg_glorad, bs = "ps") + s(avg_relhum, bs = "ps") + s(avg_soilwater, bs = "ps") + s(avg_snowstorage, bs = "ps") + s(groundwaterdepth, bs = "ps") + s(avg_infiltration, bs = "ps") + s(max_precip, bs = "ps"), data = hydro_winter_10304_kbe, family = binomial(), method = "REML", select = TRUE)
summary(gam_kbe_avg_10304_winter_select)
plot(gam_kbe_avg_10304_winter_select)

gam_kbe_avg_10304_winter_ts <- gam(lowlevel ~ YY + s(avg_precip, bs = "ts") + s(avg_airtmp, bs = "ts") + s(avg_glorad, bs = "ts") + s(avg_relhum, bs = "ts") + s(avg_soilwater, bs = "ts") + s(avg_snowstorage, bs = "ts") + s(groundwaterdepth, bs = "ts") + s(avg_infiltration, bs = "ts") + s(max_precip, bs = "ts"), data = hydro_winter_10304_kbe, family = binomial(), method = "REML")
summary(gam_kbe_avg_10304_winter_ts)
plot(gam_kbe_avg_10304_winter_ts)

# 10304, summer
gam_kbe_avg_10304_summer_select <- gam(lowlevel ~ YY + s(avg_precip) + s(avg_airtmp) + s(avg_glorad) + s(avg_relhum) + s(avg_soilwater) + s(avg_snowstorage) + s(groundwaterdepth) + s(avg_infiltration) + s(max_precip), data = hydro_summer_10304_kbe, family = binomial(), method = "REML", select = TRUE)
summary(gam_kbe_avg_10304_summer_select)
plot(gam_kbe_avg_10304_summer_select)

# 11502, winter
gam_kbe_avg_11502_winter_select <- gam(lowlevel ~ YY + s(avg_precip, bs = "ps") + s(avg_airtmp, bs = "ps") + s(avg_glorad, bs = "ps") + s(avg_relhum, bs = "ps") + s(avg_soilwater, bs = "ps") + s(avg_snowstorage, bs = "ps") + s(groundwaterdepth, bs = "ps") + s(avg_infiltration, bs = "ps") + s(max_precip, bs = "ps"), data = hydro_winter_11502_kbe, family = binomial(), method = "REML", select = TRUE)
summary(gam_kbe_avg_11502_winter_select)
plot(gam_kbe_avg_11502_winter_select)
```
## GAM for lowlevel_first with avg drivers
```{r}
library(mgcv)

# 10304, winter
gam_kbe_avg_first_10304_winter <- gam(lowlevel_first ~ YY + s(avg_precip, bs = "ps") + s(avg_airtmp, bs = "ps") + s(avg_glorad, bs = "ps") + s(avg_relhum, bs = "ps") + s(avg_soilwater, bs = "ps") + s(avg_snowstorage, bs = "ps") + s(groundwaterdepth, bs = "ps") + s(avg_infiltration, bs = "ps") + s(max_precip, bs = "ps"), data = hydro_winter_10304_kbe, family = binomial(), method = "REML")
summary(gam_kbe_avg_first_10304_winter)
plot(exp(predict(gam_kbe_avg_first_10304_winter)), residuals(gam_kbe_avg_first_10304_winter))
plot(gam_kbe_avg_first_10304_winter)

# 10304, summer
gam_kbe_avg_first_10304_summer <- gam(lowlevel_first ~ YY + s(avg_precip, bs = "ps") + s(avg_airtmp, bs = "ps") + s(avg_glorad, bs = "ps") + s(avg_relhum, bs = "ps") + s(avg_soilwater, bs = "ps") + s(avg_snowstorage, bs = "ps") + s(groundwaterdepth, bs = "ps") + s(avg_infiltration, bs = "ps") + s(max_precip, bs = "ps"), data = hydro_summer_10304_kbe, family = binomial(), method = "REML")
summary(gam_kbe_avg_first_10304_summer)
plot(exp(predict(gam_kbe_avg_first_10304_summer)), residuals(gam_kbe_avg_first_10304_summer))
plot(gam_kbe_avg_first_10304_summer)

# 11502, winter
gam_kbe_avg_first_11502_winter <- gam(lowlevel_first ~ YY + s(avg_precip, bs = "ps") + s(avg_airtmp, bs = "ps") + s(avg_glorad, bs = "ps") + s(avg_relhum, bs = "ps") + s(avg_soilwater, bs = "ps") + s(avg_snowstorage, bs = "ps") + s(groundwaterdepth, bs = "ps") + s(avg_infiltration, bs = "ps") + s(max_precip, bs = "ps"), data = hydro_winter_11502_kbe, family = binomial(), method = "REML")
summary(gam_kbe_avg_first_11502_winter)
plot(exp(predict(gam_kbe_avg_first_11502_winter)), residuals(gam_kbe_avg_first_11502_winter))
plot(gam_kbe_avg_first_11502_winter)

# 11502, summer
gam_kbe_avg_first_11502_summer <- gam(lowlevel_first ~ YY + s(avg_precip, bs = "ps") + s(avg_airtmp, bs = "ps") + s(avg_glorad, bs = "ps") + s(avg_relhum, bs = "ps") + s(avg_soilwater, bs = "ps") + s(avg_snowstorage, bs = "ps") + s(groundwaterdepth, bs = "ps") + s(avg_infiltration, bs = "ps") + s(max_precip, bs = "ps"), data = hydro_summer_11502_kbe, family = binomial(), method = "REML")
summary(gam_kbe_avg_first_11502_summer)
plot(exp(predict(gam_kbe_avg_first_11502_summer)), residuals(gam_kbe_avg_first_11502_summer))
plot(gam_kbe_avg_first_11502_summer)

# 20203, winter
gam_kbe_avg_first_20203_winter <- gam(lowlevel_first ~ YY + s(avg_precip, bs = "ps") + s(avg_airtmp, bs = "ps") + s(avg_glorad, bs = "ps") + s(avg_relhum, bs = "ps") + s(avg_soilwater, bs = "ps") + s(avg_snowstorage, bs = "ps") + s(groundwaterdepth, bs = "ps") + s(avg_infiltration, bs = "ps") + s(max_precip, bs = "ps"), data = hydro_winter_20203_kbe, family = binomial(), method = "REML")
summary(gam_kbe_avg_first_20203_winter)
plot(exp(predict(gam_kbe_avg_first_20203_winter)), residuals(gam_kbe_avg_first_20203_winter))
plot(gam_kbe_avg_first_20203_winter)

# 20203, summer
gam_kbe_avg_first_20203_summer <- gam(lowlevel_first ~ YY + s(avg_precip, bs = "ps") + s(avg_airtmp, bs = "ps") + s(avg_glorad, bs = "ps") + s(avg_relhum, bs = "ps") + s(avg_soilwater, bs = "ps") + s(avg_snowstorage, bs = "ps") + s(groundwaterdepth, bs = "ps") + s(avg_infiltration, bs = "ps") + s(max_precip, bs = "ps"), data = hydro_summer_20203_kbe, family = binomial(), method = "REML")
summary(gam_kbe_avg_first_20203_summer)
plot(exp(predict(gam_kbe_avg_first_20203_summer)), residuals(gam_kbe_avg_first_20203_summer))
plot(gam_kbe_avg_first_20203_summer)

# -> Variable lowlevel_first sehr unbalanciert -> schwer zu modellieren

```

## P Splines
```{r}
gam_1_ps <- gam(lowlevel ~ date  + s(avg_airtmp, bs = "ps") + s(avg_soilwater, bs = "ps") + s(avg_precip, bs = "ps") + avg_precip * avg_soilwater + avg_snowstorage + s(avg_glorad, bs = "ps")  + s(avg_infiltration, bs = "ps"), data = hydro_summer_10304_kbe, family = binomial(), method = "REML")
summary(gam_1_ps)

gam_2_ps <- gam(lowlevel ~ date + s(airtmp, bs = "ps") + s(soilwater, bs = "ps") + s(precip, bs = "ps") + s(soilwater, by = precip, bs = "ps"), data = hydro_summer_10304_kbe, family = binomial())

summary(gam_2_ps)
```

## Case Study on Catchment 10304 Isar Mittenwald
```{r}
glm_all_10304 <- glm(formula = lowlevel ~ YY + avg_precip + avg_airtmp+ avg_glorad + avg_relhum + avg_soilwater + avg_snowstorage + groundwaterdepth + avg_infiltration, family = binomial(), data = hydro_winter_10304_kbe)

stepAIC(object = glm_all_10304, direction = "both",
k = log(nrow(hydro_winter_10304_kbe)), trace = FALSE)
```

```{r}
gam_1_10304 <- gam(lowlevel ~ YY + avg_precip + avg_airtmp+ avg_glorad + avg_relhum + avg_soilwater + avg_snowstorage + groundwaterdepth + avg_infiltration, data = hydro_winter_10304_kbe, family = binomial())
gam_2_10304 <- gam(lowlevel ~ YY + avg_precip + avg_airtmp + avg_glorad + avg_relhum + avg_soilwater + avg30_snowstorage + groundwaterdepth + avg_infiltration, data = hydro_winter_10304_kbe, family = binomial())
gam_3_10304 <- gam(lowlevel ~ YY + s(avg_precip, bs = "ps") + s(avg_airtmp, bs = "ps") + avg_glorad + avg_relhum + avg_soilwater + avg_snowstorage + groundwaterdepth + avg_infiltration, data = hydro_winter_10304_kbe, family = binomial())
gam_4_10304 <- gam(lowlevel ~ YY + s(avg_precip, bs = "ps") + s(avg_airtmp, bs = "ps") + avg_glorad + avg_relhum + avg_soilwater + s(avg30_snowstorage, bs = "ps") + groundwaterdepth + avg_infiltration, data = hydro_winter_10304_kbe, family = binomial())
gam_5_10304 <- gam(lowlevel ~ YY + s(avg_precip, bs = "ps") + s(avg_airtmp, bs = "ps") + avg_glorad + avg_relhum + s(avg_soilwater, bs = "ps") + s(avg30_snowstorage, bs = "ps") + groundwaterdepth + avg_infiltration, data = hydro_winter_10304_kbe, family = binomial())
gam_6_10304 <- gam(lowlevel ~ YY + s(avg_precip, bs = "ps") + avg_airtmp + avg_glorad + avg_relhum + s(avg_soilwater, bs = "ps") + s(avg30_snowstorage, bs = "ps") + groundwaterdepth + avg_infiltration, data = hydro_winter_10304_kbe, family = binomial())
gam_7_10304 <- gam(lowlevel ~ YY + s(avg_precip, bs = "ps") + avg_airtmp + avg_relhum + s(avg_soilwater, bs = "ps") + s(avg30_snowstorage, bs = "ps") + groundwaterdepth + avg_infiltration, data = hydro_winter_10304_kbe, family = binomial())
gam_8_10304 <- gam(lowlevel ~ s(avg_precip, bs = "ps") + s(avg_airtmp, bs = "ps") + s(avg_glorad, bs = "ps") + s(avg_relhum, bs = "ps") + s(avg_soilwater, bs = "ps") + s(avg_snowstorage, bs = "ps") + s(groundwaterdepth, bs = "ps") + s(avg_infiltration, bs = "ps"), data = hydro_winter_10304_kbe, family = binomial())

gam_9_10304 <- gam(lowlevel ~ s(avg_precip, bs = "ps") + s(avg_glorad, bs = "ps") + s(avg_relhum, bs = "ps") + s(avg_soilwater, bs = "ps") + s(avg_snowstorage, bs = "ps") + s(groundwaterdepth, bs = "ps") + s(avg_infiltration, bs = "ps"), data = hydro_winter_10304_kbe, family = binomial())

gam_10_10304 <- gam(lowlevel ~ s(avg_precip, bs = "ps") + s(avg_glorad, bs = "ps") + s(avg_relhum, bs = "ps") + s(avg_soilwater, bs = "ps") + s(avg30_snowstorage, bs = "ps") + s(groundwaterdepth, bs = "ps") + s(avg_infiltration, bs = "ps"), data = hydro_winter_10304_kbe, family = binomial())

gam_11_10304 <- gam(lowlevel ~ s(avg_precip, bs = "ps") + s(avg_glorad, bs = "ps") + s(avg_relhum, bs = "ps") + s(avg_soilwater, bs = "ps") + s(avg30_snowstorage, bs = "ps") + s(groundwaterdepth, bs = "ps") + s(avg_infiltration, bs = "ps") + s(avg_soilwater, by = avg_precip) , data = hydro_winter_10304_kbe, family = binomial())

gam_12_10304 <- gam(lowlevel ~ s(avg_precip, bs = "ps") + s(avg_glorad, bs = "ps") + s(avg_relhum, bs = "ps") + s(avg_soilwater, bs = "ps") + s(avg30_snowstorage, bs = "ps") + s(groundwaterdepth, bs = "ps") + s(avg_infiltration, bs = "ps") + ti(avg_precip, airtmp), data = hydro_winter_10304_kbe, family = binomial())

gam_13_10304 <- gam(lowlevel ~ s(avg_precip, bs = "ps") + s(avg_glorad, bs = "ps") + s(avg_relhum, bs = "ps") + s(avg_soilwater, bs = "ps") + s(avg30_snowstorage, bs = "ps") + s(groundwaterdepth, bs = "ps") + s(avg_infiltration, bs = "ps") + s(YY, bs = "re"), data = hydro_winter_10304_kbe, family = binomial())

model_10304 <- list(model_1 = gam_1_10304, model_2 = gam_2_10304, model_3 = gam_3_10304, model_4 = gam_4_10304, model_5 = gam_5_10304, model_6 = gam_6_10304, model_7 = gam_7_10304, model_8 = gam_8_10304, model_9 = gam_9_10304, model_10 = gam_10_10304, model_11 = gam_11_10304, model_12 = gam_12_10304, model_13 = gam_13_10304)

lapply(model_10304, AIC)
```



# Autoregressive Approaches ----
## ARIMA-Models
```{r}
library(forecast)

model_Arima <- Arima(hydro_winter_10304_kbe$lowlevel,xreg = as.matrix(subset(hydro_winter_10304_kbe, select = c(precip, airtmp, glorad, relhum, soilwater, snowstorage, groundwaterdepth, infiltration))), order = c(4, 0, 7))
summary(model_Arima)

acf(resid(model_Arima))
```

## Lag DF Visualization
```{r}
lag_df <- function(dat, lag) {
  data.frame(unclass(dat[c(rep(NA, lag), 1:(nrow(dat)-lag)),]))
}

ggplot(aes(x = (lag(hydro_summer_10304_kbe$airtmp, 3)), y = hydro_summer_10304_kbe$drainage), data = hydro_summer_10304_kbe)+
  geom_point(alpha = 0.2)

ggplot(aes(x = (lag(hydro_summer_10304_kbe$soilwater, 2)), y = hydro_summer_10304_kbe$drainage), data = hydro_summer_10304_kbe)+
  geom_point(alpha = 0.2)

ggplot(aes(x = (lag(hydro_summer_10304_kbe$precip, 2)), y = hydro_summer_10304_kbe$drainage), data = hydro_summer_10304_kbe)+
  geom_point(alpha = 0.2)

ggplot(aes(x = (lag(hydro_summer_10304_kbe$avg_precip, 0)), y = hydro_summer_10304_kbe$drainage), data = hydro_summer_10304_kbe)+
  geom_point(alpha = 0.2)
```

## ARDLM
```{r}
ardlDlm_1 <- ardlDlm(formula = drainage ~ soilwater + precip + snowstorage + airtmp + glorad, data = hydro_summer_10304_kbe, p = 2, q = 2)

summary(ardlDlm_1)
```
## Fit_ardlm_across all members
```{r}
fit_ardlm_across <- function(formula, data, p, q, summary_output = TRUE) {
  results <- list()
  i <- 1
  for (member in unique(data$member)){
    data_loop <- data[data["member"] == member,]
    results[i] <- ardlDlm(formula = formula, data = data_loop, p = p, q = q)
    i <- i+ 1
  }
  
  lapply(results, summary)
}
```
##  Summer
```{r}
fit_ardlm_across(drainage ~ soilwater + precip + snowstorage + airtmp, data = hydro_summer_10304, p = 3, q = 4)
fit_ardlm_across(drainage ~ soilwater + precip + snowstorage + airtmp, data = hydro_summer_11502, p = 3, q = 4)
fit_ardlm_across(drainage ~ soilwater + precip + snowstorage + airtmp + glorad, data = hydro_summer_20203, p = 3, q = 4)


lapply(fit_ardlm_across(drainage ~ soilwater + precip + snowstorage + airtmp, data = hydro_summer_20203, p = 2, q = 2), "[[", "adj.r.squared")
```
## Winter
```{r}
fit_ardlm_across(drainage ~ soilwater + precip + snowstorage + airtmp, data = hydro_winter_10304, p = 3, q = 2)
fit_ardlm_across(drainage ~ soilwater + precip + snowstorage + airtmp, data = hydro_winter_11502, p = 3, q = 2)
fit_ardlm_across(drainage ~ soilwater + precip + snowstorage + airtmp, data = hydro_winter_20203, p = 3, q = 2)
```

## Distributed lag on Perception
```{r}
library(dlnm)
cb1.precip <- crossbasis(hydro_summer_10304_kbe$precip, lag= c(1,15), argvar=list(fun = "bs"),
                         arglag = list(fun="lin"))
gam_cb_precip <- gam(formula = lowlevel ~ soilwater + cb1.precip + snowstorage + airtmp, data = hydro_summer_10304_kbe, family = binomial())
summary(gam_cb_precip)

cb2.soilwater <- crossbasis(hydro_summer_10304_kbe$soilwater, lag= c(0,15), argvar=list(fun = "bs"),
                            arglag = list(fun="lin"))
gam_cb_soilwater <- gam(formula = lowlevel ~ cb2.soilwater + precip + snowstorage + airtmp, data = hydro_summer_10304_kbe, family = binomial())
summary(gam_cb_soilwater)

gam_cb_precip_soilwater <- gam(formula = lowlevel ~ cb2.soilwater + cb1.precip + snowstorage + s(airtmp) + precip * soilwater, data = hydro_summer_10304_kbe, family = binomial())
summary(gam_cb_precip_soilwater)

hydro_winter_10304_kbe <- hydro_winter_10304 %>% filter(member %in% "kbe")

gam_cb_precip_soilwater <- gam(formula = lowlevel ~ cb2.soilwater + cb1.precip + snowstorage + s(airtmp) + precip * soilwater, data = hydro_winter_10304_kbe, family = binomial())
summary(gam_cb_precip_soilwater)

model_list <- list()
i <- 1
for (member in hydro_winter_10304$member) {
  data <- hydro_winter_10304[hydro_winter_10304$member %in% member,]
  cb1.precip <- crossbasis(data$precip, lag= c(1,15), argvar=list(fun = "bs"), arglag = list(fun="lin"))
  cb2.soilwater <- crossbasis(data$soilwater, lag= c(0,15), argvar=list(fun = "bs"), arglag = list(fun="lin"))
  
  model_list[i] <- gam(formula = lowlevel ~ cb2.soilwater + cb1.precip + snowstorage + s(airtmp) + precip * soilwater, data = data, family = binomial())
  i= i +1
}

lapply(model_list, summary)
```

```{r}
rf_summer_11502 <- randomForest::randomForest(lowlevel ~ YY + avg_precip + avg_airtmp+ avg_glorad + avg_relhum + avg_soilwater + avg_snowstorage + groundwaterdepth + avg_infiltration + max_precip_scaled, data= hydro_summer_11502_kbe, proximity=TRUE, na.action = na.omit, importance = TRUE)

rf_winter_11502 <- randomForest::randomForest(lowlevel ~ YY + avg_precip + avg_airtmp+ avg_glorad + avg_relhum + avg_soilwater + avg_snowstorage + groundwaterdepth + avg_infiltration + max_precip_scaled, data= hydro_summer_11502_kbe, proximity=TRUE, na.action = na.omit, importance = TRUE)


rfPermute::plotImportance(rf_summer_11502)
rfPermute::plotImportance(rf_winter_11502)
```