Finalizing for zenodo

R/JSDM_convergence.R

@@ -1,5 +1,5 @@
-source(here::here("r", "utils_generic.R"))
-source(here::here("r", "utils_JSDM.R"))
+source(here::here("R", "utils_generic.R"))
+source(here::here("R", "utils_JSDM.R"))
 
 library(ape)
 library(Hmsc)

R/JSDM_format_qeq.R

@@ -1,4 +1,4 @@
-source(here::here("r", "utils_generic.R"))
+source(here::here("R", "utils_generic.R"))
 
 library(here)
 library(tidyverse)

R/JSDM_format_sort.R

@@ -1,4 +1,4 @@
-source(here::here("r", "utils_generic.R"))
+source(here::here("R", "utils_generic.R"))
 
 library(ape)
 library(Hmsc)

R/JSDM_inspect.R

@@ -1,5 +1,5 @@
-source(here::here("r", "utils_generic.R"))
-source(here::here("r", "utils_JSDM.R"))
+source(here::here("R", "utils_generic.R"))
+source(here::here("R", "utils_JSDM.R"))
 
 library(ape)
 library(Hmsc)
@@ -14,8 +14,6 @@ models <-
     "sort_mp_full_thin_1000_samples_250_chains_4.rds"
   )
 
-
-
 # Variance partitioning ---------------------------------------------------
 
 # Check how well model explains fraction of variance of each species
@@ -123,8 +121,6 @@ ptgr <- tgr %>% mutate(cil = case_when(cil == 0 ~ "None",
 ggsave(filename=here::here("figs", "JSDM_qeq_trait_gradient.svg"), plot=ptgr, device="svg",
        units="cm", height=8, width=25.8)
 
-ptgr
-
 # Combine plots
 vplot.s <- myvarplot(vpr2.s, snorcols)
 vplot.q <- myvarplot(vpr2.q, qnorcols)
@@ -140,13 +136,9 @@ fig4 <- vplot.s + bplot.s + gplot.s +
   plot_layout(nrow=2, ncol=3, widths = c(1,1,1), heights=c(1.7, 1), guides = 'collect') +
   plot_annotation(tag_levels="A")
 
-fig4
-
-ggsave(filename=here::here("figs", "fig4.svg"), plot=fig5, device="svg",
+ggsave(filename=here::here("figs", "fig4.svg"), plot=fig4, device="svg",
        units="cm", height=15, width=17.8)
 
-
-
 # Variance explained by traits --------------------------------------------
 
 # VP$R2T$Beta tells how much traits explain out of variation among species in 
@@ -170,6 +162,8 @@ r2t %>%
 
 # Phylogenetic signal -----------------------------------------------------
 
-summary(loadhmsccoda("qeq_ma_full_thin_1000_samples_250_chains_4.rds")$Rho)
+print("Phylogenetic signal sorting phase")
+print(summary(loadhmsccoda("sort_ma_full_thin_1000_samples_250_chains_4.rds")$Rho))
 
-summary(loadhmsccoda("sort_ma_full_thin_1000_samples_250_chains_4.rds")$Rho)
+print("Phylogenetic signal equilibrium phase")
+print(summary(loadhmsccoda("qeq_ma_full_thin_1000_samples_250_chains_4.rds")$Rho))

R/community_dissimilarity_plot.R

@@ -1,5 +1,5 @@
-source(here::here("r", "utils_generic.R"))
-source(here::here("r", "utils_community_composition.R"))
+source(here::here("R", "utils_generic.R"))
+source(here::here("R", "utils_community_composition.R"))
 
 # Format data -------------------------------------------------------------
 
@@ -70,9 +70,7 @@ p <- ggplot(dq1.f, aes(x=index)) +
 ggsave(here::here("figs", "figS1b.svg"), p, width=17.8, height=11.8, units="cm",
        device="svg")
 
-
-dq1.f %>%
+print(dq1.f %>%
   group_by(index) %>%
   summarize(m=mean(value),
-            sd=sd(value))
-
+            sd=sd(value)))

R/community_dissimilarity_regression.R

@@ -1,15 +1,25 @@
+library(here)
+library(withr)
+library(tidyverse)
+
+library(see)
+library(insight)
+library(bayestestR)
+library(performance)
+
 library(rstanarm)
 library(bayesplot)
 library(tidybayes)
-library(bayestestR)
-library(insight)
-library(see)
-library(performance)
+
 library(emmeans)
-library(betareg)
-library(withr)
 
-dq1.f <- read_rds( here::here("data", "community_dissimilarities.rds"))
+
+# Read data ---------------------------------------------------------------
+
+dq1.f <- readr::read_rds( here::here("data", "community_dissimilarities.rds"))
+
+
+# Beta regression ---------------------------------------------------------
 
 # Beta regression from rstanarm
 with_seed(2344523, 
@@ -19,21 +29,10 @@ with_seed(2344523,
               iter = 4000,
               cores = 4))
 
-mdraws <- as.matrix(m)
-
-
-# posterior predictive check... A bit hacky for the mgcg betar family.
-
-# m1 <- m
-# class(m1) <- c(class(m1), "betareg")
-
-# ppc_dens_overlay(y = m1$y, yrep = posterior_predict(m1, draws = 200))
-
-
+# posterior summary
 describe_posterior(m, ci = 0.95, rope_ci = 0.95,
                    test = c("p_direction", "rope"))
 
-
 # Plot the model results. Nothing in 95% ROPE
 p1 <- m %>%
   gather_draws(`(Intercept)`, indexshannon_div) %>%
@@ -48,7 +47,6 @@ p1 <- m %>%
 
 p1
 
-
 # Examine emmeans contrasts
 
 m %>%
@@ -57,18 +55,4 @@ m %>%
   median_qi() %>%
   xtable::xtable() %>%
   print() %>%
-  write_lines(here::here("tables", "table_S5.tex"))
-
-  # invlogit transform back to [0-1] interval
-  # mutate(.value=invlogit(.value),
-  #        .lower=invlogit(.lower),
-  #        .upper=invlogit(.upper)
-  # )
-
-# Pairwise contrasts (on log-odds scale)
-
-m %>%
-  emmeans( ~ index) %>%
-  contrast(method = "pairwise") %>%
-  gather_emmeans_draws() %>%
-  median_qi()
+  write_lines(here::here("tables", "table_S5.tex"))

R/competitive_LV.R

@@ -4,10 +4,7 @@ library(gauseR)
 library(withr)
 
 # Load data ---------------------------------------------------------------
-predator <- read_tsv(here("data", "formatted_predator_prey_density.tsv"), col_types="ccccdddd") %>%
-  mutate(replicate=factor(replicate, levels=c("A", "B", "C", "D"))) %>%
-  mutate(treatment=factor(treatment, levels=c("H", "HN", "HPanc", "HPevo", "HNPanc", "HNPevo")),
-         microcosmID=factor(microcosmID))
+predator <- readr::read_rds(here("data", "formatted_predator_prey_density.rds"))
 
 dat_ltv <- predator %>% filter(treatment == "HNPanc") %>%
   select(day, replicate, ciliate = ciliate_per_ml_imp, worm = worm_per_ml_imp) %>%
@@ -71,4 +68,4 @@ df_sum <- df %>% group_by(type) %>%
     sdkckn = sd(kckn, na.rm = T)
   )
 
-write_rds(df_sum, here::here("data", "competitive_LV.rds"))
+write_rds(df_sum, here::here("data", "competitive_LV.rds"))