Minor changes to get rmds to run

Author: fseaton

.gitignore

@@ -1,4 +1,9 @@
 EIDC/
 .Rhistory
 Outputs/
-
+*.html
+*.RData
+*.rdb
+*.rdx
+*.png
+*__packages

01_Summary_Calculations.R

@@ -190,7 +190,8 @@ GRFLORA_AWI_SITE <- GRFLORA_ALL %>% ungroup() %>%
 
 
 BRLEAF_SUMMARY <- left_join(BRLEAF_SUMMARY, GRFLORA_AWI) %>%
-  mutate(AWI_RICH = tidyr::replace_na(AWI_RICH, 0))
+  mutate(AWI_RICH = tidyr::replace_na(AWI_RICH, 0)) %>%
+  left_join(AWI_sites, by = "SITE_NO")
 
 regions <- select(PLDATA7101, SITE, COUNTRY) %>% filter(COUNTRY != "") %>%
   left_join(AWI_sites, by = c("SITE" = "SITE_NO")) %>%
@@ -211,7 +212,8 @@ BRLEAF_SUMM_SITE <- BLEAF_META %>%
             ASHDIEBACK = sum(ASHDIEBACK, na.rm = TRUE), .groups = "drop") %>%
   mutate(ASHDIEBACK = ASHDIEBACK/NPLOT) %>%
   inner_join(GRFLORA_SITE) %>%
