###############################################################################
# Faktory ovlivňující poptávku po bikesharingu v Brně - Eva Běhálková, 510303
###############################################################################

# install.packages("car")
# install.packages("lmtest")
# install.packages("forecast")
# install.packages("sandwich")
# install.packages("olsrr")
# install.packages("readxl")
# install.packages("lubridate")
# install.packages("mgcv")
# install.packages("gridExtra")
# install.packages(c("broom", "writexl", "dplyr"))
# install.packages("strucchange")
# install.packages("gratia")
# install.packages("corrplot")

###############################################################################

library(readxl)
library(lubridate)
library(lmtest)
library(car)
library(sandwich)
library(forecast)
library(ggplot2)
library(gridExtra)
library(broom)
library(writexl)
library(dplyr)
library(mgcv)
library(strucchange)
library(gratia)
library(corrplot)

###############################################################################

my_data <- read_excel("C:/Users/evous/OneDrive - MUNI/Diplomová práce/Dataset.xlsx")

# Převedení 'date' z číselného Excel formátu na datum:
my_data$date <- as.Date(my_data$date, origin = "1899-12-30")

#Vytvoření faktorové proměnné měsíc (month):
my_data$month <- factor(format(my_data$date_time, "%m"),
                        levels = sprintf("%02d", 5:11),
                        labels = c("Květen", "Červen", "Červenec", "Srpen", "Září", "Říjen", "Listopad"))

# Vytvoření faktorové proměnné den v týdnu (weekday) + ref. hodnota neděle:
my_data$weekday <- factor(weekdays(my_data$date))
my_data$weekday <- relevel(my_data$weekday, ref = "neděle")

# Převedení hour na faktor:
my_data$hour <- factor(my_data$hour)

# Převedení temp_interval na faktor + ref. hodnota teplota 15 - 20 stupňů:
my_data$temp_interval <- factor(my_data$temp_interval,
                             levels = c("-5 – 0", "0 – 5", "5 – 10", "10 – 15", "15 – 20", "20 – 25", "25 – 30", "30 – 35"))
my_data$temp_interval <- relevel(my_data$temp_interval, ref = "15 – 20")

###############################################################################
# (1) OLS model log + log1p(rain)
###############################################################################

log_model_ols_rainlog <- lm(
  log1p(pocet_jizd) ~ wind + humidity + log1p(rain) + temperature + spicka + 
    hour + weekday + month,
  data = my_data
)

summary(log_model_ols_rainlog)

# Uložení do Excelu
ols_tidy <- tidy(log_model_ols_rainlog) %>%
  mutate(significance = case_when(
    p.value < 0.001 ~ "***",
    p.value < 0.01  ~ "**",
    p.value < 0.05  ~ "*",
    TRUE ~ ""
  ))
write_xlsx(ols_tidy, "C:/Users/evous/OneDrive - MUNI/Diplomová práce/Výstupy excel z Rka/OLS_model_1.xlsx")

###############################################################################
# Testování předpokladů pro OLS1
###############################################################################
# 1) Multikolinearita
vif(log_model_ols_rainlog)

# 2) Náhodné složky s nulovou střední hodnotou

#Průměr reziduí
mean(residuals(log_model_ols_rainlog))

#Histogram reziduí
qqnorm(residuals(log_model_ols_rainlog))
qqline(residuals(log_model_ols_rainlog), col = "red", lwd = 2)
hist(residuals(log_model_ols_rainlog), breaks = 50, main = "Histogram reziduí", xlab = "Rezidua")

# 3) Homoskedasticita
#Breusch_Pagan test
bptest(log_model_ols_rainlog)

#White HC3
white_hc3 <- coeftest(log_model_ols_rainlog, vcov = vcovHC(log_model_ols_rainlog, type = "HC3"))

