sandhill_crane

# package names
packages<-c("tidyverse", "here", "mcp", "lubridate", "knitr", "ezknitr", "loo", "flextable")

# install any packages not previously installed
installed_packages<-packages %in% rownames(installed.packages())
if(any(installed_packages == FALSE)){
  install.packages(packages[!installed_packages])
}

# load packages
invisible(lapply(packages, library, character.only = TRUE))

Import data

df<-read_csv(here("DRUM_materials/data_files_for_DRUM/sandhill_crane.csv"))
head(df)

## # A tibble: 6 x 4
##   t_                  days_numeric nsd_daily_mean displacement
##   <dttm>                     <dbl>          <dbl>        <dbl>
## 1 2015-09-27 12:53:45           76    21150616787        145. 
## 2 2015-09-28 00:07:38           76    15976347654        126. 
## 3 2015-09-29 00:05:46           77     1127371670         33.6
## 4 2015-09-30 00:03:45           78     1142404546         33.8
## 5 2015-10-01 00:02:14           79     1130261578         33.6
## 6 2015-10-02 00:00:45           80     1140330834         33.8

I’ll use net-squared displacement back-transformed to simple displacement for easier interpretation.

m1<-list(displacement~1)
m2<-list(displacement~1,~1)
m3<-list(displacement~1,~1,~1)
m4<-list(displacement~1,~1,~1,~1)
m5<-list(displacement~1,~1,~1,~1,~1)
m6<-list(displacement~1,~1,~1,~1,~1,~1)
m7<-list(displacement~1,~1,~1,~1,~1,~1,~1)

# Fit models
int_mods<-list(m1, m2, m3, m4, m5, m6, m7)
out_mods<-list()
  
for(i in 1:length(int_mods)){
  out_mods[[i]]<-mcp(model = int_mods[[i]], data = df[,c("days_numeric", "displacement")],
      par_x = "days_numeric", adapt=5000)
      }

for(i in 1:length(out_mods)){
  out_mods[[i]]$loo<-loo(out_mods[[i]])
}

loo_compare(out_mods[[1]]$loo,
            out_mods[[2]]$loo,
            out_mods[[3]]$loo,
            out_mods[[4]]$loo,
            out_mods[[5]]$loo,
            out_mods[[6]]$loo,
            out_mods[[7]]$loo)

Due to the random nature of MCMC sampling, some run output may differ slightly. The following is saved model output.

fit1_int<-readRDS(here("other_datasets/saved_model_output/sacr_multiple_intercepts/sacr_1_intercept.rds"))
fit2_int<-readRDS(here("other_datasets/saved_model_output/sacr_multiple_intercepts/sacr_2_intercepts.rds"))
fit3_int<-readRDS(here("other_datasets/saved_model_output/sacr_multiple_intercepts/sacr_3_intercepts.rds"))
fit4_int<-readRDS(here("other_datasets/saved_model_output/sacr_multiple_intercepts/sacr_4_intercepts.rds"))
fit5_int<-readRDS(here("other_datasets/saved_model_output/sacr_multiple_intercepts/sacr_5_intercepts.rds"))
fit6_int<-readRDS(here("other_datasets/saved_model_output/sacr_multiple_intercepts/sacr_6_intercepts.rds"))
fit7_int<-readRDS(here("other_datasets/saved_model_output/sacr_multiple_intercepts/sacr_7_intercepts.rds"))

loos<-loo_compare(fit1_int$loo,
            fit2_int$loo,
            fit3_int$loo,
            fit4_int$loo,
            fit5_int$loo,
            fit6_int$loo,
            fit7_int$loo)

loo_df<-as.data.frame(loos[1:7, 1:2])
loo_df[,1:2]<-lapply(loo_df[,1:2], round, 2)
loo_df<-cbind.data.frame('Model Syntax'=c("Six Intercepts", 
                                          "Five Intercepts",
                                          "Seven Intercepts",
                                          "Four Intercepts",
                                          "Two Intercepts",
                                          "Three Intercepts",
                                          "One Intercept"), loo_df)

