rm(list = ls())

require(glmmTMB)
require(DHARMa)
require(performance)
require(ggplot2)
require(dplyr)
require(tidyr)
require(car)
require(ggeffects)
require(scales)
require(cowplot)

df<-read.csv("C:/path/data_plot.csv")

df$lgtsr = log2(df$tsr)
df$lgfwd = log10(df$fwd_ha)
df$lgcwd = ifelse(df$cwd_vol_tot_ha>0, log10(df$cwd_vol_tot_ha), NA)
df$cwd_yesno = ifelse(df$cwd_vol_tot>0, 1, 0)
df$lgt4decay = log(df$t4_decay)
df$lgcwddecay = log(df$cwd_decay)

eps <- min(df$MPD[df$MPD > 0], na.rm=TRUE) / 100  # offset for log transformation
eps

df = df %>%
  mutate(sc.sl = scale(slope),
         sc.tsr = scale(log2(tsr)),
         sc.cc = scale(canopy),
         sc.north = scale(north),
         sc.east = scale(east),
         sc.slope = scale(slope),         
         sc.prod = scale(productivity),
         sc.mpd_l = scale(log10(MPD+3.25)),
         sc.wdens = scale(w_density),
         sc.cwd=scale(lgcwd),
         )

## FWD ####

m1 = glmmTMB(lgfwd ~ sc.cc + sc.tsr + sc.north + sc.east +
               sc.prod + sc.mpd_l + sc.wdens + sc.slope +site, df)

summary(m1)
testResiduals(m1)
plot(simulateResiduals(m1))
Anova(m1)
check_collinearity(m1)

## CWD ####

m3 = glmmTMB(lgcwd ~ sc.tsr + sc.prod + sc.mpd_l + sc.cc +
             sc.wdens + sc.north + sc.east + sc.slope + (1|site), df[!is.na(df$lgcwd),])
testResiduals(m3)
plot(simulateResiduals(m3))
summary(m3)
Anova(m3)
check_collinearity(m3)

## DECAY ####

df_clean <- df[is.finite(df$lgt4decay), ]
m4 = glmmTMB(lgt4decay ~ sc.tsr + sc.cc + sc.north + lgcwd + sc.east + sc.slope + site, data=df_clean)

testResiduals(m4)
plot(simulateResiduals(m4))
summary(m4)
Anova(m4)
check_collinearity(m4)


## CV % ####

# mean and sd values for ms text
df %>%group_by(site, tsr)%>%
  summarise(
    #sum_cwd=sum(cwd_vol_tot_ha, na.rm=TRUE),
    mean_cwd = mean(cwd_vol_tot_ha, na.rm = TRUE),
    sd_cwd = sd(cwd_vol_tot_ha, na.rm = TRUE),
    #sum_t4=sum(t4_vol_tot_ha, na.rm=TRUE),
    #mean_t4 = mean(t4_vol_tot_ha, na.rm = TRUE),
    #sd_t4 = sd(t4_vol_tot_ha, na.rm = TRUE),
    #sum_fwd=sum(fwd_ha, na.rm=TRUE),
    mean_fwd = mean(fwd_ha, na.rm = TRUE),
    sd_fwd = sd(fwd_ha, na.rm = TRUE),
    .groups = "drop")



# cv%
library(dplyr)
library(tidyr)
library(tibble)

cv_boot_std <- function(x, k, nboot = 1000, replace = FALSE) {
  x <- na.omit(x)
  n <- length(x)
  k <- as.integer(floor(k))
  
  empty_out <- tibble(
    cv_obs = NA_real_,
    boot_median = NA_real_,
    ci_lower = NA_real_,
    ci_upper = NA_real_,
    boot_se = NA_real_,
    n = n,
    k = k
  )
  
  if (n < 2L || is.na(k) || k < 2L) return(empty_out)
  if (!replace && n < k) stop("k is larger than the available sample size and replace = FALSE.")
  
  x_mean <- mean(x)
  if (is.na(x_mean) || x_mean == 0) return(empty_out)
  
  cv_obs <- sd(x) / x_mean * 100
  
  boot_vals <- replicate(nboot, {
    xb <- sample(x, size = k, replace = replace)
    xb_mean <- mean(xb)
    
    if (length(xb) < 2L || is.na(xb_mean) || xb_mean == 0) {
      NA_real_
    } else {
      sd(xb) / xb_mean * 100
    }
  })
  
  boot_vals_ok <- boot_vals[!is.na(boot_vals)]
  
  if (length(boot_vals_ok) == 0L) {
    return(tibble(
      cv_obs = cv_obs,
      boot_median = NA_real_,
      ci_lower = NA_real_,
      ci_upper = NA_real_,
      boot_se = NA_real_,
      n = n,
      k = k
    ))
  }
  
  tibble(
    cv_obs = cv_obs,
    boot_median = median(boot_vals_ok),
    ci_lower = unname(quantile(boot_vals_ok, 0.025)),
    ci_upper = unname(quantile(boot_vals_ok, 0.975)),
    boot_se = sd(boot_vals_ok),
    n = n,
    k = k
  )
}