# Uložení do Excelu
white_df <- data.frame(
  Proměnná = rownames(white_hc3),
  Odhad = white_hc3[, 1],
  Std_chyba = white_hc3[, 2],
  t_hodnota = white_hc3[, 3],
  p_hodnota = white_hc3[, 4]
)
white_df$Hvezdicky <- dplyr::case_when(
  white_df$p_hodnota < 0.001 ~ "***",
  white_df$p_hodnota < 0.01  ~ "**",
  white_df$p_hodnota < 0.05  ~ "*",
  white_df$p_hodnota < 0.1   ~ ".",
  TRUE ~ ""
)
write_xlsx(white_df, "C:/Users/evous/OneDrive - MUNI/Diplomová práce/Výstupy excel z Rka/white_hc3.xlsx")

#Newey-West standard errors
newey_west <- coeftest(log_model_ols_rainlog, vcov = NeweyWest(log_model_ols_rainlog, lag = 1, prewhite = FALSE))

# Uložení do Excelu
newey_df <- data.frame(
  Proměnná = rownames(newey_west),
  Odhad = newey_west[, 1],
  Std_chyba = newey_west[, 2],
  t_hodnota = newey_west[, 3],
  p_hodnota = newey_west[, 4]
)
newey_df$Hvezdicky <- dplyr::case_when(
  newey_df$p_hodnota < 0.001 ~ "***",
  newey_df$p_hodnota < 0.01  ~ "**",
  newey_df$p_hodnota < 0.05  ~ "*",
  newey_df$p_hodnota < 0.1   ~ ".",
  TRUE ~ ""
)
write_xlsx(newey_df, "C:/Users/evous/OneDrive - MUNI/Diplomová práce/Výstupy excel z Rka/newey_west.xlsx")

#Graf rezidua vs predikované hodnoty
df_ols <- data.frame(fitted = fitted(log_model_ols_rainlog), residuals = residuals(log_model_ols_rainlog))

ggplot(df_ols, aes(x = fitted, y = residuals)) +
  geom_point(color = "darkgreen", alpha = 0.5) +
  geom_smooth(method = "loess", color = "red") +
  geom_hline(yintercept = 0, linetype = "dashed", color = "black") +
  labs(title = "Rezidua vs Predikované hodnoty (OLS model)",
       x = "Predikované hodnoty",
       y = "Rezidua") +
  theme_minimal()

# 4) Nekorelovanost
#Durbin-Watson test
dwtest(log_model_ols_rainlog)

# 5) Konstantní parametry
#F-testy porovnávajících modely s a bez interakce sezónního členění
# Vytvoření nové proměnné sezona na základě month
my_data$sezona <- dplyr::case_when(
  my_data$month %in% c("Květen", "Červen") ~ "jaro",
  my_data$month %in% c("Červenec", "Srpen") ~ "leto",
  my_data$month %in% c("Září", "Říjen", "Listopad") ~ "podzim",
  TRUE ~ NA_character_
)
my_data$sezona <- factor(my_data$sezona, levels = c("jaro", "leto", "podzim"))

# Podmnožiny pro dvojice sezón
data_jaro_leto <- subset(my_data, sezona %in% c("jaro", "leto"))
data_jaro_podzim <- subset(my_data, sezona %in% c("jaro", "podzim"))
data_leto_podzim <- subset(my_data, sezona %in% c("leto", "podzim"))

# Modely bez interakce
model_spojeny_jaro_leto <- lm(pocet_jizd ~ wind + humidity + rain + temperature + spicka, data = data_jaro_leto)
model_spojeny_jaro_podzim <- lm(pocet_jizd ~ wind + humidity + rain + temperature + spicka, data = data_jaro_podzim)
model_spojeny_leto_podzim <- lm(pocet_jizd ~ wind + humidity + rain + temperature + spicka, data = data_leto_podzim)

