Update:

- recalculate all analyses using R-Ratepol 0.6.0
- change the 'early' and 'late' to follow the paleo nomenclature
- add 3 new sensitivity analyses
- recreated all data and all figures

Author: OndrejMottl

R/00_config.R

@@ -63,18 +63,18 @@ high_diversity <-  50
 
 # position of the enviromental change in the sequence 
 breaks_late <-  c(2000, 3000)
-breaks_recent <-  c(5500, 6500)
+breaks_early <-  c(5500, 6500)
 
 # Number of simulated datasest of pollen data
 N_rep <-  100
 
 # template of time sequence with uneven distribution of points
-time_seq <-  data_example$list_ages[[4]]$ages$age
+time_seq <-  data_example$sample_age[[4]]$age
 
 
 # 3.2. ROC Rratepol setting -----
 age_lim <- 8e3
-roc_n_rand <- 10e3
+roc_n_rand <- 1e3
 pollen_grains <- 150
 
 

R/01_data_creation.R

@@ -17,15 +17,15 @@ source("R/00_config.R")
 # 1. Simulate datasets -----
 #----------------------------------------------------------#
 
-# low diversity recent
+# low diversity early
 
-sim_ld_recent <- 
+sim_ld_early <- 
   .simulate.pollen.data.in.multiple.datasets(
     time = time_seq, 
     nforc = N_env, 
     nprox = low_diversity, 
     manual_edit = TRUE,
-    breaks = breaks_recent,
+    breaks = breaks_early,
     jitter = TRUE,
     rarity = TRUE,
     N_datasets = N_rep)
@@ -43,13 +43,13 @@ sim_ld_late <-
     N_datasets=N_rep)
 
 
-sim_hd_recent <- 
+sim_hd_early <- 
   .simulate.pollen.data.in.multiple.datasets(
     time=time_seq, 
     nforc=N_env, 
     nprox=high_diversity, 
     manual_edit = T,
-    breaks=breaks_recent,
+    breaks=breaks_early,
     jitter = T,
     rarity=T,
     N_datasets=N_rep)
@@ -75,19 +75,19 @@ simulated_dataset <-
   dpyr::bind_rows(
 
     tibble::tibble(
-      sim_ld_recent,
+      sim_ld_early,
       diversity = "low_diversity",
-      position = "breaks_recent"),
+      position = "breaks_early"),
 
     tibble::tibble(
       sim_ld_late,
       diversity = "low_diversity",
       position = "breaks_late"),
 
     tibble::tibble(
-      sim_ld_recent,
+      sim_ld_early,
       diversity = "high_diversity",
-      position = "breaks_recent"),
+      position = "breaks_early"),
 
     tibble::tibble(
       sim_ld_late,
@@ -108,4 +108,4 @@ simulated_dataset
 #----------------------------------------------------------#
 write_rds(
   simulated_dataset,
-  "data/output/datasets/simulated/simulated_dataset.rds")
+  "data/output/datasets/simulated/simulated_dataset.rds")

R/02_ROC_of_simulated_data.R

@@ -27,7 +27,7 @@ simulated_dataset <-
 
 ROC_levels <- 
   .estimate.RoC.by.all.methods(
-    simulated_dataset,
+    random_data = simulated_dataset,
     Working_Unit = "levels", 
     interest_threshold = 8000)
 

R/03_model_building.R

@@ -31,9 +31,9 @@ data_sum <-
 
     position = factor(
       position,
-      levels = c("breaks_recent", "breaks_late")) %>%
+      levels = c("breaks_early", "breaks_late")) %>%
       fct_recode(
-        "low density level" = "breaks_recent",
+        "low density level" = "breaks_early",
         "high density level" = "breaks_late"),
     diversity = factor(
       diversity,
@@ -80,7 +80,6 @@ data_sum <-
     RoC_setting = paste0(smooth," - ",DC) %>% 
       as.factor()) 
 
-
 # explore data
 data_sum %>% 
   summary()
@@ -326,28 +325,38 @@ mod_detail_correct_dd %>%
   View()
 
 mod_detail_correct_m1 <-  
