## R code for generating Figures 4, 5, and 6
## Fixed effects plot, random intercept ridge plot, and random slope ridge plot

library(dplyr)
library(ggplot2)
library(ggridges)
library(tidyverse)
library(coda)
cbPalette <- viridisLite::viridis(5)

# specify a relative path containing the Rds and Rdata files
path <- "D:/Thermal_UAS/code_for_publishing"

# Random effect plots
modl <- list.files(path, pattern = ".*iter_lasso.*", full.names = TRUE) #load in pixel-level models
load(modl[[1]]); load(modl[[2]]); load(modl[[3]]);

int <- as.data.frame(mJune, pars = "^r_plantType") %>% 
  mutate(Month = "June") %>% select(-contains("CHM")) %>% 
  bind_rows(as.data.frame(mJuly, pars = "^r_plantType") %>% 
              mutate(Month = "July") %>% select(-contains("CHM"))) %>% 
  bind_rows(as.data.frame(mAugust, pars = "^r_plantType") %>% 
              mutate(Month = "August") %>% select(-contains("CHM")))

slo <- as.data.frame(mJune, pars = "^r_plantType") %>% 
  select(contains("CHM")) %>% mutate(Month = "June") %>% 
  bind_rows(as.data.frame(mJuly, pars = "^r_plantType") %>% 
              select(contains("CHM")) %>% mutate(Month = "July")) %>% 
  bind_rows(as.data.frame(mAugust, pars = "^r_plantType") %>% 
              select(contains("CHM")) %>% mutate(Month = "August"))

intercepts <- int %>% 
  pivot_longer(cols = 1:5) %>% 
  mutate(name = substr(name, 13, 16)) %>% 
  rename("EffectSize" = "value", "PlantType" = "name")

slopes <- slo %>% 
  pivot_longer(cols = 1:5) %>% 
  mutate(name = substr(name, 13, 16)) %>% 
  rename("EffectSize" = "value", "PlantType" = "name")

# find POD for each subspecies/month
int_pr <- intercepts %>% 
  mutate(PlantType = as.factor(PlantType), Month = as.factor(Month)) %>%
  group_by(Month, PlantType) %>%
  summarize(n_iter = n(),
            ngt0 = sum(EffectSize > 0),
            p_gt0 = ngt0 / n_iter) %>%
  select(Month, PlantType, p_gt0) %>%
  mutate(p_gt0 = ifelse(p_gt0 > 0.945 | p_gt0 < 0.055, 1, 0)) 

sl_pr <- slopes %>% 
  mutate(PlantType = as.factor(PlantType), Month = as.factor(Month)) %>%
  group_by(Month, PlantType) %>%
  summarize(n_iter = n(),
            ngt0 = sum(EffectSize > 0),
            p_gt0 = ngt0 / n_iter) %>%
  select(Month, PlantType, p_gt0) %>%
  mutate(p_gt0 = ifelse(p_gt0 > 0.945 | p_gt0 < 0.055, 1, 0)) 

# cols <- c(rgb(0,0,0,.5), rgb(1,1,1,.9)) # for significant
tcol <- c(viridisLite::plasma(16))
cols <- c(tcol[4], tcol[8], tcol[8], tcol[12], tcol[12])

# ridge-density plot ####
vline.dat <- data.frame(z=levels(as.factor(intercepts$Month)), vl=0) 

#### Figure 5 ####
intercepts %>% 
  mutate(PlantType = as.factor(PlantType), Month = as.factor(Month)) %>%
  left_join(int_pr, by = c("Month", "PlantType"), keep = FALSE) %>%
  mutate(PlantType = fct_recode(PlantType, "T-2x"="2n_T",
                                "T-4x"="4n_T","V-2x"="2n_V",
                                "V-4x"="4n_V","W-4x"="4n_W")) %>% 
  mutate(Month = factor(Month, levels = c("June", "July", "August")),
         PlantType = factor(PlantType, levels = c("T-2x", "T-4x", "V-2x", "V-4x", "W-4x"))) %>%
  group_by(PlantType, Month) %>%
  filter(EffectSize > coda::HPDinterval(as.mcmc(EffectSize), prob = 0.89)[1,][1],
         EffectSize < coda::HPDinterval(as.mcmc(EffectSize), prob = 0.89)[1,][2]) %>%
  ungroup() %>%
  ggplot(aes(y=PlantType, x=EffectSize, fill=as.factor(PlantType), height=..density..))+ 
  geom_density_ridges(scale=1, stat="density", panel_scaling = T, trim=1)+
  scale_fill_manual(values=rev(cols), name="Subspceis:cytotype") +
  scale_y_discrete(limits=rev) +
  labs(x=expression(paste("Relative foliar crown temperature (°C)")), y="Sagebrush Type")+
  guides(fill = FALSE) + theme_bw(base_size = 14) +
  geom_vline(aes(xintercept=vl), data=vline.dat, colour = "black", linetype = "dashed") +
  facet_grid(. ~ Month)

