-
Notifications
You must be signed in to change notification settings - Fork 0
/
02_examine_fits.R
120 lines (88 loc) · 3.45 KB
/
02_examine_fits.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
library(tidyverse)
library(spOccupancy)
#library(coda)
library(bayesplot)
library(ggpubr)
rm(list=ls())
select <- dplyr::select
theme_set(theme_classic())
#-------------
# Examine fits
#-------------
# MODEL 3: 2022_04_11 --------------------------------------------------
#occ.ms.formula <- ~ .
#det.ms.formula <- ~ .
#species_sub1 <- c("Tasmanian Devil", "cat", "Bennett's Wallaby",
# "Tasmanian Pademelon", "Spotted-tail Quoll")
# directory stuff
run_date <- "2022_04_11"
run_id <- ""
save_dir <- paste0("results/m_",run_date,"_",run_id) # git ignored
plot_dir <- paste0("plots/model_runs/",run_date) # git visible
if(!dir.exists(plot_dir)) dir.create(plot_dir)
# load fit
fit <- readRDS(paste0(save_dir,"/fit_",run_date,".rds"))
summary(fit)
# select model: alpha (detection) or beta (occupancy) model
var_type <- "alpha"
c("Devil","Cat", "Wallaby","Pademelon","Quoll") %>%
map(~ fit[[paste0(var_type,".samples")]] %>%
as.data.frame %>%
select(contains(.x)) %>%
rename_with(~ .x %>% str_remove("\\-.*")) %>%
mcmc_areas() +
xlim(-20,20) +
labs(title = str_remove(.x,".* "))) %>%
ggarrange(plotlist = .) %>%
annotate_figure(top = paste(if_else(var_type=="alpha", "Detection","Occupancy"), "model") %>%
text_grob(size = 14, face = "bold"))
#ggsave(paste0(plot_dir,"/est_",if_else(var_type=="alpha", "Detection","Occupancy"),".jpg"))
unmarked::occuMulti()
unmarked::
# MODEL 2: 2022_04_07 --------------------------------------------------
#occ.ms.formula <- ~ . - site
#det.ms.formula <- ~ . - site
#species_sub1 <- c("Tasmanian Devil", "cat", "Bennett's Wallaby",
# "Tasmanian Pademelon", "Spotted-tail Quoll")
run_date <- "2022_04_07"
run_id <- ""
save_dir <- paste0("results/m_",run_date,"_",run_id);save_dir
fit <- readRDS(paste0(save_dir,"/fit_",run_date,".rds"))
# community-level estimates
fit$beta.comm.samples %>% mcmc_areas # occupancy
fit$alpha.comm.samples %>% mcmc_areas() # detection
# species-level estimates
fit$alpha.samples %>% mcmc_areas
fit$beta.samples %>% colnames
c("Devil","Wallaby","Pademelon","Quoll") %>%
walk(~ {fit$beta.samples %>% as.data.frame %>%
select(contains(.x)) %>%
mcmc_areas() + labs(title = .x)} %>%
ggsave(paste0(save_dir,"/beta_",.x,".pdf"), plot = .))
# MODEL 1: 2022_03_28 --------------------------------------------------
#occ.ms.formula <- ~ . - site
#det.ms.formula <- ~ . - site
#species_sub1 <- c("Tasmanian Devil", "Bennett's Wallaby",
# "Tasmanian Pademelon", "Spotted-tail Quoll")
run_date <- "2022_03_28"
run_id <- ""
save_dir <- paste0("results/m_",run_date,"_",run_id);save_dir
fit <- readRDS(paste0(save_dir,"/fit_",run_date,".rds"))
# community-level estimates
fit$beta.comm.samples %>% mcmc_areas # occupancy
fit$alpha.comm.samples %>% mcmc_areas() # detection
# species-level estimates
fit$alpha.samples %>% mcmc_areas
fit$beta.samples %>% colnames
fit$beta.samples %>% as.data.frame %>% colnames()
rename_with(~ .x %>% str_remove("-Tasmanian Devil"))
plots <- c("Devil","Cat", "Wallaby","Pademelon","Quoll") %>%
map(~ fit$beta.samples %>% as.data.frame %>%
select(contains(.x)) %>%
rename_with(~ .x %>% str_remove("\\-.*")) %>%
mcmc_areas() +
xlim(-20,20) +
labs(title = str_remove(.x,".* ")))
#ggsave(paste0(save_dir,"/beta_",.x,".pdf"), plot = .))
ggpubr::ggarrange(plotlist = plots)
fit$X.p