# ANALYSE ECOSYSTEM SERVICES

################################################################################
# load libraries

library(ggplot2)
library(nlme)
library(gamm4)
library(mgcv)
library(dplyr)
# install.packages('gamm4')
# library(dplyr)

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

# 1. LOAD DATA ----

ES_df <- read.csv("ECOSER_yearly.csv")
SCA_ES_df <- read.csv("SCALED_ECOSER_yearly.csv")
SCA_ES_df_rolling <- read.csv("SCALED_ECOSER_yearly_rolling.csv")
SCA_ES_df2 <- read.csv("SCALED_ECOSER_mean.csv")

WOOD_AGB <- as.numeric(ES_df$WoodAGB)
date <- as.numeric(ES_df$date)
PLOT <- as.factor(ES_df$PLOT)
PROJECTION <- as.factor(ES_df$PROJECTION)
MANAG_INTENSITY <- as.factor(ES_df$MANAG_INTENSITY)
dat <- data.frame(WOOD_AGB, date, PLOT, PROJECTION, MANAG_INTENSITY)
dat <- subset(dat, WOOD_AGB != 0)

# WOOD_AGB <- as.numeric(SCA_ES_df$SCA_WoodAGB)
# date <- as.numeric(SCA_ES_df$date)
# PLOT <- as.factor(SCA_ES_df$PLOT)
# PROJECTION <- as.factor(SCA_ES_df$PROJECTION)
# MANAG_INTENSITY <- as.factor(SCA_ES_df$MANAG_INTENSITY)
# dat <- data.frame(WOOD_AGB, date, PLOT, PROJECTION, MANAG_INTENSITY)

# 2. MODELS ----
# Test the effect of Management Intensity and Climate over ES

# TEST GAMM

gamm1 <- mgcv::gamm(WOOD_AGB ~
                      s(date, by = MANAG_INTENSITY, bs = "sz"), # +
                      # s(date, by = PROJECTION, bs = "fs"),
              random = list(PLOT = ~ 1),
              data = dat,
              family = Gamma(link = "log"))
test <- mgcv::gam(WOOD_AGB ~
                    s(date, by = MANAG_INTENSITY, bs = "sz"), # +
                  # s(date, by = PROJECTION, bs = "fs"),
                    # random = list(PLOT = ~ 1),
                    data = dat,
                    family = Gamma(link = "log"))
AIC(gamm1, test)
summary(gamm1$gam, cor = F)
anova(gamm1$gam)
plot(gamm1$gam)

summary(gamm1$lme, cor = F)
plot(gamm1$lme)

fit <- sqrt(fitted(gamm1$lme))
obs <- sqrt(dat$WOOD_AGB)
plot(fit, obs);abline() # must be a straigth line

E4 <- resid(gamm1$lme, type = "pearson")
plot(E4, obs);abline() # must be a band with no patterns

residuals <- resid(gamm1$lme, type = "normalized")
plot(residuals, fit);abline() # must be a clear cone

# continue with gamms # https://bart-larsen.github.io/GAMM-Tutorial/
'install.packages('gratia')'

gratia::smooth_estimates(gamm1, "date")
gratia::appraise(gamm1$gam)

# Calculate the derivative
d <- gratia::derivatives(gamm1, n = 1000)

# Take a look at it using gratia::draw
d_plot <- gratia::draw(d)
print(d_plot) # this shows at with values of the smooth term there is significance

# Identify significant areas # Not working
d <- d |>
  mutate(sig = !(0 > .lower_ci & 0 < .upper_ci)) #Ages where the CI does not include zero
# cat(sprintf(
#   "\nSignificant change: %1.2f - %1.2f\n",
#   min(d$data[d$sig == T]),
#   max(d$data[d$sig == T])
# )) # this only works properly if there is one contiguous area of significant change (sorry)

# Add some ornaments to the plot
d_plot <- ggplot(data = d,
                 aes(
                   x = date,
                   y = .derivative,
                   color = sig,
                   ymin = .lower_ci,
                   ymax = .upper_ci
                 )) +
  facet_grid(. ~ MANAG_INTENSITY) +
  geom_ribbon(fill = "black",
                                  alpha = .3,
                                  color = NA) +
  geom_line(size = 1, show.legend = F) +
  scale_color_manual(values = c("TRUE" = "firebrick", "FALSE" = "black")) +
  geom_hline(yintercept = 0, linetype = 2)

# Let's visualize using our visualizer function

big_plot <- gratia::visualize_model(
  modobj = gamm1,
  smooth_var = "date",
  int_var = "MANAG_INTENSITY",
  group_var = "bblid",
  plabels = "scatter + derivative",
  derivative_plot = T
)





## Test -------

lme1 <- lme(WoodAGB ~ date + PROJECTION * MANAG_INTENSITY,
            random = ~ 1 | PLOT,
            # correlation = corARMA(p = 1, form = ~ date | PLOT, q = 0),
            data = ES_df)
plot(lme1)
mod_final_res <- residuals(lme1)
qqnorm(mod_final_res)
qqline(mod_final_res, lty=2, col="red")
summary(lme1)

# ----

lme1 <- lme(WOOD_AGB ~ date + PROJECTION + MANAG_INTENSITY,
    random = ~1 | PLOT,
    # correlation = corARMA(p=1,form = ~date | PLOT, q = 0),
    data = dat)
plot(lme1)
mod_final_res<-residuals(lme1)
qqnorm(mod_final_res)
qqline(mod_final_res, lty=2, col="red")

# -----

beta1 <- glmmTMB::glmmTMB(WOOD_AGB ~ date + PROJECTION + MANAG_INTENSITY + (1 | PLOT),
                 family = glmmTMB::beta_family(link="logit"),
                 data = dat)
summary(beta1)

res1 <- glmmTMB::simulateResiduals(beta1)
plot(res1)

mod_final_res<-residuals(beta1)
plot(mod_final_res)
qqnorm(mod_final_res)
qqline(mod_final_res, lty=2, col="red")

#










# lme(SCA_WoodAGB ~ date + PROJECTION*MANAG_INTENSITY,
#     random = ~ 1 | PLOT,
#     correlation = corARMA(form = 1 | PLOT, p = 1, q = 0),
#     data = dat)

lme(BAI~date,
    random=~1|SP/PLOT/ID,
    correlation=corARMA(p=3,form=~YEAR|SP/PLOT/ID,q=0),data=pinus_data)