First commit with the code Hector shared with me

10_EXPORTS/03_analyse-ecoser.R

@@ -0,0 +1,186 @@
+# 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)
+
+
+