#### Figure 6 ####
slopes %>% 
  mutate(PlantType = as.factor(PlantType), Month = as.factor(Month)) %>%
  left_join(sl_pr, by = c("Month", "PlantType"), keep = FALSE) %>%
  mutate(PlantType = fct_recode(PlantType, "T-2x"="2n_T",
                                 "T-4x"="4n_T","V-2x"="2n_V",
                                 "V-4x"="4n_V","W-4x"="4n_W")) %>% 
  mutate(Month = factor(Month, levels = c("June", "July", "August")),
         PlantType = factor(PlantType, levels = c("T-2x", "T-4x", "V-2x", "V-4x", "W-4x"))) %>%
  group_by(PlantType, Month) %>%
  filter(EffectSize > coda::HPDinterval(as.mcmc(EffectSize), prob = 0.89)[1,][1],
         EffectSize < coda::HPDinterval(as.mcmc(EffectSize), prob = 0.89)[1,][2]) %>%
  ungroup() %>%
  ggplot(aes(y=PlantType,x=EffectSize,fill=as.factor(PlantType),height=..density..))+ #fill=as.factor(p_gt0),
  geom_density_ridges(scale=1,stat="density", panel_scaling = T,trim=1)+
  scale_fill_manual(values=rev(cols),name="Subspecies:cytotype")+
  scale_y_discrete(limits=rev) +
  labs(x=expression(paste("Change in crown temperature (°C) with pixel height")), y="Sagebrush Type")+
  guides(fill = FALSE) + theme_bw(base_size = 14) +
  geom_vline(aes(xintercept=vl), data=vline.dat, colour = "black", linetype = "dashed") +
  facet_grid(. ~ Month)

# Save active plot 
# ggsave("rSlope.png",
#      width=180,height=120,units=c("mm"),dpi=600,path = path)

#### Figure 4 ####
# fixed effect plot 
load(paste0(path, "mJune_final_6000iter_lasso_priors.RData"))
load(paste0(path, "mJuly_final_6000iter_lasso_priors.RData"))
load(paste0(path, "mAugust_final_10000iter_lasso_priors.RData"))

colmo <- c("cadetblue1","red3","pink3")

jul <- as.data.frame(mJuly) %>% 
  select("b_Intercept", "b_CHM_2019", "b_DistToPlantCenter", "b_Radiation_2019") %>% 
  rename(Intercept = b_Intercept, PixelHeight = b_CHM_2019, 
         DistanceToCenter = b_DistToPlantCenter, 
         Radiation = b_Radiation_2019) %>%
  pivot_longer(cols = 1:4) %>% 
  mutate(name = as.factor(name)) %>% 
  filter(name != "Intercept") %>%
  mutate(Month = as.factor("July"))

aug <- as.data.frame(mAugust) %>% 
  select("b_Intercept", "b_CHM_2019", "b_DistToPlantCenter", "b_Radiation_2019") %>% 
  rename(Intercept = b_Intercept, PixelHeight = b_CHM_2019, 
         DistanceToCenter = b_DistToPlantCenter, 
         Radiation = b_Radiation_2019) %>%
  pivot_longer(cols = 1:4) %>% 
  mutate(name = as.factor(name)) %>% 
  filter(name != "Intercept") %>%
  mutate(Month = as.factor("August"))


as.data.frame(mJune) %>% 
  select("b_Intercept", "b_CHM_2019", "b_DistToPlantCenter", "b_Radiation_2019") %>% 
  rename(Intercept = b_Intercept, PixelHeight = b_CHM_2019, 
         DistanceToCenter = b_DistToPlantCenter, 
         Radiation = b_Radiation_2019) %>%
  pivot_longer(cols = 1:4) %>% 
  mutate(name = as.factor(name)) %>% 
  filter(name != "Intercept") %>% 
  mutate(Month = as.factor("June")) %>% bind_rows(jul, aug) %>% 
  group_by(name, Month) %>%
  mutate(lower = coda::HPDinterval(as.mcmc(value), prob = 0.89)[1,][1],
         upper = coda::HPDinterval(as.mcmc(value), prob = 0.89)[1,][2],
         mean = mean(value)) %>%
  ungroup() %>% select(-value) %>%
  distinct(lower, upper, mean, .keep_all = TRUE) %>% 
  ggplot(aes(x = name, y = mean, colour = Month)) +
  geom_errorbar(aes(ymin=upper, ymax=lower), width=0.5, size=.5, position = position_dodge(width = .75)) + 
  geom_point(position = position_dodge(width = .75)) + 
  scale_colour_manual(values = rep(colmo, 3)) + 
  lims(y = c(-15, 5)) +
  geom_hline(yintercept = 0, linetype = "dashed", colour = "darkgray") +
  coord_flip() +
  labs(x=expression(paste("")), y=expression(paste("Relative foliar-crown temperature (°C)")))+
  theme_bw(base_size = 12) + theme(axis.text.y = element_text(size = 12)) 

# Save active plot 
# ggsave("fEff.png",
#        width=15, height=9, units="cm", dpi=300, path = path)