loo_df %>%  flextable() %>% 
  align(part = "all") %>% # left align
  set_header_labels(values=list("Model Syntax" = "Model Syntax", 
                       "elpd_diff" = "ELPD Difference", 
                      "se_diff" = "SE of Difference")) %>% 
  set_caption(caption = "Leave-one-out Cross Validation was used to compute the Estimated Log Predictive Density (ELPD) for 7 models, each with an increasing number of intercepts. The highest ELPD for the model with 6 intercepts indicates that it has the highest predictive accuracy.") %>% 
  # font(fontname = "Calibri (Body)", part = "all") %>%
  fontsize(size = 12, part = "body") %>% 
  # add footer if you want
  # add_footer_row(values = "* p < 0.05. ** p < 0.01. *** p < 0.001.", 
  #                colwidths = 4) %>% 
  theme_booktabs() %>% # default theme
  border_outer() %>% 
  vline(part="header") %>% 
  align(align ='center', part = 'all') %>% 
  bg(part="header", bg="#dbdbdb") %>% 
  bold(part="header") %>% 
  colformat_num(digits = 2) %>% 
  autofit()

Leave-one-out Cross Validation was used to compute the Estimated Log Predictive Density (ELPD) for 7 models, each with an increasing number of intercepts. The highest ELPD for the model with 6 intercepts indicates that it has the highest predictive accuracy.
Model Syntax	ELPD Difference	SE of Difference
Six Intercepts	0.00	0.00
Five Intercepts	-0.48	0.83
Seven Intercepts	-19.59	32.82
Four Intercepts	-24.73	9.52
Two Intercepts	-57.39	14.01
Three Intercepts	-324.69	16.93
One Intercept	-336.17	16.73

Plot the best model.

crane<-readRDS(here("other_datasets/saved_model_output/sacr_migration_staging.rds"))
plot(crane, q_predict=T, q_fit=T)+
  xlab("Julian Date")+ylab("Displacement from breeding territory (in km)\n")+
  theme_bw()+
  theme(plot.title=element_text(hjust=0.5))

Footer

sessionInfo()

## R version 4.0.2 (2020-06-22)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 19043)
## 
## Matrix products: default
## 
## locale:
## [1] LC_COLLATE=English_United States.1252 
## [2] LC_CTYPE=English_United States.1252   
## [3] LC_MONETARY=English_United States.1252
## [4] LC_NUMERIC=C                          
## [5] LC_TIME=English_United States.1252    
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] flextable_0.6.8 loo_2.4.1       ezknitr_0.6     knitr_1.30     
##  [5] lubridate_1.7.9 mcp_0.3.0.9000  here_0.1        forcats_0.5.0  
##  [9] stringr_1.4.0   dplyr_1.0.2     purrr_0.3.4     readr_1.3.1    
## [13] tidyr_1.1.2     tibble_3.0.3    ggplot2_3.3.2   tidyverse_1.3.0
## 
## loaded via a namespace (and not attached):
##  [1] Rcpp_1.0.8         assertthat_0.2.1   rprojroot_1.3-2    digest_0.6.28     
##  [5] R6_2.4.1           cellranger_1.1.0   backports_1.1.10   reprex_0.3.0      
##  [9] evaluate_0.14      httr_1.4.2         pillar_1.4.6       gdtools_0.2.2     
## [13] rlang_0.4.10       uuid_0.1-4         readxl_1.3.1       data.table_1.13.0 
## [17] rstudioapi_0.11    jquerylib_0.1.4    blob_1.2.1         checkmate_2.0.0   
## [21] rmarkdown_2.11     munsell_0.5.0      broom_0.7.0        compiler_4.0.2    
## [25] modelr_0.1.8       xfun_0.26          systemfonts_0.2.3  base64enc_0.1-3   
## [29] pkgconfig_2.0.3    htmltools_0.5.2    tidyselect_1.1.0   matrixStats_0.57.0
## [33] fansi_0.4.1        crayon_1.3.4       dbplyr_1.4.4       withr_2.3.0       
## [37] grid_4.0.2         jsonlite_1.7.1     gtable_0.3.0       lifecycle_0.2.0   
## [41] DBI_1.1.1          magrittr_2.0.1     scales_1.1.1       zip_2.1.1         
## [45] cli_2.0.2          stringi_1.5.3      fs_1.5.0           xml2_1.3.2        
## [49] bslib_0.3.1        ellipsis_0.3.1     generics_0.0.2     vctrs_0.3.4       
## [53] tools_4.0.2        glue_1.4.2         officer_0.4.0      hms_0.5.3         
## [57] parallel_4.0.2     fastmap_1.1.0      yaml_2.2.1         colorspace_1.4-1  
## [61] rvest_0.3.6        haven_2.3.1        patchwork_1.1.1    sass_0.4.0

sandhill_crane

David

3/11/2022