k_use <- df %>%
  filter(tsr_24 != 24) %>%
  count(tsr_24) %>%
  summarise(k = max(2L, floor(0.9 * min(n))), .groups = "drop") %>%
  pull(k)

k_use

k_site <- df %>%
  filter(tsr_24 != 24) %>%
  count(site, tsr_24) %>%
  group_by(site) %>%
  summarise(k = max(2L, floor(0.9 * min(n))), .groups = "drop")

k_site

var_fwd_t_eq <- df %>%
  filter(tsr_24 != 24) %>%
  group_by(tsr_24) %>%
  group_modify(~ cv_boot_std(.x$fwd_ha, k = k_use, nboot = 1000, replace = TRUE)) %>%
  ungroup() %>%
  arrange(tsr_24) %>%
  mutate(display = sprintf("%.1f%% (%.1f–%.1f%%)", boot_median, ci_lower, ci_upper))

tab_var_fwd_t_eq <- var_fwd_t_eq %>%
  select(tsr_24, display) %>%
  pivot_wider(names_from = tsr_24, values_from = display)

var_cwd_t_eq <- df %>%
  filter(tsr_24 != 24) %>%
  group_by(tsr_24) %>%
  group_modify(~ cv_boot_std(.x$cwd_vol_tot_ha, k = k_use, nboot = 1000, replace = TRUE)) %>%
  ungroup() %>%
  arrange(tsr_24) %>%
  mutate(display = sprintf("%.1f%% (%.1f–%.1f%%)", boot_median, ci_lower, ci_upper))

tab_var_cwd_t_eq <- var_cwd_t_eq %>%
  select(tsr_24, display) %>%
  pivot_wider(names_from = tsr_24, values_from = display)


var_fwd_s_eq <- df %>%
  filter(tsr_24 != 24) %>%
  left_join(k_site, by = "site") %>%
  group_by(site, tsr_24) %>%
  group_modify(~ cv_boot_std(.x$fwd_ha, k = first(.x$k), nboot = 1000, replace = TRUE)) %>%
  ungroup() %>%
  arrange(site, tsr_24) %>%
  mutate(display = sprintf("%.1f%% (%.1f–%.1f%%)", boot_median, ci_lower, ci_upper))

tab_var_fwd_s_eq <- var_fwd_s_eq %>%
  select(site, tsr_24, display) %>%
  pivot_wider(id_cols = site, names_from = tsr_24, values_from = display)

var_cwd_s_eq <- df %>%
  filter(tsr_24 != 24) %>%
  left_join(k_site, by = "site") %>%
  group_by(site, tsr_24) %>%
  group_modify(~ cv_boot_std(.x$cwd_vol_tot_ha, k = first(.x$k), nboot = 1000, replace = TRUE)) %>%
  ungroup() %>%
  arrange(site, tsr_24) %>%
  mutate(display = sprintf("%.1f%% (%.1f–%.1f%%)", boot_median, ci_lower, ci_upper))

tab_var_cwd_s_eq <- var_cwd_s_eq %>%
  select(site, tsr_24, display) %>%
  pivot_wider(id_cols = site, names_from = tsr_24, values_from = display)

# Inspect tables
tab_var_fwd_t_eq
tab_var_cwd_t_eq
tab_var_fwd_s_eq
tab_var_cwd_s_eq


# Figures  ####

unscale <- function(scaled, orig) {
  scaled * sd(orig, na.rm=TRUE) + mean(orig, na.rm=TRUE)
}

pc1 <- ggpredict(m3, terms=c("sc.tsr[all]")) %>%
  ggplot(aes(x = unscale(x, log2(df$tsr)), y = predicted)) +
  geom_jitter(df, mapping = aes(x = unscale(sc.tsr, log2(df$tsr)), y = lgcwd),
              size = 0.8, shape = 18, width = 0.1, height = 0.05, alpha=0.3) +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high), alpha = 0.2) +
  geom_line(size = 0.5,linetype = "dashed")+
  scale_y_continuous(name = "Coarse woody debris [m³/ha]",
                     breaks=c(-1.4, 0, 0.8, 1.4, 1.65), labels=c("0", "5", "10", "25", "50"))+ 
  scale_x_continuous(name="Tree species richness",
                     breaks=c(0,1,2,3,4), labels=c("1", "2", "4", "8", "16"))+ 
  theme_classic()+
  theme(
    legend.position = "none",
    axis.title = element_text(size = 7),
    axis.text = element_text(size = 7))
pc1

pc2 <- ggpredict(m3, terms=c("sc.prod[all]")) %>%
  ggplot(aes(x = unscale(x, df$productivity), y = predicted)) +
  geom_jitter(df, mapping = aes(x = unscale(sc.prod, df$productivity), y = lgcwd),
              size = 0.8, shape = 18, width = 0.1, height = 0.05, alpha=0.3) +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high), alpha = 0.2) +
  geom_line(size = 0.5) +
  scale_y_continuous(name = "Coarse woody debris [m³/ha]",
                     breaks=c(-1.4, 0, 0.8, 1.4, 1.65), labels=c("0", "5", "10", "25", "50"))+ 
  scale_x_continuous(name="Stand productivity [m³/ha/yr]")+ 
  theme_classic()+
  theme(
    legend.position = "none",
    axis.title = element_text(size = 7),
    axis.text = element_text(size = 7))