# Modely s interakcí se sezónou
model_interakce_jaro_leto <- lm(pocet_jizd ~ sezona * (wind + humidity + rain + temperature + spicka), data = data_jaro_leto)
model_interakce_jaro_podzim <- lm(pocet_jizd ~ sezona * (wind + humidity + rain + temperature + spicka), data = data_jaro_podzim)
model_interakce_leto_podzim <- lm(pocet_jizd ~ sezona * (wind + humidity + rain + temperature + spicka), data = data_leto_podzim)

# F-testy porovnávající modely bez a s interakcí
cat("=== Jaro vs Léto ===\n")
print(anova(model_spojeny_jaro_leto, model_interakce_jaro_leto))

cat("\n=== Jaro vs Podzim ===\n")
print(anova(model_spojeny_jaro_podzim, model_interakce_jaro_podzim))

cat("\n=== Léto vs Podzim ===\n")
print(anova(model_spojeny_leto_podzim, model_interakce_leto_podzim))

# 6) Lineární vztah mezi vysvětlující a vysvětlovanou proměnnou 
#Ramsey RESET test
resettest(log_model_ols_rainlog)

#Grafy parciálních regresí
crPlots(log_model_ols_rainlog)

# 7) Normalita reziduí
#Kolmogorov-Smirnovův test
ks.test(residuals(log_model_ols_rainlog), "pnorm", mean = mean(residuals(log_model_ols_rainlog)), sd = sd(residuals(log_model_ols_rainlog)))

#QQ plot
qqnorm(residuals(log_model_ols_rainlog))
qqline(residuals(log_model_ols_rainlog), col = "red")

#Dodatečná diagnostika pro detekci vlivných bodů

#Cookova vzdálenost
plot(log_model_ols_rainlog, which = 4)
influence.measures(log_model_ols_rainlog)


###############################################################################
# (2) OLS model s log + rain_binary
###############################################################################
log_model_ols_rainbinary <- lm(
  log1p(pocet_jizd) ~ wind + humidity + rain_binary + temperature + spicka + 
    hour + weekday + month,
  data = my_data
)

summary(log_model_ols_rainbinary)

ols_tidy2 <- tidy(log_model_ols_rainbinary) %>%
  mutate(significance = case_when(
    p.value < 0.001 ~ "***",
    p.value < 0.01  ~ "**",
    p.value < 0.05  ~ "*",
    TRUE ~ ""
  ))
write_xlsx(ols_tidy2, "C:/Users/evous/OneDrive - MUNI/Diplomová práce/Výstupy excel z Rka/OLS_model_2.xlsx")

###############################################################################
# Testování předpokladů pro OLS2
###############################################################################
# 1) Multikolinearita
vif(log_model_ols_rainbinary)

# 2) Náhodné složky s nulovou střední hodnotou

#Průměr reziduí
mean(residuals(log_model_ols_rainbinary))

#Histogram reziduí
qqnorm(residuals(log_model_ols_rainbinary))
qqline(residuals(log_model_ols_rainbinary), col = "red", lwd = 2)
hist(residuals(log_model_ols_rainbinary), breaks = 50, main = "Histogram reziduí", xlab = "Rezidua")

# 3) Homoskedasticita
#Breusch_Pagan test
bptest(log_model_ols_rainbinary)

#White HC3
white_hc3 <- coeftest(log_model_ols_rainbinary, vcov = vcovHC(log_model_ols_rainbinary, type = "HC3"))

# Uložení do Excelu
white_df <- data.frame(
  Proměnná = rownames(white_hc3),
  Odhad = white_hc3[, 1],
  Std_chyba = white_hc3[, 2],
  t_hodnota = white_hc3[, 3],
  p_hodnota = white_hc3[, 4]
)
white_df$Hvezdicky <- dplyr::case_when(
  white_df$p_hodnota < 0.001 ~ "***",
  white_df$p_hodnota < 0.01  ~ "**",
  white_df$p_hodnota < 0.05  ~ "*",
  white_df$p_hodnota < 0.1   ~ ".",
  TRUE ~ ""
)
write_xlsx(white_df, "C:/Users/evous/OneDrive - MUNI/Diplomová práce/Výstupy excel z Rka/white_hc3_2.xlsx")

