-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathFit_example_species_SLOPE.R
129 lines (86 loc) · 3.58 KB
/
Fit_example_species_SLOPE.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
121
122
123
124
125
126
## building a BYM route-level trend model for the BBS
setwd("C:/GitHub/BBS_iCAR_route_trends")
library(tidyverse)
library(bbsBayes)
# library(cmdstanr)
# # library(rstan)
# # rstan_options(auto_write = TRUE, javascript = FALSE)
# # library(shinystan)
# library(sf)
# library(spdep)
# # library(doParallel)
# # library(foreach)
# library(ggforce)
# #library(tidybayes)
# #source("functions/mungeCARdata4stan.R")
# source("functions/neighbours_define.R") ## function to define neighbourhood relationships
# source("functions/prepare-jags-data-alt.R") ## small alteration of the bbsBayes function
# source("functions/get_basemap_function.R") ## loads one of the bbsBayes strata maps
# source("functions/posterior_summary_functions.R") ## functions similar to tidybayes that work on cmdstanr output
# ## changes captured in a commit on Nov 20, 2020
# load and stratify CASW data ---------------------------------------------
#species = "Pacific Wren"
#species = "Barred Owl"
strat = "bbs_usgs"
model = "slope"
firstYear = 2004
lastYear = 2019
scope = "RangeWide"
species = "Blue-headed Vireo"
#species = "Dickcissel"
species_f <- gsub(species,pattern = " ",replacement = "_",fixed = T)
spp <- "SLOPE"
# out_base <- paste0(species_f,spp,firstYear,"_",lastYear)
#
#
# # SPECIES LOOP ------------------------------------------------------------
#
# output_dir <- "output"
#
#
# sp_file <- paste0(output_dir,"/",species_f,"_",scope,"_",firstYear,"_",lastYear,"_slope_route_iCAR.RData")
#
# load(sp_file)
# output_dir <- "output"
# out_base <- paste0(species_f,spp,firstYear,"_",lastYear)
#
# stan_data$ncounts
# jags_data$ncounts
strat_data = stratify(by = strat)
jags_data = bbsBayes::prepare_data(strat_data = strat_data,
species_to_run = species,
model = model,
#n_knots = 10,
min_year = firstYear,
max_year = lastYear,
min_n_routes = 1)# spatial neighbourhood define --------------------------------------------
jags_mod = run_model(jags_data = jags_data,
parameters_to_save = c("nslope"),
model_file_path = "temp.txt",
parallel = TRUE,
n_chains = 3)
inds_n <- generate_indices(jags_mod = jags_mod,
jags_data = jags_data)
inds_slope <- generate_indices(jags_mod = jags_mod,
jags_data = jags_data,
alternate_n = "nslope")
t_slope <- generate_trends(inds_slope)
n_slope <- generate_trends(inds_n)
write.csv(t_slope,file = "data/bbsBayes_slope_trends_example.csv",
row.names = FALSE)
write.csv(inds_slope$data_summary,file = "data/bbsBayes_slope_indices_example.csv",
row.names = FALSE)
# model with no year-effects ----------------------------------------------
jags_mod = run_model(jags_data = jags_data,
parameters_to_save = c("nslope"),
model_file_path = "temp_slope.txt",
parallel = TRUE,
n_chains = 3)
inds_slope <- generate_indices(jags_mod = jags_mod,
jags_data = jags_data,
alternate_n = "nslope")
t_slope <- generate_trends(inds_slope)
write.csv(t_slope,file = "data/bbsBayes_slope_only_trends_example.csv",
row.names = FALSE)
write.csv(inds_slope$data_summary,file = "data/bbsBayes_slope_only_indices_example.csv",
row.names = FALSE)