-  inner_join(GRFLORA_AWI_SITE)
+  inner_join(GRFLORA_AWI_SITE) %>%
+  left_join(AWI_sites, by = "SITE_NO")
 
 
 write.csv(BRLEAF_SUMM_SITE, "Outputs/Site level richness.csv",
@@ -459,6 +461,8 @@ Climate_data <- pivot_longer(Clim_data,
 ash_plots <- filter(DBH22, SITE_NO < 200 & 
                       SPECIES == "Fraxinus excelsior") %>% 
   select(SITE_NO, PLOT_NO) %>% unique() %>% mutate(ASHPRESENT = TRUE)
+write.csv(ash_plots, "Outputs/ash_plots.csv",
+          row.names = FALSE)
 ash_sites <- unique(ash_plots$SITE_NO)
 
 # deer risk

02_ChangeOverTime.Rmd

@@ -26,6 +26,7 @@ source("Helper_functions.R")
 
 
 ```{r data read}
+BRLEAF_SUMMARY <- read.csv("Outputs/Broadleaf_summary_metrics.csv")
 BRLEAF_SUMMARY <- BRLEAF_SUMMARY %>%
   mutate(YR = c("1971","2001","2022")[YEAR],
          YDAY = (DAYOFYEAR - 200)/30,
@@ -34,8 +35,9 @@ BRLEAF_SUMMARY <- BRLEAF_SUMMARY %>%
 regions <- read.csv("Outputs/regions.csv")
 ```
 
-# Check priors for count models
+# Count models
 
+### Check priors
 
 ```{r count model sample prior}
 mod_pr <- prior(normal(2,1), class = "b", coef = "YR1971") +
@@ -50,7 +52,6 @@ test_mod <- brm(SPECIES_RICH ~ -1 + YR + YDAY + (1|SITE:PLOT) + (YR|SITE),
                 data = BRLEAF_SUMMARY, family = "poisson",
                 prior = mod_pr, cores = 4, sample_prior = "only")
 plot(test_mod)
-pp_check(test_mod)
 pp_check(test_mod, "ecdf_overlay", ndraws = 20) +
   scale_x_continuous(limits = c(0,200))
 ```
@@ -286,7 +287,7 @@ basal_emm %>%
 # Cover models
 
 
-## Check priors for cover/gamma models
+### Check priors for cover/gamma models
 
 
 ```{r prop model sample prior}
@@ -538,7 +539,8 @@ awi_emm %>%
 # Binary models - Regeneration
 
 ```{r regen data prep}
-REGEN <- REGEN %>% select(-YEAR) %>% rename(YEAR = YR) %>%
+REGEN <- read.csv("Outputs/REGEN.csv")
+REGEN <- REGEN %>% select(-YEAR) %>% 
   inner_join(select(BRLEAF_SUMMARY, SITE_NO, PLOT_NO, YEAR, SITE, PLOT, YDAY, YR)) 
 
 ```
@@ -680,22 +682,15 @@ siterich_emm %>%
 
 # Site level AWI proportion
 
-```{r site level data prep awi}
-BRLEAF_SUMMARY_SITE <- GRFLORA_SITE %>%
-  mutate(YR = c("1971","2001","2022")[YEAR],
-         YDAY = (DAYOFYEAR - 200)/30,
-         SITE = as.character(SITE_NO)) %>%
-  left_join(select(BRLEAF_SUMMARY, SITE_NO, AWI_region) %>% distinct())
-```
 
-```{r count model sample prior v2}
+```{r prop model sample prior v2}
 mod_pr <- prior(normal(0,1), class = "b") +
   prior(student_t(5, 0, 1), class = "sd")
 ```
 
 ```{r awi site rich mod run}
 awi_site_rich_mod <- brm(AWI_RICH | trials(SPECIES_RICH) ~ -1 + YR + YDAY + (1|SITE) + (1|AWI_region),
-                         data = BRLEAF_SUMMARY_SITE, family = binomial(),
+                         data = BRLEAF_SUMM_SITE, family = binomial(),
                          prior = mod_pr, cores = 4, warmup = 2000, iter = 6000, thin = 4,
                          control = list(adapt_delta = 0.95),
                          file = paste0(model_loc, "SITERICH_AWI_BIN"))
@@ -728,22 +723,3 @@ awisite_emm %>%
 ```
 
 
-
-
-```{r emmeans site beta pairwise comparison table}
-sitebeta_emm <- emmeans(site_beta_mod, ~ YR)
-pairs(sitebeta_emm, type = "response") %>%
-  knitr::kable(digits = 3)
-```
-
-```{r emmeans site beta plot}
-sitebeta_emm %>%
-  gather_emmeans_draws() %>%
-  mutate(Year = as.numeric(as.character(YR)),
-         .value = exp(.value)) %>%
-  ggplot(aes(x = Year, y = .value)) +
-  stat_lineribbon(alpha = 1/4, fill = teal) +
-  theme(axis.text = element_text(size = 12), axis.title = element_text(size = 14)) +
-  labs(y = "Site level beta diversity")
-```
-

03_AshDieback.Rmd

@@ -32,6 +32,8 @@ source("Helper_functions.R")
 
 
 ```{r data read}
+ash_plots <- read.csv("Outputs/ash_plots.csv")
+ash_sites <- unique(ash_plots$SITE_NO)
 BRLEAF_SUMMARY <- read.csv("Outputs/Broadleaf_summary_metrics.csv")
 BRLEAF_SUMMARY <- BRLEAF_SUMMARY %>%
   full_join(ash_plots) %>%
@@ -782,15 +784,6 @@ siterich_emm %>%
 
 # Site level AWI proportion
 
-```{r site level awi data prep}
-BRLEAF_SUMMARY_SITE <- GRFLORA_SITE %>%
-  mutate(YR = c("1971","2001","2022")[YEAR],
-         YDAY = (DAYOFYEAR - 200)/30,
-         SITE = as.character(SITE_NO)) %>%
-  left_join(select(BRLEAF_SUMMARY, SITE_NO, AWI_region) %>% distinct()) %>%
-  filter(SITE_NO %in% ash_sites)
-
-```
 
 ```{r binomial model awi sample prior v2}
 mod_pr <- prior(normal(0,1), class = "b") +
@@ -799,7 +792,7 @@ mod_pr <- prior(normal(0,1), class = "b") +
 
 ```{r awi site rich mod run}
 awisite_mod <- brm(AWI_RICH | trials(SPECIES_RICH) ~ -1 + YR*ASHDIEBACK + YDAY + (1|SITE) + (1|AWI_region),
-                   data = BRLEAF_SUMMARY_SITE, family = binomial(),
+                   data = BRLEAF_SUMM_SITE, family = binomial(),
                    prior = mod_pr, cores = 4, warmup = 2000, iter = 6000, thin = 4,
                    control = list(adapt_delta = 0.95),
                    file = paste0(model_loc, "SITERICH_AWI_BIN"))

04_Deer.Rmd

@@ -41,6 +41,7 @@ BRLEAF_SUMMARY <- BRLEAF_SUMMARY %>%
          PLOT = paste(SITE_NO, PLOT_NO, sep = "_"),
          DEER = factor(DEER, c("Low","Moderate","High"),
                        ordered = TRUE))
+regions <- read.csv("Outputs/regions.csv")
 ```
 
 
@@ -802,17 +803,6 @@ siterich_emm %>%
 
 # Site level AWI proportion
 
-```{r site level awi data prep}
-GRFLORA_SITE <- read.csv("Outputs/Site level richness with AWI.csv")
-BRLEAF_SUMMARY_SITE <- GRFLORA_SITE %>%
-  mutate(YR = c("1971","2001","2022")[YEAR],
-         YDAY = (DAYOFYEAR - 200)/30,
-         SITE = as.character(SITE_NO)) %>%
-  left_join(select(BRLEAF_SUMMARY, SITE_NO, AWI_region) %>% distinct()) %>%
-  inner_join(deerrisk) %>%
-  mutate(DEER = factor(DEER, c("Low","Moderate","High"),
-                       ordered = TRUE))
-```
 
 ```{r binomial model sample prior site level}
 mod_pr <- prior(normal(0,1), class = "b") +
@@ -821,7 +811,7 @@ mod_pr <- prior(normal(0,1), class = "b") +
 
 ```{r awi site rich mod run}
 awisite_mod <- brm(AWI_RICH | trials(SPECIES_RICH) ~ -1 + YR*mo(DEER) + YDAY + (1|SITE) + (1|AWI_region),
-                   data = BRLEAF_SUMMARY_SITE, family = binomial(),
+                   data = BRLEAF_SUMM_SITE, family = binomial(),
                    prior = mod_pr, cores = 4, warmup = 2000, iter = 6000, thin = 4,
                    control = list(adapt_delta = 0.95),
                    file = paste0(model_loc, "SITERICH_AWI_BIN"))

05_AshDiebackandDeer.Rmd

@@ -35,6 +35,8 @@ source("Helper_functions.R")
 deerrisk <- read.csv("Metadata/DEER_RISK.csv") %>%
   mutate(DEER = ifelse(DEER == "None", "Low", DEER))
 regions <- read.csv("Outputs/regions.csv")
+ash_plots <- read.csv("Outputs/ash_plots.csv")
+ash_sites <- unique(ash_plots$SITE_NO)
 ```
 
 ```{r data read}
@@ -804,16 +806,6 @@ siterich_emm %>%
 
 # Site level AWI proportion
 
-```{r awi site level data prep}
-BRLEAF_SUMM_SITE <- BRLEAF_SUMM_SITE %>%
-  mutate(YR = c("1971","2001","2022")[YEAR],
-         YDAY = (DAYOFYEAR - 200)/30,
-         SITE = as.character(SITE_NO)) %>%
-  inner_join(select(BRLEAF_SUMMARY, SITE_NO, AWI_region) %>% distinct()) %>%
-  inner_join(deerrisk) %>% 
-  mutate(DEER = factor(DEER, c("Low","Moderate","High"),
-                       ordered = TRUE))
-```
 
 ```{r binomial model sample prior site level}
 mod_pr <- prior(normal(0,1), class = "b") +

06_Multivariate_Models.Rmd

@@ -22,32 +22,11 @@ theme_set(theme_classic())
 library(patchwork)
 teal <- "#0383a4"
 two_cols <- c("#999999","#0072B2","#E69F00")
+model_loc <- "Outputs/Models/Multivariate_Models/"
 ```
 
 ```{r, include = FALSE}
-par_summary <- function(x, digits = 3, ...){
-  summ <- summary(x)
-  fixef <- summ$fixed
-  rand1 <- summ$random[[1]]
-  rownames(rand1) <- paste0("SITE__",rownames(rand1))
-  if(length(summ$random)>1){
-    rand2 <- summ$random[[2]]
-    rownames(rand2) <- paste0("SITE:PLOT__",rownames(rand2))
-  } else{
-    rand2 <- rand1[0,]
-  }
-  cor_pars <- summ$cor_pars
-  spec_pars <- summ$spec_pars
-  if("mo" %in% names(summ)){
-    mo <- summ$mo
-    rownames(mo) <- paste0("smplx__", rownames(mo))
-    all <- do.call(rbind, list(fixef, rand1, rand2, spec_pars, mo))
-  } else{
-    all <- do.call(rbind, list(fixef, rand1, rand2, spec_pars))
-  }
-  all[,c("Bulk_ESS","Tail_ESS")] <- floor(all[,c("Bulk_ESS","Tail_ESS")])
-  knitr::kable(all, digits = digits, ...)
-}
+source("Helper_functions.R")
 ```
 
 
@@ -57,6 +36,9 @@ par_summary <- function(x, digits = 3, ...){
 deerrisk <- read.csv("Metadata/DEER_RISK.csv") %>%
   mutate(DEER = ifelse(DEER == "None", "Low", DEER)) %>%
   rename(SITE_NO = SITE_NUMBER)
+regions <- read.csv("Outputs/regions.csv")
+ash_plots <- read.csv("Outputs/ash_plots.csv")
+ash_sites <- unique(ash_plots$SITE_NO)
 ```
 
 ```{r data read}
@@ -88,17 +70,13 @@ BRLEAF_SUMMARY <- inner_join(BRLEAF_SUMMARY, Climate_data) %>%
          Winter_rainfall = (Winter_rainfall - 300)/100,
          Summer_tasmax = (Summer_tasmax - 19),
          Winter_tasmin = (Winter_tasmin - 1))
-REGEN <- read.csv("Outputs/REGEN_ASHDEER.csv") %>%
+REGEN <- read.csv("Outputs/REGEN.csv") %>%
   select(SITE_NO, PLOT_NO, YEAR, ends_with("REGEN")) %>% distinct()
 BRLEAF_SUMMARY <- inner_join(BRLEAF_SUMMARY, REGEN)
 ```
 
 # Multivariate Model
 
-```{r mult model loc, include = FALSE}
-model_loc <- "Outputs/Models/Multivariate_Models/"
-```
-
 ## Full model
 
 
@@ -195,6 +173,11 @@ pp_check(sem_allclim, resp = "SPECIESRICH", ndraws = 20, type = "ecdf_overlay")
 pp_check(sem_allclim, resp = "REGEN", ndraws = 20, type = "ecdf_overlay") 
 ```
 
+```{r}
+resid_bysite_mult(sem_allclim)
+```
+
+
 
 ```{r full model kfold}
 set.seed(1)
@@ -562,7 +545,7 @@ model_structure <-
 ```
 
 
-```{r run winter sr med model}
+```{r run winter regen med model}
 sem_win_regmed <- brm(model_structure, data = BRLEAF_SUMMARY, 
                      prior = mod_pr_nodeer +
                        prior(dirichlet(1), class = "simo", resp = "SPECIESRICH", coef = "moDEER1"),
@@ -571,7 +554,7 @@ sem_win_regmed <- brm(model_structure, data = BRLEAF_SUMMARY,
 par_summary(sem_win_regmed)
 ```
 
-```{r winter model kfold}
+```{r winter regen med model kfold}
 sem_win_regmed <- add_criterion(sem_win_regmed, "kfold", folds = folds,
                              file = paste0(model_loc, "SEM_WinterClim_RegenMed"))
 ```

Helper_functions.R

@@ -34,11 +34,40 @@ resid_bysite <- function(x, coldat = regions, ...){
     geom_hline(yintercept = 0) + labs(x = "Site", y = "Estimated Residual") +
     theme(axis.text.x = element_text(angle = 90, vjust = 0))
   p2 <- ggplot(plotdat, aes(y = Estimate)) +
-    stat_summary(aes(x = REGION), fun.data = "mean_se") +
+    stat_summary(aes(x = REGION), fun.data = "mean_se", 
+                 fun.args = list(mult = 1.96)) +
     geom_hline(yintercept = 0) + 
     labs(x = "Region", y = "Estimated Residual") + 
     theme(axis.text.x = element_text(angle = 90, vjust = 0.5))
 
   gridExtra::grid.arrange(p1,p2, nrow= 1, widths = c(2,1))
 
 }
+
+resid_bysite_mult <- function(x, coldat = regions, ...){
+  res <- resid(x)
+  var_names <- dimnames(res)[[3]]
+  pl_list <- lapply(var_names, function(i){
+    res_l <- res[,,i]
+    dat <- x$data
+    plotdat <- cbind(dat, res_l) %>%
+      left_join(mutate(coldat, SITE = as.character(SITE)),
+                by = "SITE")
+    p1 <- ggplot(plotdat, aes(y = Estimate)) +
+      geom_boxplot(aes(x = SITE)) +
+      facet_wrap(~REGION, scales = "free_x") +
+      geom_hline(yintercept = 0) + 
+      labs(x = "Site", y = "Estimated Residual",
+           title = i) +
+      theme(axis.text.x = element_text(angle = 90, vjust = 0))
+    p2 <- ggplot(plotdat, aes(y = Estimate)) +
+      stat_summary(aes(x = REGION), fun.data = "mean_se", 
+                   fun.args = list(mult = 1.96)) +
+      geom_hline(yintercept = 0) + 
+      labs(x = "Region", y = "Estimated Residual") + 
+      theme(axis.text.x = element_text(angle = 90, vjust = 0.5))
+    pa <- gridExtra::grid.arrange(p1,p2, nrow= 1, widths = c(2,1))
+  })
+  
+  
+}