#Newey-West standard errors
newey_west <- coeftest(log_model_ols_rainbinary, vcov = NeweyWest(log_model_ols_rainbinary, lag = 1, prewhite = FALSE))

# Uložení do Excelu
newey_df <- data.frame(
  Proměnná = rownames(newey_west),
  Odhad = newey_west[, 1],
  Std_chyba = newey_west[, 2],
  t_hodnota = newey_west[, 3],
  p_hodnota = newey_west[, 4]
)
newey_df$Hvezdicky <- dplyr::case_when(
  newey_df$p_hodnota < 0.001 ~ "***",
  newey_df$p_hodnota < 0.01  ~ "**",
  newey_df$p_hodnota < 0.05  ~ "*",
  newey_df$p_hodnota < 0.1   ~ ".",
  TRUE ~ ""
)
write_xlsx(newey_df, "C:/Users/evous/OneDrive - MUNI/Diplomová práce/Výstupy excel z Rka/newey_west_2.xlsx")

#Graf rezidua vs predikované hodnoty
df_ols <- data.frame(fitted = fitted(log_model_ols_rainbinary), residuals = residuals(log_model_ols_rainbinary))

ggplot(df_ols, aes(x = fitted, y = residuals)) +
  geom_point(color = "darkgreen", alpha = 0.5) +
  geom_smooth(method = "loess", color = "red") +
  geom_hline(yintercept = 0, linetype = "dashed", color = "black") +
  labs(title = "Rezidua vs Predikované hodnoty (OLS model)",
       x = "Predikované hodnoty",
       y = "Rezidua") +
  theme_minimal()

# 4) Nekorelovanost
#Durbin-Watson test
dwtest(log_model_ols_rainbinary)

# 5) Konstantní parametry
#F-testy porovnávajících modely s a bez interakce sezónního členění
# Vytvoření nové proměnné sezona na základě month
my_data$sezona <- dplyr::case_when(
  my_data$month %in% c("Květen", "Červen") ~ "jaro",
  my_data$month %in% c("Červenec", "Srpen") ~ "leto",
  my_data$month %in% c("Září", "Říjen", "Listopad") ~ "podzim",
  TRUE ~ NA_character_
)
my_data$sezona <- factor(my_data$sezona, levels = c("jaro", "leto", "podzim"))

# Podmnožiny pro dvojice sezón
data_jaro_leto <- subset(my_data, sezona %in% c("jaro", "leto"))
data_jaro_podzim <- subset(my_data, sezona %in% c("jaro", "podzim"))
data_leto_podzim <- subset(my_data, sezona %in% c("leto", "podzim"))

# Modely bez interakce
model_spojeny_jaro_leto <- lm(pocet_jizd ~ wind + humidity + rain_binary + temperature + spicka, data = data_jaro_leto)
model_spojeny_jaro_podzim <- lm(pocet_jizd ~ wind + humidity + rain_binary + temperature + spicka, data = data_jaro_podzim)
model_spojeny_leto_podzim <- lm(pocet_jizd ~ wind + humidity + rain_binary + temperature + spicka, data = data_leto_podzim)

# Modely s interakcí se sezónou
model_interakce_jaro_leto <- lm(pocet_jizd ~ sezona * (wind + humidity + rain_binary + temperature + spicka), data = data_jaro_leto)
model_interakce_jaro_podzim <- lm(pocet_jizd ~ sezona * (wind + humidity + rain_binary + temperature + spicka), data = data_jaro_podzim)
model_interakce_leto_podzim <- lm(pocet_jizd ~ sezona * (wind + humidity + rain_binary + temperature + spicka), data = data_leto_podzim)