-  glmmTMB(success ~ DC + position +  smooth + 
-            DC:position + DC:smooth + position:smooth + 
-            DC:position:smooth +  # fixed effets
+  glmmTMB(success ~ position + smooth +
+            position:smooth + # fixed effets
             (1|dataset_ID), # random effects
           data = data_detail_correct,
           family = beta_family(link = "logit")
   )
 
 
 mod_detail_correct_m2 <-  
-  glmmTMB(success ~ DC + diversity + position +  smooth + 
-            DC:position + DC:smooth + position:smooth + 
-            DC:position:smooth + # fixed effets
+  glmmTMB(success ~ DC + position + smooth +
+            position:smooth + # fixed effets
+            (1|dataset_ID), # random effects
+          data = data_detail_correct,
+          family = beta_family(link = "logit")
+  )
+
+mod_detail_correct_m3 <-  
+  glmmTMB(success ~ diversity + position + smooth +
+            position:smooth + # fixed effets
             (1|dataset_ID), # random effects
           data = data_detail_correct,
           family = beta_family(link = "logit")
   )
 
+
+
+
 mod_detail_correct_final_comp <-
   compare_performance(
     mod_detail_correct_m1,
     mod_detail_correct_m2,
+    mod_detail_correct_m3,
     rank = T)
 
 mod_detail_correct_final_comp
@@ -356,16 +365,17 @@ mod_detail_correct_final_comp %>%
   as_tibble() %>% 
   left_join(.,
             tibble(
-              Model = paste0("mod_detail_correct_m",1:2),
-              formula = lapply(paste0("mod_detail_correct_m",1:2), FUN = function(x) formula(get(x))) %>% 
+              Model = paste0("mod_detail_correct_m",1:nrow(mod_detail_correct_final_comp)),
+              formula = lapply(paste0("mod_detail_correct_m",1:nrow(mod_detail_correct_final_comp)),
+                               FUN = function(x) formula(get(x))) %>% 
                 as.character()),
             by = "Model") %>% 
   dplyr::select(formula, everything()) %>% 
   write_csv(
     .,"data/output/result_tables/mod_detail_correct_final_comp.csv")
 
-# -> model 1 is the best
-mod_detail_correct_select <-  mod_detail_correct_m1
+# -> model 2 is the best
+mod_detail_correct_select <-  mod_detail_correct_m2
 write_rds(mod_detail_correct_select,"data/output/models/mod_detail_correct_select.rds")
 
 #--------------------------------#
@@ -459,89 +469,33 @@ mod_detail_false_dd %>%
 
 mod_detail_false_m1 <- 
   glmmTMB(success ~ DC + position +  smooth + 
-            DC:smooth + position:smooth + 
+            position:smooth +
             (1|dataset_ID), # random effects
           data = data_detail_false,
           family = beta_family(link = "logit"))
 
 mod_detail_false_m2 <- 
   glmmTMB(success ~ DC + position +  smooth + 
-            DC:position + DC:smooth + position:smooth + 
+            DC:position + position:smooth + 
             DC:position:smooth + 
             (1|dataset_ID), # random effects
           data = data_detail_false,
           family = beta_family(link = "logit"))
 
 mod_detail_false_m3 <- 
   glmmTMB(success ~ DC + diversity + position +  smooth + 
-            DC:diversity + DC:smooth + position:smooth + 
-            (1|dataset_ID), # random effects
-          data = data_detail_false,
-          family = beta_family(link = "logit"))
-
-mod_detail_false_m4 <- 
-  glmmTMB(success ~ DC + diversity + position +  smooth + 
-            DC:diversity + DC:position + DC:smooth + position:smooth + 
-            DC:position:smooth +
-            (1|dataset_ID), # random effects
-          data = data_detail_false,
-          family = beta_family(link = "logit"))
-
-mod_detail_false_m5 <- 
-  glmmTMB(success ~ DC + position +  smooth + 
-            DC:position + DC:smooth + position:smooth + 
-            (1|dataset_ID), # random effects
-          data = data_detail_false,
-          family = beta_family(link = "logit"))
-
-mod_detail_false_m6 <- 
-  glmmTMB(success ~ DC + diversity + position +  smooth + 
-            DC:diversity + DC:position + DC:smooth + position:smooth + 
-            (1|dataset_ID), # random effects
-          data = data_detail_false,
-          family = beta_family(link = "logit"))
-
-mod_detail_false_m7 <- 
-  glmmTMB(success ~ DC + diversity + position +  smooth + 
-            DC:smooth + position:smooth + 
-            (1|dataset_ID), # random effects
-          data = data_detail_false,
-          family = beta_family(link = "logit"))
-
-mod_detail_false_m8 <- 
-  glmmTMB(success ~ DC + diversity + position +  smooth + 
-            DC:diversity + DC:smooth + diversity:position + position:smooth + 
-            (1|dataset_ID), # random effects
-          data = data_detail_false,
-          family = beta_family(link = "logit"))
-
-mod_detail_false_m9 <- 
-  glmmTMB(success ~  position +  smooth + 
             position:smooth + 
             (1|dataset_ID), # random effects
           data = data_detail_false,
           family = beta_family(link = "logit"))
 
-mod_detail_false_m10 <- 
-  glmmTMB(success ~ DC + position +  smooth + 
-            position:smooth + 
-            (1|dataset_ID), # random effects
-          data = data_detail_false,
-          family = beta_family(link = "logit"))
 
 
 mod_detail_false_final_comp <-
   compare_performance(
     mod_detail_false_m1,
     mod_detail_false_m2,
     mod_detail_false_m3,
-    mod_detail_false_m4,
-    mod_detail_false_m5,
-    mod_detail_false_m6,
-    mod_detail_false_m7,
-    mod_detail_false_m8,
-    mod_detail_false_m9,
-    mod_detail_false_m10,
     rank = T)
 
 mod_detail_false_final_comp
@@ -550,8 +504,9 @@ mod_detail_false_final_comp %>%
   as_tibble() %>% 
   left_join(.,
             tibble(
-              Model = paste0("mod_detail_false_m",1:10),
-              formula = lapply(paste0("mod_detail_false_m",1:10), FUN = function(x) formula(get(x))) %>% 
+              Model = paste0("mod_detail_false_m",1:nrow(mod_detail_false_final_comp)),
+              formula = lapply(paste0("mod_detail_false_m",1:nrow(mod_detail_false_final_comp)),
+                               FUN = function(x) formula(get(x))) %>% 
                 as.character()),
             by = "Model") %>% 
   dplyr::select(formula, everything()) %>%