pc3 <- ggpredict(m3, terms=c("sc.wdens[all]")) %>%
  ggplot(aes(x = unscale(x, df$w_density), y = predicted)) +
  geom_jitter(df, mapping = aes(x = unscale(sc.wdens, df$w_density), y = lgcwd),
              size = 0.8, shape = 18, width = 0.1, height = 0.05, alpha=0.3) +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high), alpha = 0.2) +
  geom_line(size = 0.5) +
  scale_y_continuous(name = "Coarse woody debris [m³/ha]",
                     breaks=c(-1.4, 0, 0.8, 1.4, 1.65), labels=c("0", "5", "10", "25", "50"))+ 
  scale_x_continuous(name="Wood density [g/m³]", limits = c(0.35, 0.75), breaks = seq(0.35, 0.75, 0.1))+
  theme_classic()+
  theme(
    legend.position = "none",
    axis.title = element_text(size = 7),
    axis.text = element_text(size = 7))

pf1 <- ggpredict(m1, terms=c("sc.tsr[all]")) %>%
  ggplot(aes(x = unscale(x, log2(df$tsr)), y = predicted)) +
  geom_jitter(df, mapping = aes(x = unscale(sc.tsr, log2(df$tsr)), y = lgfwd),
              size = 0.8, shape = 18, width = 0.1, height = 0.05, alpha=0.3) +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high), alpha = 0.2) +
  geom_line(size = 0.5) +
  scale_y_continuous(name = "Fine woody debris [kg/ha]",
                     breaks=c(2,2.7, 3,3.7, 4), labels=c("100", "500", "1k", "5k", "10k"))+  
  scale_x_continuous(name="Tree species richness",
                     breaks=c(0,1,2,3,4), labels=c("1", "2", "4", "8", "16"))+  
  theme_classic() +
  theme(
    legend.position = "none",
    axis.title = element_text(size = 7),
    axis.text = element_text(size = 7))
pf1

pf2 <- ggpredict(m1, terms=c("sc.prod[all]")) %>%
  ggplot(aes(x = unscale(x, df$productivity), y = predicted)) +
  geom_jitter(df, mapping = aes(x = unscale(sc.prod, df$productivity), y = lgfwd),
              size = 0.8, shape = 18, width = 0.1, height = 0.05, alpha=0.3) +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high), alpha = 0.2) +
  geom_line(size = 0.5) +
  scale_y_continuous(name = "Fine woody debris [kg/ha]",
                     breaks=c(2,2.7, 3,3.7, 4), labels=c("100", "500", "1k", "5k", "10k"))+  
  scale_x_continuous(name="Stand productivity [m³/ha/yr]")+ 
  theme_classic() +
  theme(
    legend.position = "none",
    axis.title = element_text(size = 7),
    axis.text = element_text(size = 7))

pf3 <- ggpredict(m1, terms=c("sc.cc")) %>%
  ggplot(aes(x = unscale(x, df$canopy), y = predicted)) +
  geom_jitter(df, mapping = aes(x = unscale(sc.cc, df$canopy), y = lgfwd),
              size = 0.8, shape = 18, width = 0.1, height = 0.05, alpha=0.3) +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high), alpha = 0.2) +
  geom_line(size = 0.5) +
  scale_y_continuous(name = "Fine woody debris [kg/ha]",
                     breaks=c(2,2.7, 3,3.7, 4), labels=c("100", "500", "1k", "5k", "10k"))+  
  scale_x_continuous(name="Canopy cover [%]", breaks = seq(0, 1, 0.25))+ 
  theme_classic() +
  theme(
    legend.position = "none",
    axis.title = element_text(size = 7),
    axis.text = element_text(size = 7))

figure_3 <- plot_grid(plotlist = list(pc1, pf1, pc2, pf2,NULL, pf3),
                      labels = c("(a)", "(b)", "(c)","(d)","", "(e)"),
                      label_size = 8,
                      ncol = 2,
                      label_x = 0.25,   # (0 = sx, 1 =  dx)
                      label_y = 0.98)   # (0 = bottom, 1 = top)

figure_3

empty <- ggplot() +
  theme_void() +
  theme(
    plot.background = element_rect(fill = "white", colour = NA)
  )

figure_3 <- plot_grid(
  pc1, pf1, pc2, pf2, empty, pf3,
  labels = c("(a)", "(b)", "(c)", "(d)", "", "(e)"),
  label_size = 8,
  ncol = 2,
  label_x = 0.25,
  label_y = 0.98
)

figure_3
ggsave("Figure_3.tiff", figure_3, unit="mm", dpi=800, width=90, height=135)