# F-testy porovnávající modely bez a s interakcí
cat("=== Jaro vs Léto ===\n")
print(anova(model_spojeny_jaro_leto, model_interakce_jaro_leto))

cat("\n=== Jaro vs Podzim ===\n")
print(anova(model_spojeny_jaro_podzim, model_interakce_jaro_podzim))

cat("\n=== Léto vs Podzim ===\n")
print(anova(model_spojeny_leto_podzim, model_interakce_leto_podzim))


# 6) Lineární vztah mezi vysvětlující a vysvětlovanou proměnnou 
#Ramsey RESET test
resettest(log_model_ols_rainbinary)

#Grafy parciálních regresí
crPlots(log_model_ols_rainbinary)

# 7) Normalita reziduí
#Kolmogorov-Smirnovův test
ks.test(residuals(log_model_ols_rainbinary), "pnorm", mean = mean(residuals(log_model_ols_rainbinary)), sd = sd(residuals(log_model_ols_rainbinary)))

#QQ plot
qqnorm(residuals(log_model_ols_rainbinary))
qqline(residuals(log_model_ols_rainbinary), col = "red")

#Dodatečná diagnostika pro detekci vlivných bodů

#Cookova vzdálenost
plot(log_model_ols_rainbinary, which = 4)
influence.measures(log_model_ols_rainbinary)


###############################################################################
# (3) OLS model s temp_interval
###############################################################################
log_model_ols_temp_int <- lm(
  log1p(pocet_jizd) ~ wind + humidity + log1p(rain) + temp_interval + spicka + 
    hour + weekday + month,
  data = my_data
)

summary(log_model_ols_temp_int)

# Uložení do Excelu
ols_tidy <- tidy(log_model_ols_temp_int) %>%
  mutate(significance = case_when(
    p.value < 0.001 ~ "***",
    p.value < 0.01  ~ "**",
    p.value < 0.05  ~ "*",
    p.value < 0.1   ~ ".",
    TRUE ~ ""
  ))
write_xlsx(ols_tidy, "C:/Users/evous/OneDrive - MUNI/Diplomová práce/Výstupy excel z Rka/OLS_model_3.xlsx")

###############################################################################
# Testování předpokladů OLS3
###############################################################################
# 1) Multikolinearita
vif(log_model_ols_temp_int)

# 2) Náhodné složky s nulovou střední hodnotou
#Průměr reziduí
mean(residuals(log_model_ols_temp_int))

#Histogram reziduí
qqnorm(residuals(log_model_ols_temp_int))
qqline(residuals(log_model_ols_temp_int), col = "red", lwd = 2)
hist(residuals(log_model_ols_temp_int), breaks = 50, main = "Histogram reziduí", xlab = "Rezidua")

# 3) Homoskedasticita
#Breusch_Pagan test
bptest(log_model_ols_temp_int)

#White HC3
white_hc3 <- coeftest(log_model_ols_temp_int, vcov = vcovHC(log_model_ols_temp_int, type = "HC3"))

# Uložení do Excelu
white_df <- data.frame(
  Proměnná = rownames(white_hc3),
  Odhad = white_hc3[, 1],
  Std_chyba = white_hc3[, 2],
  t_hodnota = white_hc3[, 3],
  p_hodnota = white_hc3[, 4]
)
white_df$Hvezdicky <- dplyr::case_when(
  white_df$p_hodnota < 0.001 ~ "***",
  white_df$p_hodnota < 0.01  ~ "**",
  white_df$p_hodnota < 0.05  ~ "*",
  white_df$p_hodnota < 0.1   ~ ".",
  TRUE ~ ""
)
write_xlsx(white_df, "C:/Users/evous/OneDrive - MUNI/Diplomová práce/Výstupy excel z Rka/white_hc3_3.xlsx")

#Newey-West standard errors
newey_west <- coeftest(log_model_ols_temp_int, vcov = NeweyWest(log_model_ols_temp_int, lag = 1, prewhite = FALSE))

