--- title: <center><font size="6"><b>Multi-level inference to predict functional diversity</b></font></center> subtitle: <center><font size="4"><b>Coral reef functional diversity</b></font></center> author: <center><font size="4"><a href="https://marce10.github.io/">Marcelo Araya-Salas, Andrea Arriaga-Madrigal</a></font></center> date: "`r Sys.Date()`" toc: true toc-depth: 2 toc-location: left number-sections: true highlight-style: pygments format: html: df-print: kable code-fold: show code-tools: true css: qmd.css --- <!-- this code add line numbers to code blocks --> <!-- only works when code folding is not used in yaml (code_folding: show) --> ```{r set root directory, echo = FALSE} # set working directory knitr::opts_knit$set(root.dir = "..") ``` ```{r add link to github repo, echo = FALSE, results='asis'} # print link to github repo if any if (file.exists("./.git/config")){ config <- readLines("./.git/config") url <- grep("url", config, value = TRUE) url <- gsub("\\turl = |.git$", "", url) cat("\nSource code and data found at [", url, "](", url, ")", sep = "") } ``` ```{r setup style, echo = FALSE, message = FALSE, warning=FALSE} # options to customize chunk outputs knitr::opts_chunk$set( tidy.opts = list(width.cutoff = 65), tidy = TRUE, message = FALSE ) ``` <!-- skyblue box --> # Purpose {.unnumbered .unlisted} - Format data- Run stats# Analysis flowchart {.unnumbered .unlisted} ```{mermaid,eval = FALSE, echo = FALSE, fig.align = "center"} flowchart A[Read data] --> B(Format data) B --> C(Graphs) C --> D(Statistical analysis) D --> E(Model summary) style A fill:#44015466 style B fill:#3E4A894D style C fill:#26828E4D style D fill:#6DCD594D ``` # Load packages {.unnumbered .unlisted} ```{r load packages} # install/ load packages sketchy::load_packages( packages = c( "knitr", "formatR", "pbapply", "ggbeeswarm", "ggplot2", "viridis", "cluster", "dplyr", "magrittr", "reshape", "plyr", "ade4", "svMisc", "geometry", "ape", "zoo", "MuMIn", "phangorn", "brms", "lme4", "vegan", "gamlss", "gamlss.dist", "gamlss.add", "tidybayes", "ggridges", "HDInterval", "kableExtra", "ggcorrplot", "cmdstanr", "cowplot", "bayestestR", "warbleR", "cluster" ) ) ``` ```{r parameters functions and global options} theme_set(theme_classic(base_size = 20)) knitr::opts_knit$set(dpi = 100, fig.width = 16, fig.height = 8, warning = FALSE, root.dir = "..", message = FALSE) cols <- viridis(10, alpha = 0.7) ``` ```{r parameters functions and global options 2, eval = FALSE} source("~/Dropbox/Projects/reef_functional_diversity/scripts/assemblage_random_subsets.R") source("./scripts/multidimFD.R") # source("~/Dropbox/R_package_testing/warbleR/warbleR/R/pblapply_wrblr_int.R") pblapply_wrblr_int <- warbleR:::pblapply_wrblr_int ``` ```{r load and prepare data} # species_abundances <- readRDS("./data/species_abundances.rds") abundance_and_func_traits <- readRDS("~/Dropbox/Projects/eastern_tropical_pacific_herbivory/data_intermediate/fish/regional/herbfishes_regional_ecoregion_traits_biomass.rds") abundance_and_func_traits$site_id[abundance_and_func_traits$site_id == "nicaragua-toro_miscal/toro_marsella"] <- "nicaragua-toro_miscal-toro_marsella" # transect_areas1 <- readRDS("./data/transect_areas.RDS") transect_areas <- abundance_and_func_traits[!duplicated(paste(abundance_and_func_traits$ID_transect,abundance_and_func_traits$site_id)), c("ID_transect", "site_id", "area_uvc", "region", "date_y_m")] transect_areas$ID_transect_num <- 1:nrow(transect_areas) # Max number of sites a transect is found in (should be 1!) max(table(transect_areas$ID_transect, transect_areas$site_id)) uvc_sst_df <-readRDS("~/Dropbox/Projects/eastern_tropical_pacific_herbivory/results/biophy/uvc_sst_fishes_sites_data_frame.RDS") uvc_chl_df <-readRDS("~/Dropbox/Projects/eastern_tropical_pacific_herbivory/results/biophy/uvc_chl_fishes_sites_data_frame.RDS") predictors <- readRDS("~/Dropbox/Projects/eastern_tropical_pacific_herbivory/data_intermediate/predictors.RDS") uvc_sst_df$site_id[uvc_sst_df$site_id == "nicaragua-toro_miscal/toro_marsella"] <- uvc_chl_df$site_id[uvc_chl_df$site_id == "nicaragua-toro_miscal/toro_marsella"] <- predictors$site_id[predictors$site_id == "nicaragua-toro_miscal/toro_marsella"] <- "nicaragua-toro_miscal-toro_marsella" # these ones don't have coordinates (should be 4 sites from mexico) setdiff(transect_areas$site_id, predictors$site_id) transect_areas <- merge(transect_areas, predictors[!duplicated(predictors$site_id), ], by = "site_id") transect_areas$years_protected[is.na(transect_areas$years_protected)] <- 0 transect_areas$ID_transect_num <- 1:nrow(transect_areas) # load fish taxa to have "taxonomy" load("~/Dropbox/Projects/eastern_tropical_pacific_herbivory/data_intermediate/taxa/fish_taxa.rda") # sampling coverage sc <- readRDS("~/Dropbox/Projects/eastern_tropical_pacific_herbivory/results/mixed/sc.RDS") rownames(sc)[rownames(sc) == "nicaragua-toro_miscal/toro_marsella"] <- "nicaragua-toro_miscal-toro_marsella" # remove those with low coverage transect_areas <- transect_areas[transect_areas$site_id %in% rownames(sc)[sc$sc >= 0.95], ] ``` ```{r test with 1 run, eval = FALSE, include = FALSE} # get assemblages from random transect subsets random_assemblages <- assemblage_random_subsets(transc_metadata = transect_areas, sp_abund = abundance_and_func_traits, reps = 39, buffer = 0.1, parallel = 18, target.area = 400, site.specific.biomass = TRUE) assemblage_area <- unlist(pbapply::pblapply(1:nrow(random_assemblages), cl = 15, function(e){ w <- strsplit(rownames(random_assemblages)[e], "/", fixed = TRUE)[[1]] site_id <- w[1] ID_transect_nums <- strsplit(w[2], "-", fixed = TRUE)[[1]] aa <- sum(transect_areas$area_uvc[transect_areas$ID_transect_num %in% ID_transect_nums & transect_areas$site_id == site_id]) return(aa) })) # double check range of areas must be around target area +/- buffer range(assemblage_area) fds <- FD_random_assemblages(random_assemblages = random_assemblages, sp_cutoff = 3, prop_cutoff = 0.5, abundance_and_func_traits = abundance_and_func_traits, taxonomy = fish_taxa, transect_areas = transect_areas) fds_predictors <- add_predictors(X = fds, transect_areas = transect_areas, uvc_chl_df = uvc_chl_df, uvc_sst_df = uvc_sst_df, parallel = 1) # find best distribution # https://stats.stackexchange.com/questions/132652/how-to-determine-which-distribution-fits-my-data-best fit <- fitDist(fds_predictors$tax_distinctnes, trace = FALSE, try.gamlss = TRUE, type = "realAll") summary(fit) ``` # Create 100 random subsets ```{r run procedure several times, eval = FALSE} random_FD_list <- warbleR:::pblapply_wrblr_int(X = 1:110, FUN = function(x) { # print(x) random_assemblages <- assemblage_random_subsets(transc_metadata = transect_areas, sp_abund = abundance_and_func_traits, reps = 100, buffer = 0.1, parallel = 1, target.area = 400, pb = FALSE, seed = x, site.specific.biomass = TRUE) fds <- FD_random_assemblages(random_assemblages = random_assemblages, sp_cutoff = 3, prop_cutoff = 0.5, abundance_and_func_traits = abundance_and_func_traits, taxonomy = fish_taxa, transect_areas = transect_areas) fds_predictors <- try(add_predictors(X = fds, transect_areas = transect_areas, uvc_chl_df = uvc_chl_df, uvc_sst_df = uvc_sst_df, parallel = 1, pb = FALSE), silent = TRUE) return(fds_predictors) }, cl = 1) # remove errors random_FD_list <- random_FD_list[sapply(random_FD_list, class) == "data.frame"] # keep 100 random_FD_list <- random_FD_list[1:100] saveRDS(random_FD_list, "./data/100_replicated_random_assemblages_data_set_400m2_target_area.RDS") ``` # Check collinearity ## Predictors ```{r colinearity of predictors, fig.width= 12, fig.height=24} random_FD_list <- readRDS("./data/100_replicated_random_assemblages_data_set_400m2_target_area.RDS") # dev.off() # par(mfrow = c(4, 4)) # for(i in 1:16){ # print(head(random_FD_list[[i]])) ggcorrplots <- lapply(1:15, function(i){ cormat <- cor(random_FD_list[[i]][, c("sst_range", "chl_range", "gravity", "years_protected")]) ggcorrplot(cormat, type = "lower", outline.col = "black", lab=TRUE, ggtheme = ggplot2::theme_classic, colors = c("#6D9EC1", "white", "#E46726")) # # print(cp <- corrplot::corrplot.mixed(cormat, upper = "ellipse", lower.col = "black")) # } }) cowplot::plot_grid(plotlist = ggcorrplots, ncol = 3) ``` ## Responses ```{r colinearity of responses, fig.width= 12, fig.height=24} random_FD_list <- readRDS("./data/100_replicated_random_assemblages_data_set_400m2_target_area.RDS") # dev.off() # par(mfrow = c(4, 4)) # for(i in 1:16){ # print(head(random_FD_list[[i]])) ggcorrplots <- lapply(1:15, function(i){ cormat <- cor(random_FD_list[[i]][, c("FRic", "FDiv", "FEve", "FDis", "tax_distinctnes", "corrected_biomass", "density")]) ggcorrplot(cormat, type = "lower", outline.col = "black", lab=TRUE, ggtheme = ggplot2::theme_classic, colors = c("#6D9EC1", "white", "#E46726")) # # print(cp <- corrplot::corrplot.mixed(cormat, upper = "ellipse", lower.col = "black")) # } }) cowplot::plot_grid(plotlist = ggcorrplots, ncol = 3) ``` ```{r model averaging using MuMin proof of concept NOT USED, eval = FALSE, include = FALSE} # create 2 models with slightly different data # extract data Cement1 <- Cement[1:12,] fm1 <- lm(y ~ X1 + X2, data = Cement1) Cement1 <- Cement[2:13,] set.seed(103) Cement1$X2 <- sample(Cement1$X2) fm2 <- lm(y ~ X1 + X2, data = Cement1) fm1 fm2 options(na.action = "na.fail") dd <- dredge(fm2) # plot(dd) sdd <- subset(dd, delta < 5) mod.list <- list(fm1, fm2) # # class(sdd) <- c("model.selection", "data.frame") attrs <- attributes(sdd) model.call.names <- names(attrs$model.calls) attrs$model.calls <- replicate(n = length(attrs$model.calls), attrs$model.calls[[2]]) names(attrs$model.calls) <- model.call.names attrs$coefTables <- lapply(mod.list, function(x) { sm <- summary(x) tab <- sm$coefficients[, c("Estimate", "Std. Error")] tab <- cbind(tab, df = max(sm$df)) return(tab) }) # attrs$coefTables <- replicate(n = length(attrs$coefTables), attrs$coefTables[[2]], simplify = FALSE) sdd <- as.data.frame(sdd) sdd <- sdd[c(2,2), ] sdd$weight <- c(0.5, 0.5) attributes(sdd) <- attrs # model.avg(list(fm2, fm1), rank = "AICc") ma <- model.avg(sdd) summary(ma) coef(ma) confint(ma) # fm1 # fm2 mns <- t(sapply(mod.list, function(x) { sm <- summary(x) sm$coefficients[1:3,1] })) mns colMeans(mns) coef(ma) # summary(fm1) # summary(fm2) # mns <- t(sapply(mod.list, function(x) { # sm <- summary(x) # sm$coefficients[1:3,1] # })) # # confint.default(fm1) # summary(fm1) # # sm <- summary(fm1) # # sm$coefficients[,1] - 1.96*sm$coefficients[,2] ``` ```{r function for model averaging using MuMin NOT USED, eval = FALSE, include = FALSE} dredge_model_average <- function(mod.list){ dd <- dredge(mod.list[[1]]) # plot(dd) sdd <- subset(dd, which(1:nrow(dd) <= length(mod.list))) # mod.list <- list(fm1, fm2) # # class(sdd) <- c("model.selection", "data.frame") attrs <- attributes(sdd) model.call.names <- names(attrs$model.calls) attrs$model.calls <- replicate(n = length(attrs$model.calls), attrs$model.calls[[2]]) names(attrs$model.calls) <- model.call.names attrs$coefTables <- lapply(mod.list, function(x) { sm <- summary(x) tab <- sm$coefficients[, c("Estimate", "Std. Error")] # tab <- cbind(tab, df = max(sm$df)) tab <- cbind(tab, df = df.residual(x)) return(tab) }) # attrs$coefTables <- replicate(n = length(attrs$coefTables), attrs$coefTables[[2]], simplify = FALSE) sdd <- as.data.frame(sdd) sdd <- sdd[rep(2, length(mod.list)), ] sdd$weight <- rep(1 / length(mod.list), length(mod.list)) attrs$row.names <- 1:length(mod.list) attributes(sdd) <- attrs # model.avg(list(fm2, fm1), rank = "AICc") ma <- model.avg(sdd) # summary(ma) coefs <- coef(ma) conf_ints <- confint(ma) out_mat <- cbind(coefs, conf_ints) return(out_mat) } # create 2 models with slightly different data # extract data Cement1 <- Cement[1:12,] fm1 <- lm(y ~ X1 + X2, data = Cement1) Cement1 <- Cement[2:13,] set.seed(103) Cement1$X2 <- sample(Cement1$X2) fm2 <- lm(y ~ X1 + X2, data = Cement1) # 2 models dredge_model_average(mod.list = list(fm1, fm2)) set.seed(200) Cement1$X2 <- sample(Cement1$X2) fm3 <- lm(y ~ X1 + X2, data = Cement1) # 3 models dredge_model_average(mod.list = list(fm1, fm2, fm3)) # 5 models dredge_model_average(mod.list = list(fm1, fm2, fm3, fm3, fm2, fm1, fm1, fm2, fm1, fm1 , fm1)) ### using categorical variable iris1 <- iris fm1 <- lm(Sepal.Length ~ Sepal.Width + Species, data = iris1) iris1$Sepal.Width <- sample(iris1$Sepal.Width) fm2 <- lm(Sepal.Length ~ Sepal.Width + Species, data = iris1) dredge_model_average(mod.list = list(fm1, fm2)) ### using a random effect iris1 <- iris fm1 <- lmer(Sepal.Length ~ Sepal.Width + (1 | Species), data = iris1, REML = FALSE) iris1$Sepal.Width <- sample(iris1$Sepal.Width) fm2 <- lmer(Sepal.Length ~ Sepal.Width + (1 | Species), data = iris1, REML = FALSE) dredge_model_average(mod.list = list(fm1, fm1)) ### using a random effect iris1 <- iris iris1$Sepal.Length2 <- iris1$Sepal.Length iris1$Sepal.Length2 <- iris1$Sepal.Length2 / max(iris1$Sepal.Length2) fm1 <- lmer(Sepal.Length2 ~ Sepal.Width + (1 | Species), data = iris1, REML = FALSE) iris1$Sepal.Width <- sample(iris1$Sepal.Width) fm2 <- lmer(Sepal.Length ~ Sepal.Width + (1 | Species), data = iris1, REML = FALSE) dredge_model_average(mod.list = list(fm1, fm1)) ``` # Models ```{r list models, eval = TRUE} # modelos ## FDiv # modelos ## FDiv mod_all_noiucn_FDiv <- "FDiv ~ sst_range + chl_range + gravity + years_protected + geomorphology + (1|Ecoregion)" mod_all_noyears_FDiv <- "FDiv ~ sst_range + chl_range + gravity + iucn_cat + geomorphology + (1|Ecoregion)" mod_all_intr_FDiv <- "FDiv ~ sst_range + chl_range + gravity + years_protected : iucn_cat + geomorphology + (1|Ecoregion)" mod_envir_FDiv <- "FDiv ~ sst_range + chl_range + geomorphology + (1|Ecoregion)" mod_antrop_noiucn_FDiv <- "FDiv ~ gravity + years_protected + (1|Ecoregion)" mod_antrop_noyears_FDiv <- "FDiv ~ gravity + iucn_cat + (1|Ecoregion)" mod_antrop_intr_FDiv <- "FDiv ~ gravity + years_protected : iucn_cat + (1|Ecoregion)" ### FEve mod_all_noiucn_FEve <- "FEve ~ sst_range + chl_range + gravity + years_protected + geomorphology + (1|Ecoregion)" mod_all_noyears_FEve <- "FEve ~ sst_range + chl_range + gravity + iucn_cat + geomorphology + (1|Ecoregion)" mod_all_intr_FEve <- "FEve ~ sst_range + chl_range + gravity + years_protected : iucn_cat + geomorphology + (1|Ecoregion)" mod_envir_FEve <- "FEve ~ sst_range + chl_range + geomorphology + (1|Ecoregion)" mod_antrop_noiucn_FEve <- "FEve ~ gravity + years_protected + (1|Ecoregion)" mod_antrop_noyears_FEve <- "FEve ~ gravity + iucn_cat + (1|Ecoregion)" mod_antrop_intr_FEve <- "FEve ~ gravity + years_protected : iucn_cat + (1|Ecoregion)" ### FDis mod_all_noiucn_FDis <- "FDis ~ sst_range + chl_range + gravity + years_protected + geomorphology + (1|Ecoregion)" mod_all_noyears_FDis <- "FDis ~ sst_range + chl_range + gravity + iucn_cat + geomorphology + (1|Ecoregion)" mod_all_intr_FDis <- "FDis ~ sst_range + chl_range + gravity + years_protected : iucn_cat + geomorphology + (1|Ecoregion)" mod_envir_FDis <- "FDis ~ sst_range + chl_range + geomorphology + (1|Ecoregion)" mod_antrop_noiucn_FDis <- "FDis ~ gravity + years_protected + (1|Ecoregion)" mod_antrop_noyears_FDis <- "FDis ~ gravity + iucn_cat + (1|Ecoregion)" mod_antrop_intr_FDis <- "FDis ~ gravity + years_protected : iucn_cat + (1|Ecoregion)" ### tax dis mod_all_noiucn_taxdis <- "tax_distinctnes ~ sst_range + chl_range + gravity + years_protected + geomorphology + (1|Ecoregion)" mod_all_noyears_taxdis <- "tax_distinctnes ~ sst_range + chl_range + gravity + iucn_cat + geomorphology + (1|Ecoregion)" mod_all_intr_taxdis <- "tax_distinctnes ~ sst_range + chl_range + gravity + years_protected : iucn_cat + geomorphology + (1|Ecoregion)" mod_envir_taxdis <- "tax_distinctnes ~ sst_range + chl_range + geomorphology + (1|Ecoregion)" mod_antrop_noiucn_taxdis <- "tax_distinctnes ~ gravity + years_protected + (1|Ecoregion)" mod_antrop_noyears_taxdis <- "tax_distinctnes ~ gravity + iucn_cat + (1|Ecoregion)" mod_antrop_intr_taxdis <- "tax_distinctnes ~ gravity + years_protected : iucn_cat + (1|Ecoregion)" ### FRiq mod_all_noiucn_FRic <- "FRic ~ sst_range + chl_range + gravity + years_protected + geomorphology + (1|Ecoregion)" mod_all_noyears_FRic <- "FRic ~ sst_range + chl_range + gravity + iucn_cat + geomorphology + (1|Ecoregion)" mod_all_intr_FRic <- "FRic ~ sst_range + chl_range + gravity + years_protected : iucn_cat + geomorphology + (1|Ecoregion)" mod_envir_FRic <- "FRic ~ sst_range + chl_range + geomorphology + (1|Ecoregion)" mod_antrop_noiucn_FRic <- "FRic ~ gravity + years_protected + (1|Ecoregion)" mod_antrop_noyears_FRic <- "FRic ~ gravity + iucn_cat + (1|Ecoregion)" mod_antrop_intr_FRic <- "FRic ~ gravity + years_protected : iucn_cat + (1|Ecoregion)" ### Density mod_all_noiucn_dens <- "density ~ sst_range + chl_range + gravity + years_protected + geomorphology + (1|Ecoregion)" mod_all_noyears_dens <- "density ~ sst_range + chl_range + gravity + iucn_cat + geomorphology + (1|Ecoregion)" mod_all_intr_dens <- "density ~ sst_range + chl_range + gravity + years_protected : iucn_cat + geomorphology + (1|Ecoregion)" mod_envir_dens <- "density ~ sst_range + chl_range + geomorphology + (1|Ecoregion)" mod_antrop_noiucn_dens <- "density ~ gravity + years_protected + (1|Ecoregion)" mod_antrop_noyears_dens <- "density ~ gravity + iucn_cat + (1|Ecoregion)" mod_antrop_intr_dens <- "density ~ gravity + years_protected : iucn_cat + (1|Ecoregion)" ### Biomass mod_all_noiucn_biom <- "corrected_biomass ~ sst_range + chl_range + gravity + years_protected + geomorphology + (1|Ecoregion)" mod_all_noyears_biom <- "corrected_biomass ~ sst_range + chl_range + gravity + iucn_cat + geomorphology + (1|Ecoregion)" mod_all_intr_biom <- "corrected_biomass ~ sst_range + chl_range + gravity + years_protected : iucn_cat + geomorphology + (1|Ecoregion)" mod_envir_biom <- "corrected_biomass ~ sst_range + chl_range + geomorphology + (1|Ecoregion)" mod_antrop_noiucn_biom <- "corrected_biomass ~ gravity + years_protected + (1|Ecoregion)" mod_antrop_noyears_biom <- "corrected_biomass ~ gravity + iucn_cat + (1|Ecoregion)" mod_antrop_intr_biom <- "corrected_biomass ~ gravity + years_protected : iucn_cat + (1|Ecoregion)" ### Richness mod_all_noiucn_rich <- "Nb_sp ~ sst_range + chl_range + gravity + years_protected + geomorphology + (1|Ecoregion)" mod_all_noyears_rich <- "Nb_sp ~ sst_range + chl_range + gravity + iucn_cat + geomorphology + (1|Ecoregion)" mod_all_intr_rich <- "Nb_sp ~ sst_range + chl_range + gravity + years_protected : iucn_cat + geomorphology + (1|Ecoregion)" mod_envir_rich <- "Nb_sp ~ sst_range + chl_range + geomorphology + (1|Ecoregion)" mod_antrop_noiucn_rich <- "Nb_sp ~ gravity + years_protected + (1|Ecoregion)" mod_antrop_noyears_rich <- "Nb_sp ~ gravity + iucn_cat + (1|Ecoregion)" mod_antrop_intr_rich <- "Nb_sp ~ gravity + years_protected : iucn_cat + (1|Ecoregion)" models <- list( mod_all_noiucn_FDiv = mod_all_noiucn_FDiv, mod_all_noyears_FDiv = mod_all_noyears_FDiv, mod_all_intr_FDiv = mod_all_intr_FDiv, mod_envir_FDiv = mod_envir_FDiv, mod_antrop_noiucn_FDiv = mod_antrop_noiucn_FDiv, mod_antrop_noyears_FDiv = mod_antrop_noyears_FDiv, mod_antrop_intr_FDiv = mod_antrop_intr_FDiv, mod_all_noiucn_FEve = mod_all_noiucn_FEve, mod_all_noyears_FEve = mod_all_noyears_FEve, mod_all_intr_FEve = mod_all_intr_FEve, mod_envir_FEve = mod_envir_FEve, mod_antrop_noiucn_FEve = mod_antrop_noiucn_FEve, mod_antrop_noyears_FEve = mod_antrop_noyears_FEve, mod_antrop_intr_FEve = mod_antrop_intr_FEve, mod_all_noiucn_FDis = mod_all_noiucn_FDis, mod_all_noyears_FDis = mod_all_noyears_FDis, mod_all_intr_FDis = mod_all_intr_FDis, mod_envir_FDis = mod_envir_FDis, mod_antrop_noiucn_FDis = mod_antrop_noiucn_FDis, mod_antrop_noyears_FDis = mod_antrop_noyears_FDis, mod_antrop_intr_FDis = mod_antrop_intr_FDis, mod_all_noiucn_taxdis = mod_all_noiucn_taxdis, mod_all_noyears_taxdis = mod_all_noyears_taxdis, mod_all_intr_taxdis = mod_all_intr_taxdis, mod_envir_taxdis = mod_envir_taxdis, mod_antrop_noiucn_taxdis = mod_antrop_noiucn_taxdis, mod_antrop_noyears_taxdis = mod_antrop_noyears_taxdis, mod_antrop_intr_taxdis = mod_antrop_intr_taxdis, mod_all_noiucn_FRic = mod_all_noiucn_FRic, mod_all_noyears_FRic = mod_all_noyears_FRic, mod_all_intr_FRic = mod_all_intr_FRic, mod_envir_FRic = mod_envir_FRic, mod_antrop_noiucn_FRic = mod_antrop_noiucn_FRic, mod_antrop_noyears_FRic = mod_antrop_noyears_FRic, mod_antrop_intr_FRic = mod_antrop_intr_FRic, mod_all_noiucn_dens = mod_all_noiucn_dens, mod_all_noyears_dens = mod_all_noyears_dens, mod_all_intr_dens = mod_all_intr_dens, mod_envir_dens = mod_envir_dens, mod_antrop_noiucn_dens = mod_antrop_noiucn_dens, mod_antrop_noyears_dens = mod_antrop_noyears_dens, mod_antrop_intr_dens = mod_antrop_intr_dens, mod_all_noiucn_biom = mod_all_noiucn_biom, mod_all_noyears_biom = mod_all_noyears_biom, mod_all_intr_biom = mod_all_intr_biom, mod_envir_biom = mod_envir_biom, mod_antrop_noiucn_biom = mod_antrop_noiucn_biom, mod_antrop_noyears_biom = mod_antrop_noyears_biom, mod_antrop_intr_biom = mod_antrop_intr_biom, mod_all_noiucn_rich = mod_all_noiucn_rich, mod_all_noyears_rich = mod_all_noyears_rich, mod_all_intr_rich = mod_all_intr_rich, mod_envir_rich = mod_envir_rich, mod_antrop_noiucn_rich = mod_antrop_noiucn_rich, mod_antrop_noyears_rich = mod_antrop_noyears_rich, mod_antrop_intr_rich = mod_antrop_intr_rich ) # temporarily remove those including "iucn_cat" models <- models[grep(pattern = "iucn", x = models, invert = TRUE)] print(models) ``` ## Bayesian mixed models - 12 models for each data subset```{r statistical analysis BAYESIAN brms no spatial autocorrelation propagating uncertainty in FD, eval = FALSE} random_FD_list <- readRDS("./data/100_replicated_random_assemblages_data_set_400m2_target_area.RDS") # check sites on all data sets sites <- unlist(sapply(random_FD_list, function(x) unique(x$site_id))) sites_count <- table(sites) shared_sites <- names(sites_count)[sites_count == length(random_FD_list)] length(shared_sites) # predictors to scale vars <- c("sst_range", "chl_range", "gravity", "years_protected") iters <- 30000 # z-transform predictors and keep only shared sites random_FD_list <- lapply(random_FD_list, function(W){ # scale parameters for(i in vars) W[, i] <- scale(W[, i]) # remove non-shared sites W <- W[W$site_id %in% shared_sites, ] #make taxonomic distinctness a proportion W$tax_distinctnes <- W$tax_distinctnes / (max(W$tax_distinctnes) * 1.0001) # set base level for iucn categories W$iucn_cat <- as.character(W$iucn_cat) return(W) }) priors <- c(set_prior("normal(0, 1.5)", class = "Intercept"), set_prior("normal(0, 2)", class = "b")) # loop to run bayesian mcmc glmmm in brms (using stan) controlling for spatial autocorrelation # remove iucn models <- grep("iucn", models, value = TRUE, invert = TRUE) shuffled_models <- sample(1:length(models)) brms_model_out <- warbleR:::pblapply_wrblr_int(shuffled_models, function(u){ x <- models[[u]] print(x) print(names(models)[u]) rds_name <- file.path("./data/statistical_models_2024", paste0(names(models)[u], "-10-data-sets_brm_multiple.RDS")) # check Rhats # if (file.exists(rds_name)){ # mod <- readRDS(rds_name) # rhats <- rhat(mod[[1]][[1]]) # rhats <- rhats[-length(rhats)] # iters <- attr(mod[[1]][[1]]$fit, "stan_args")[[1]]$iter + 10000 # # go <- if (any(rhats > 1.05)) TRUE else # FALSE # } else go <- TRUE # # if (go) { if (TRUE) { # get priors for model sub_priors <- priors[sapply(priors$coef, function(x) grepl(x, as.character(x)[3])), ] # mod <- try(brm(formula = x, iter = 0000, thin = 1, warmup = 5000, # data = random_FD_list[[1]], # family = Beta(), # control = list(adapt_delta = 0.999), # silent = 2, # chains = 1, # prior = sub_priors # ), silent = TRUE) mod <- try(brm_multiple(formula = x, iter = iters, thin = 1, data = random_FD_list[1:10], family = Beta(), cores = 10, control = list(adapt_delta = 0.999), silent = 2, chains = 1, combine = FALSE, prior = sub_priors), silent = TRUE) # add loo if (!is(mod, "try-error")) for (i in 1:length(mod)) mod[[i]] <- add_criterion(mod[[i]], c("loo")) if (!is(mod, "try-error")) saveRDS(mod, file = rds_name) rm(mod) } return(NULL) }) ``` ```{r statistical analysis BAYESIAN brms siteid as random effect no spatial autocorrelation propagating uncertainty in FD, eval = FALSE} random_FD_list <- readRDS("./data/100_replicated_random_assemblages_data_set_400m2_target_area.RDS") random_FD_df <- do.call(rbind, random_FD_list) # predictors to scale vars <- c("sst_range", "chl_range", "gravity", "years_protected") iters <- 50000 for (i in vars) random_FD_df[, i] <- scale(random_FD_df[, i]) #make taxonomic distinctness, Nb_sp density and corrected_biomass a proportion random_FD_df$tax_distinctnes <- random_FD_df$tax_distinctnes / (max(random_FD_df$tax_distinctnes) * 1.0001) random_FD_df$Nb_sp <- random_FD_df$Nb_sp / (max(random_FD_df$Nb_sp) * 1.0001) random_FD_df$density <- random_FD_df$density / (max(random_FD_df$density) * 1.0001) random_FD_df$corrected_biomass <- random_FD_df$corrected_biomass / (max(random_FD_df$corrected_biomass) * 1.0001) priors <- c( set_prior("normal(0, 1.5)", class = "Intercept"), set_prior("normal(0, 2)", class = "b") ) # loop to run bayesian mcmc glmmm in brms (using stan) controlling for spatial autocorrelation # remove iucn models <- grep("iucn", models, value = TRUE, invert = TRUE) shuffled_models <- sample(1:length(models)) brms_model_out <- warbleR:::pblapply_wrblr_int(shuffled_models, function(u) { x <- models[[u]] print(x) print(names(models)[u]) rds_name <- file.path("./data/statistical_models_2024", paste0(names(models)[u], "-single_combined_100_data_sets.RDS")) # check Rhats if (!file.exists(rds_name)){ # mod <- readRDS(rds_name) # rhats <- rhat(mod[[1]][[1]]) # rhats <- rhats[-length(rhats)] # iters <- attr(mod[[1]][[1]]$fit, "stan_args")[[1]]$iter + 10000 # # go <- if (any(rhats > 1.05)) TRUE else # FALSE # } else go <- TRUE # # if (go) { # if (TRUE) { # get priors for model sub_priors <- priors[sapply(priors$coef, function(x) grepl(x, as.character(x)[3])),] mod <- try(brm( formula = paste(x, "+ (1|site_id)"), iter = iters, thin = 1, data = random_FD_df, family = Beta(), control = list(adapt_delta = 0.999), silent = 2, chains = 2, cores = 2, backend = "cmdstanr", threads = threading(2), prior = sub_priors, save_pars = save_pars(all = TRUE) ), silent = TRUE) if (!is(mod, "try-error")) mod2 <- try(add_criterion(mod, c("loo")), silent = TRUE) # if (!is(mod2, "try-error")) # mod <- mod2 if (!is(mod, "try-error")) saveRDS(mod, file = rds_name) rm(mod) } return(NULL) }) ``` ```{r print results statistical analysis BAYESIAN brms no spatial autocorrelation, eval = FALSE} brms_model_rds <- list.files("./data/statistical_models_2024", pattern = "combined_100_data_sets.RDS", full.names = TRUE) brms_model_rds <- brms_model_rds[brms_model_rds != "./data/statistical_models/models_combined_100_data_sets.RDS" ] # temporarily remove iucn # brms_model_rds <- grep(pattern = "iucn", x = brms_model_rds, invert = TRUE, value = TRUE) combined_mod_results <- pblapply(brms_model_rds, cl = 1, function(e) { mod <- readRDS(e) # comb_mod <- combine_models(mlist = mods, check_data = FALSE) # mods[[1]]$criteria$loo$estimates # print(x$formula) summ <- summary(mod)$fixed fit <- mod$fit betas <- grep("^b_", names(fit@sim$samples[[1]]), value = TRUE) hdis <- t(sapply(betas, function(y) hdi(fit@sim$samples[[1]][[y]])) ) summ_table <- data.frame(summ, hdis, iterations = attr(fit, "stan_args")[[1]]$iter, chains = length(attr(fit, "stan_args"))) summ_table <- summ_table[rownames(summ_table) != "Intercept", c("Estimate", "Rhat", "Bulk_ESS", "CI_low", "CI_high", "iterations", "chains")] summ_table <- as.data.frame(summ_table) summ_table$Rhat <- round(summ_table$Rhat, digits = 3) summ_table$CI_low <- round(unlist(summ_table$CI_low), digits = 3) summ_table$CI_high <- round(unlist(summ_table$CI_high), digits = 3) out <- lapply(betas, function(y) data.frame(variable = y, value = sort(fit@sim$samples[[1]][[y]], decreasing = FALSE))) posteriors <- do.call(rbind, out) posteriors <- posteriors[posteriors$variable != "b_Intercept", ] gg <- ggplot(data = posteriors, aes(y = variable, x = value, fill = stat(quantile))) + geom_density_ridges_gradient(quantile_lines = TRUE, quantile_fun = HDInterval::hdi, vline_linetype = 2) + theme_classic() + xlim(c(min(summ_table[ , "CI_low"]), max(summ_table[ , "CI_high"]))) + geom_vline(xintercept = 0, col = "red") + scale_fill_manual(values = c("transparent", "lightblue", "transparent"), guide = "none") + ggtitle(mod$formula) # print(gg) summ_table$Rhat <- ifelse(summ_table$Rhat > 1.05, cell_spec(summ_table$Rhat, "html", color ="white", background = "red", bold = TRUE, font_size = 12), cell_spec(summ_table$Rhat, "html")) signif <- summ_table[,"CI_low"] * summ_table[,"CI_high"] > 0 # summ_table$CI_low <- ifelse(signif, cell_spec(summ_table$CI_low, "html", color ="black", background = adjustcolor(cols[9], alpha.f = 0.3), bold = TRUE, font_size = 12), cell_spec(summ_table$CI_low, "html")) # summ_table$CI_high <- ifelse(signif, cell_spec(summ_table$CI_high, "html", color ="black", background = adjustcolor(cols[9], alpha.f = 0.3), bold = TRUE, font_size = 12), cell_spec(summ_table$CI_high, "html")) df1 <- kbl(summ_table, row.names = TRUE, escape = FALSE, format = "html", digits = 3) df1 <- row_spec(kable_input = df1,row = which(summ_table$CI_low * summ_table$CI_high > 0), background = adjustcolor(cols[9], alpha.f = 0.3)) df1 <- kable_styling(df1, bootstrap_options = c("striped", "hover", "condensed", "responsive"), full_width = FALSE, font_size = 12) # print(df1) results <- list(table = df1, gg = gg) return(results) return(cm) } ) names(combined_mod_results) <- gsub("-single_combined_100_data_sets.RDS", "", basename(brms_model_rds)) saveRDS(combined_mod_results, "./data/statistical_models_2024/models_combined_100_data_sets_2024.RDS") ``` # Results for a single combined data set per model (100 subsets, site as random factor) ```{r plots and tables single model site random, eval = F, fig.width= 12, results = 'asis', message=FALSE} models_combined_100_data_set <- readRDS("./data/statistical_models_2024/models_combined_100_data_sets_2024.RDS") for(i in 1:length(models_combined_100_data_set)){ print(names(models_combined_100_data_set)[i]) print(models_combined_100_data_set[[i]]$table) print(models_combined_100_data_set[[i]]$gg) } ``` ## Model comparison ```{r calculate model comparison 100 data sets and site id random, eval = FALSE} brms_model_rds <- list.files(path = "./data/statistical_models_2024/", pattern = "-single_combined_100_data_sets.RDS", full.names = TRUE) looic_list <- pblapply(brms_model_rds, cl = 1, function(e) { mod <- readRDS(e) if (is.null(mod$criteria$loo$estimates)) mod <- try(add_criterion(mod, c("loo")), silent = TRUE) loos <- mod$criteria$loo$estimates rm(mod) out_df <- data.frame( models = gsub("-single_combined_100_data_sets.RDS", "", basename(e)), looic = loos[3, 1], SE = loos[3, 2] ) return(out_df) }) looic_df <- do.call(rbind, looic_list) saveRDS(looic_df, "./data/looic_by_model_site_id_as_random_effect_2024.RDS") ``` ```{r model comparison 100 data set site id random, eval = FALSE, fig.width= 12, results = 'asis', message=FALSE} looic_df <- readRDS("./data/looic_by_model_site_id_as_random_effect_2024.RDS") brms_model_rds <- list.files("./data/statistical_models_2024/", pattern = "-single_combined_100_data_sets.RDS", full.names = TRUE) mod_groups <- c("FDis", "FDiv", "FEve", "taxdis", "biom", "dens", "rich") loo_comparisons <- lapply(mod_groups, function(i){ print(i) sub_looic <- looic_df[grep(i, looic_df$models), ] sub_looic$delta.looic <- sub_looic$looic - min(sub_looic$looic) mods <- lapply(file.path("./data/statistical_models_2024/", paste0(sub_looic$models, "-single_combined_100_data_sets.RDS")), FUN = readRDS) weights <- brms::loo_model_weights(mods[[1]], mods[[2]], mods[[3]], model_names = sub_looic$models, ndraws = 5000) rm(mods) sub_looic$looic.weights <- sapply(sub_looic$models, function(x) as.numeric(weights[names(weights) == x])) sub_looic <- sub_looic[order(sub_looic$delta.looic, decreasing = TRUE), ] sub_looic$cum.weight <- cumsum(sub_looic$looic.weights) sub_looic <- sub_looic[order(sub_looic$delta.looic, decreasing = FALSE), ] df1 <- knitr::kable(sub_looic, row.names = TRUE, escape = FALSE, format = "html", digits = 3, caption = i) df1 <- row_spec(kable_input = df1,row = 1:which(sub_looic$cum.weight >= 0.95)[1], background = adjustcolor(cols[9], alpha.f = 0.3)) df1 <- kable_styling(df1, bootstrap_options = c("striped", "hover", "condensed", "responsive"), full_width = FALSE, font_size = 12) return(df1) }) names(loo_comparisons) <- mod_groups saveRDS(loo_comparisons, "./data/statistical_models_2024/models_comparison_results_100_data_sets_2024.RDS") ``` ```{r model comparison 100 data set site id random print, eval = TRUE, fig.width= 12, results = 'asis', message=FALSE} loo_comparisons <- readRDS("./data/statistical_models_2024/models_comparison_results_100_data_sets_2024.RDS") for(i in 1:length(loo_comparisons)){ print(names(loo_comparisons)[i]) print(loo_comparisons[[i]]) } ``` ## Model averaging ```{r model averaging, eval = FALSE, echo = TRUE} brms_model_rds <- list.files("./data/statistical_models_2024", pattern = "-single_combined_100_data_sets.RDS", full.names = TRUE) average_model_list <- pblapply(c("FDis", "FDiv", "FEve", "taxdis", "biom", "dens", "rich"), function(x){ mods <- lapply(grep(x, brms_model_rds, value = TRUE), readRDS) params <- lapply(mods, function(x) rownames(summary(x)$fixed)[-1]) combs <- combn(1:length(mods), 2) post_avrg_list <- lapply(1:ncol(combs), function(y){ prm <- paste0("b_", intersect(params[[combs[1, y]]], params[[combs[2, y]]])) if (!prm[1] == "b_") # calculate average posteriors and model post_avrg <- posterior_average(mods[[combs[1, y]]], mods[[combs[1, y]]], weights = "loo", pars = prm) else post_avrg <- NULL return(post_avrg) }) post_avrg_list <- post_avrg_list[sapply(post_avrg_list, is.data.frame)] post_avrgs <- do.call(cbind, post_avrg_list) average_model <- describe_posterior(as.data.frame(post_avrgs), ci_method = "HDI") average_model <- as.data.frame(average_model) average_model_raw <- average_model average_model$CI_low <- round(average_model$CI_low, digits = 3) average_model$CI_high <- round(average_model$CI_high, digits = 3) signif <- average_model[, "CI_low"] * average_model[, "CI_high"] > 0 # average_model$CI_low <- ifelse(signif, cell_spec(average_model$CI_low, "html", color ="black", background = adjustcolor(cols[9], alpha.f = 0.3), bold = TRUE, font_size = 12), cell_spec(average_model$CI_low, "html")) # # average_model$CI_high <- ifelse(signif, cell_spec(average_model$CI_high, "html", color ="black", background = adjustcolor(cols[9], alpha.f = 0.3), bold = TRUE, font_size = 12), cell_spec(average_model$CI_high, "html")) # df1 <- kbl(average_model, row.names = TRUE, escape = FALSE, format = "html", digits = 3) df1 <- row_spec(kable_input = df1,row = which(signif), background = adjustcolor(cols[9], alpha.f = 0.3)) df1 <- kable_styling(df1, bootstrap_options = c("striped", "hover", "condensed", "responsive"), full_width = FALSE, font_size = 12) average_model_raw$response <- x out <- list(average_model = average_model_raw, table = df1) return(out) }) names(average_model_list) <- c("FDis", "FDiv", "FEve", "taxdis", "biom", "dens", "rich") averaged_model_table <- do.call(rbind, lapply(average_model_list, "[[", 'average_model')) saveRDS(averaged_model_table, "./data/statistical_models_2024/averaged_model_table.RDS") saveRDS(average_model_list, "./data/statistical_models_2024/average_models.RDS") ``` ```{r average model results, eval = TRUE, fig.width= 12, results = 'asis', message=FALSE} average_model_list <- readRDS("./data/statistical_models_2024/average_models.RDS") for(i in 1:length(average_model_list)){ print(names(average_model_list)[i]) print(average_model_list[[i]]$table) } ``` # Results for a single data set per model ```{r plots and tables single model, eval = FALSE, fig.width= 12, results = 'asis', message=FALSE} models_1_data_set <- readRDS("./data/statistical_models/combined_models_results_1_data_set.RDS") for(i in 1:length(models_1_data_set)){ print(names(models_1_data_set)[i]) print(models_1_data_set[[i]]$table) print(models_1_data_set[[i]]$gg) } ``` # Results for 10 data set models (propagating FD uncertainty) ```{r plots and tables 10 data set models, eval = FALSE, fig.width= 12, results = 'asis', message=FALSE} models_10_data_set <- readRDS("./data/statistical_models/combined_models_results_10_data_set_no_iucn.RDS") for(i in 1:length(models_10_data_set)){ print(names(models_10_data_set)[i]) print(models_10_data_set[[i]]$table) print(models_10_data_set[[i]]$gg) } ``` ## Model comparison ```{r calculate model comparison 10 data set, eval = FALSE} brms_model_rds <- list.files("./data/statistical_models", pattern = "10-data", full.names = TRUE) mean_looic <- pblapply(brms_model_rds, cl = 1, function(e) { mods <- readRDS(e) loos <- sapply(mods, function(x) x$criteria$loo$estimates)[c(3, 6), ] out_df <- data.frame( models = gsub("-10-data-sets_brm_multiple.RDS", "", basename(e)), looic = mean(loos[1 , ]), min.looic = min(loos[1 , ]), max.looic = max(loos[1 , ]), SE = mean(loos[2 , ]) ) return(out_df) }) mean_looic <- do.call(rbind, mean_looic) saveRDS(mean_looic, "mean_looic_by_model.RDS") ``` ```{r model comparison 10 data set, eval = FALSE, fig.width= 12, results = 'asis', message=FALSE} mean_looic <- readRDS("mean_looic_by_model.RDS") for(i in c("FDis", "FDiv", "FEve", "taxdis")){ sub_looic <- mean_looic[grep(i, mean_looic$models), ] sub_looic <- sub_looic[order(sub_looic$looic, decreasing = FALSE), ] sub_looic$delta.looic <- sub_looic$looic - min(sub_looic$looic) df1 <- knitr::kable(sub_looic, row.names = TRUE, escape = FALSE, format = "html", digits = 3, caption = i) df1 <- kable_styling(df1, bootstrap_options = c("striped", "hover", "condensed", "responsive"), full_width = FALSE, font_size = 12) print(df1) } ``` # Results for 100 data set models (propagating FD uncertainty) ```{r plots and tables 100 data set models, eval = FALSE, fig.width= 12, results = 'asis', message=FALSE} models_100_data_set <- readRDS("./data/statistical_models/combined_models_results_10_data_sets.RDS") for(i in 1:length(models_100_data_set)){ print(names(models_100_data_set)[i]) print(models_100_data_set[[i]]$table) print(models_100_data_set[[i]]$gg) } ``` # Print same model with different data set numbers together ```{r print output of models with different data set numbers together, eval = FALSE, fig.width= 12, results = 'asis'} all_models <- c(models_100_data_set, models_10_data_set, models_1_data_set) for(i in sort(names(all_models))){ print(names(all_models)[i]) print(all_models[[i]]$table) print(all_models[[i]]$gg) } ``` ```{r testing of statistical analysis BAYESIAN brms accounting for spatial autocorrelation, eval = FALSE, include = FALSE} random_FD_list <- readRDS("./data/100_replicated_random_assemblages_data_set_400m2_target_area.RDS") # loop to run bayesian mcmc glmmm in brms (using stan) controlling for spatial autocorrelation brms_model_list <- lapply(random_FD_list[1:2], function(Y){ # scale parameters vars <- c("sst_range", "chl_range", "gravity", "years_protected") for(i in vars) Y[, i] <- scale(Y[, i]) Y$tax_distinctnes <- Y$tax_distinctnes / (max(Y$tax_distinctnes) * 1.0001) # get distance matrix dst <- dist(Y[, c("lat", "long")]) distance_matrix <- as.matrix(dst) # convert to "correlation" matrix distance_matrix <- distance_matrix / max(distance_matrix) corr_matrix <- 1 - distance_matrix # add site_id labels colnames(corr_matrix) <- rownames(corr_matrix) <- Y$site_id # run model eval_models <- warbleR:::pblapply_wrblr_int(models, function(y){ y <- as.character(y) y <- as.formula(paste(y[2], y[1], y[3], "+ (1|gr(site_id, cov = A))")) print(y) brm(formula = y, iter = 10000, thin = 10, data = Y, family = Beta(), control = list(adapt_delta = 0.95), silent = 2, data2 = list(A = corr_matrix), chains = 2 #, # prior = c( # prior(normal(0, 10), "b"), # prior(normal(0, 50), "Intercept"), # prior(student_t(3, 0, 20), "sd"), # prior(student_t(3, 0, 20), "sigma") # ) ) } ) return(eval_models) }) saveRDS(brms_model_list, "./data/DELETE_test_models_from_2_random_assemblages_accounting_for_spatial_autocorrelation.RDS") # 4 h 29 min in 1 iteration ``` ```{r testing of statistical analysis BAYESIAN brms, eval = FALSE, include = FALSE} comb_models <- combine_models(mlist = out, check_data = FALSE) waic(comb_models) out2 <- lapply(out, add_criterion, "waic") loo_compare(out2[[1]], out2[[2]], criterion = "waic") brms::loo_compare(x = out2, criterion = "waic") print(out2[[1]]$criteria$) print(out2[[2]]$criteria$loo) # comb_models <- add_criterion(comb_models, criterion = "waic") summary(comb_models) ``` # Takeaways {.unnumbered .unlisted} <!-- '---' adds a gray vertical line --> <!-- add packages used, system details and versions --> # Session information {.unnumbered .unlisted} ```{r session info, echo=F} sessionInfo() ``` 