# Uložení do Excelu
newey_df <- data.frame(
  Proměnná = rownames(newey_west),
  Odhad = newey_west[, 1],
  Std_chyba = newey_west[, 2],
  t_hodnota = newey_west[, 3],
  p_hodnota = newey_west[, 4]
)
newey_df$Hvezdicky <- dplyr::case_when(
  newey_df$p_hodnota < 0.001 ~ "***",
  newey_df$p_hodnota < 0.01  ~ "**",
  newey_df$p_hodnota < 0.05  ~ "*",
  newey_df$p_hodnota < 0.1   ~ ".",
  TRUE ~ ""
)
write_xlsx(newey_df, "C:/Users/evous/OneDrive - MUNI/Diplomová práce/Výstupy excel z Rka/newey_west_3.xlsx")

#Graf rezidua vs predikované hodnoty
df_ols <- data.frame(fitted = fitted(log_model_ols_temp_int), residuals = residuals(log_model_ols_temp_int))

ggplot(df_ols, aes(x = fitted, y = residuals)) +
  geom_point(color = "darkgreen", alpha = 0.5) +
  geom_smooth(method = "loess", color = "red") +
  geom_hline(yintercept = 0, linetype = "dashed", color = "black") +
  labs(title = "Rezidua vs Predikované hodnoty (OLS model)",
       x = "Predikované hodnoty",
       y = "Rezidua") +
  theme_minimal()

# 4) Nekorelovanost
#Durbin-Watson test
dwtest(log_model_ols_temp_int)

# 5) Konstantní parametry
#F-testy porovnávajících modely s a bez interakce sezónního členění
# Vytvoření nové proměnné sezona na základě month
my_data$sezona <- dplyr::case_when(
  my_data$month %in% c("Květen", "Červen") ~ "jaro",
  my_data$month %in% c("Červenec", "Srpen") ~ "leto",
  my_data$month %in% c("Září", "Říjen", "Listopad") ~ "podzim",
  TRUE ~ NA_character_
)
my_data$sezona <- factor(my_data$sezona, levels = c("jaro", "leto", "podzim"))

# Podmnožiny pro dvojice sezón
data_jaro_leto <- subset(my_data, sezona %in% c("jaro", "leto"))
data_jaro_podzim <- subset(my_data, sezona %in% c("jaro", "podzim"))
data_leto_podzim <- subset(my_data, sezona %in% c("leto", "podzim"))

# Modely bez interakce
model_spojeny_jaro_leto <- lm(pocet_jizd ~ wind + humidity + rain + temp_interval + spicka, data = data_jaro_leto)
model_spojeny_jaro_podzim <- lm(pocet_jizd ~ wind + humidity + rain + temp_interval + spicka, data = data_jaro_podzim)
model_spojeny_leto_podzim <- lm(pocet_jizd ~ wind + humidity + rain + temp_interval + spicka, data = data_leto_podzim)

# Modely s interakcí se sezónou
model_interakce_jaro_leto <- lm(pocet_jizd ~ sezona * (wind + humidity + rain + temp_interval + spicka), data = data_jaro_leto)
model_interakce_jaro_podzim <- lm(pocet_jizd ~ sezona * (wind + humidity + rain + temp_interval + spicka), data = data_jaro_podzim)
model_interakce_leto_podzim <- lm(pocet_jizd ~ sezona * (wind + humidity + rain + temp_interval + spicka), data = data_leto_podzim)

# F-testy porovnávající modely bez a s interakcí
cat("=== Jaro vs Léto ===\n")
print(anova(model_spojeny_jaro_leto, model_interakce_jaro_leto))

cat("\n=== Jaro vs Podzim ===\n")
print(anova(model_spojeny_jaro_podzim, model_interakce_jaro_podzim))

cat("\n=== Léto vs Podzim ===\n")
print(anova(model_spojeny_leto_podzim, model_interakce_leto_podzim))

# 6) Lineární vztah mezi vysvětlující a vysvětlovanou proměnnou 
#Ramsey RESET test
resettest(log_model_ols_temp_int)

#Grafy parciálních regresí
crPlots(log_model_ols_temp_int)

# 7) Normalita reziduí
#Kolmogorov-Smirnovův test
ks.test(residuals(log_model_ols_temp_int), "pnorm", mean = mean(residuals(log_model_ols_temp_int)), sd = sd(residuals(log_model_ols_temp_int)))

#QQ plot
qqnorm(residuals(log_model_ols_temp_int))
qqline(residuals(log_model_ols_temp_int), col = "red")

#Dodatečná diagnostika pro detekci vlivných bodů

#Cookova vzdálenost
plot(log_model_ols_temp_int, which = 4)
influence.measures(log_model_ols_temp_int)


###############################################################################
# MODEL GAM
###############################################################################

# Testování rozdílných modelů (FINÁLNÍ MODEL = GAM2)
gam_model1 <- gam(
  log1p(pocet_jizd) ~ 
    s(temperature) + 
    s(humidity) + 
    s(wind) + 
    s(rain) + 
    spicka +
    factor(hour) + 
    factor(weekday) +
    factor(month),
  data = my_data
)

gam_model2 <- gam(
  log1p(pocet_jizd) ~ 
    s(temperature) + 
    s(humidity) + 
    s(wind) + 
    s(log1p(rain)) + 
    spicka +
    factor(hour) + 
    factor(weekday)+
    factor(month),
  data = my_data
)

gam_model3 <- gam(
  log1p(pocet_jizd) ~ 
    s(temperature) + 
    s(humidity) + 
    s(wind) + 
    rain_binary + 
    spicka +
    factor(hour) + 
    factor(weekday) +
    factor(month),
  data = my_data
)

summary(gam_model1)
summary(gam_model2)
summary(gam_model3)
AIC(gam_model1, gam_model2, gam_model3)

s <- summary(gam_model2)
param_table <- as.data.frame(s$p.table)
param_table$Term <- rownames(s$p.table)
get_signif <- function(p) {
  if (is.na(p)) return(" ")
  else if (p < 0.001) return("***")
  else if (p < 0.01) return("**")
  else if (p < 0.05) return("*")
  else if (p < 0.1) return(".")
  else return(" ")
}
param_table$Signif <- sapply(param_table$`Pr(>|t|)`, get_signif)
param_table <- param_table[, c("Term", "Estimate", "Std. Error", "t value", "Pr(>|t|)", "Signif")]
rownames(param_table) <- NULL
smooth_terms <- as.data.frame(s$s.table)
smooth_terms$Term <- rownames(s$s.table)
smooth_terms$Signif <- sapply(smooth_terms[, "p-value"], get_signif)
rownames(smooth_terms) <- NULL

# Uložit do Excelu
write_xlsx(list(
    "Parametric_coefficients" = param_table,
    "Smooth_terms" = smooth_terms), path = "C:/Users/evous/OneDrive - MUNI/Diplomová práce/Výstupy excel z Rka/gam_model2_summary.xlsx")

# Grafy hladkých funkcí pro spline proměnné
plot(gam_model2, pages = 1, shade = TRUE, se = TRUE)

#Diagnostika
gam.check(gam_model2)

###############################################################################
# Testování předpokladů GAM
###############################################################################
# Autokorelace
dwtest(gam_model3$residuals ~ 1)
acf(residuals(gam_model3), main = "ACF reziduí GAM modelu")

#Multikolinearita

#Heatmapa korelační matice surových prediktorů
prediktory <- my_data[, c("temperature", "humidity", "wind", "rain")]
cor_matrix <- cor(prediktory, use = "complete.obs")
print(cor_matrix)

corrplot(cor_matrix, method = "color", addCoef.col = "black",
         tl.col = "black", tl.srt = 45, number.cex = 0.8)

#VIF
lm_temp <- lm(log1p(pocet_jizd) ~ temperature + humidity + wind + rain, data = my_data)
vif(lm_temp)

###############################################################################
# Rozšířená analýza
###############################################################################
# Interakce temperature a month
gam_model_inter_temp_month <- gam(
  log1p(pocet_jizd) ~ 
    ti(temperature, month, bs = c("tp", "re")) +
    s(humidity) +
    s(wind) +
    s(log1p(rain)) +
    spicka +
    factor(hour) +
    factor(weekday),
  data = my_data
)

plot(gam_model_inter_temp_month, pages = 1, shade = TRUE)
summary(gam_model_inter_temp_month)

my_data$month_num <- as.numeric(my_data$month)
gam_model_inter_temp_month2 <- gam(
  log1p(pocet_jizd) ~ 
    ti(temperature, month_num, bs = c("tp", "tp")) +
    s(humidity) +
    s(wind) +
    s(log1p(rain)) +
    spicka +
    factor(hour) +
    factor(weekday),
  data = my_data
)

draw(gam_model_inter_temp_month2, select = 1) +
  scale_y_continuous(
    breaks = 1:7,
    labels = c("Květen", "Červen", "Červenec", "Srpen", "Září", "Říjen", "Listopad")
  ) +
  labs(
    title = "Interakce mezi teplotou a měsícem",
    x = "Teplota (°C)",
    y = "Měsíc"
  )

# Interakce rain a den v týdnu
my_data$weekday_num <- as.numeric(factor(my_data$weekday, 
                                         levels = c("pondělí", "úterý", "středa", "čtvrtek", "pátek", "sobota", "neděle")))

gam_model_inter_rain_weekday2 <- gam(
  log1p(pocet_jizd) ~ 
    ti(log1p(rain), weekday_num, bs = c("tp", "re")) +
    s(temperature) +
    s(humidity) +
    s(wind) +
    spicka +
    factor(hour) +
    factor(month),
  data = my_data
)

draw(gam_model_inter_rain_weekday2, select = 1) +
  scale_y_continuous(
    breaks = 1:7,
    labels = c("Pondělí", "Úterý", "Středa", "Čtvrtek", "Pátek", "Sobota", "Neděle")
  ) +
  labs(
    title = "Interakce srážek a dne v týdnu",
    x = "log(1 + srážky v mm)",
    y = "Den v týdnu"
  )

# Derivative analysis
## temperature
deriv_temp <- derivatives(gam_model2,
                          term = "s(temperature)", 
                          interval = "confidence")
draw(deriv_temp) +
  labs(
    title = "Derivace vlivu teploty na počet jízd",
    x = "Teplota (°C)",
    y = "Změna logaritmu počtu jízd při změně teploty"
  )


#Parciální efekty - temperature, humidity, rain, wind
draw(gam_model2, select = "s(temperature)") +
  labs(
    title = "Parciální efekt teploty na počet jízd",
    x = "Teplota (°C)",
    y = "Parciální efekt na logaritmus počtu jízd"
  )

draw(gam_model2, select = "s(humidity)") +
  labs(
    title = "Parciální efekt vlhkosti na počet jízd",
    x = "Vlhkost vzduchu (%)",
    y = "Parciální efekt na logaritmus počtu jízd"
  )


draw(gam_model2, select = "s(log1p(rain))") +
  labs(
    title = "Parciální efekt srážek na počet jízd",
    x = "log(1 + srážky) [mm]",
    y = "Parciální efekt na logaritmus počtu jízd"
  )

draw(gam_model2, select = "s(wind)") +
  labs(
    title = "Parciální efekt rychlosti větru na počet jízd",
    x = "Rychlost větru (m/s)",
    y = "Parciální efekt na logaritmus počtu jízd"
  )