-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy path3_Bering_sea_indices_with_CPI
108 lines (98 loc) · 4.62 KB
/
3_Bering_sea_indices_with_CPI
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
library(devtools)
library(TMB)
#devtools::install_local("/FishStatsUtils-2.7.0/FishStatsUtils-2.7.0")
#devtools::install_local("/VAST-3.5.0/VAST-3.5.0", ref='--no-staged-install')
library(VAST)
library(dplyr)
gc()
#load Data_Geostat & Extrapolation list
working_directory <- getwd()
Species = c("Gadus chalcogrammus", "Gadus macrocephalus","Pleuronectes quadrituberculatus")[1] #c(pollock, Pcod,APlaice)
fish_name = c("pollock", "cod","plaice")[1] #c(pollock, Pcod,APlaice)
n_x= 500
Data_Geostat <- readRDS(file = paste0(working_directory,"/Data_Geostat_",Species,".rds"))
Extrapolation_List <- readRDS(file = paste0(working_directory,"/Extrapolation_List.rds"))
Extrapolation_List$projargs="+proj=longlat +ellps=WGS84"
Region="User"
strata.limits <- data.frame('STRATA'="All_areas")
FieldConfig = matrix( c("IID","IID","IID","IID","IID","IID"), ncol=2, nrow=3, dimnames=list(c("Omega","Epsilon","Beta"),c("Component_1","Component_2")) )
RhoConfig = c("Beta1" = 0, "Beta2" = 0, "Epsilon1" = 4, "Epsilon2" = 4)
if( Species=="Gadus chalcogrammus" ) ObsModel = c(2,1) #ObsModel = c(2,0)
if( Species=="Gadus macrocephalus" ) ObsModel = c(2,1)#c(2,3)
if( Species=="Pleuronectes quadrituberculatus" ) ObsModel = c(2,1) #ObsModel = c(2,0)
#create catchability matrix
#fixlambda <- -1
fixlambda <- 0
##annual
if(fixlambda == -1){
message('turning on annual catchability')
## Annual coefficient, create model matrix with sum to zero contrasts
## to avoid confounding with betas
yearf <- factor(Data_Geostat$Year)
Q_ik <- model.matrix(~yearf, contrasts=list(yearf='contr.sum'))
## Zero out non-WBS rows so they are unaffected by the lambdas
Q_ik[which(Data_Geostat$Center!='TINRO'),] <- 0
} else {
## Constant over time
Q_ik <- matrix(ifelse(Data_Geostat$Center =='TINRO', 1, 0), ncol=1)
}
settings = make_settings( Version = "VAST_v9_2_0",
FieldConfig = FieldConfig,
RhoConfig = RhoConfig,
n_x= 500,
Region= Region,
purpose="index2",
strata.limits=data.frame(STRATA = c("All_areas","EBS","NBS","WBS")),
#zone = Extrapolation_List$zone,
bias.correct = TRUE,
use_anisotropy = TRUE,
fine_scale = TRUE,
max_cells = 5000,
knot_method = "grid")
source(paste0(getwd(),'/fit_model.R'))
fit = fit_model( "settings" = settings,
"Lat_i" = Data_Geostat[,'Lat'],
"Lon_i" = Data_Geostat[,'Lon'],
"t_i" = Data_Geostat[,'Year'],
"c_i" = rep(0,nrow(Data_Geostat)),
"b_i" = Data_Geostat[,'Catch_KG'],
"a_i" = Data_Geostat[,'AreaSwept_km2'],
"v_i" = Data_Geostat[,'Vessel']-1,
"Q_ik" = Q_ik,
"extrapolation_list" = Extrapolation_List,
"run_model" = FALSE,
"optimize_args" =list("lower"=-Inf,"upper"=Inf),
"newtonsteps" = 0,
working_dir = paste0(getwd(), "/",fish_name,"_CPE_results/"))
#load & center cold pool index
Covariate = c("None", "Bottom", "BotVar", "Cold_pool")[4]
if( Covariate=="Cold_pool" ){
CP = read.csv( paste0(getwd(),"/data/cpa_areas2018.csv") )
CP = CP[ which(CP[,'YEAR'] %in% unique(Data_Geostat[,'Year'])), 'AREA_SUM_LTE2' ]
f = function(a) a # Works better using natural NOT log link
CP = ( f(CP) - mean(f(CP)) ) / sd(f(CP))
X_xtp = rep(1,n_x) %o% CP %o% c(1)
Xconfig_zcp = c(2,2) %o% 1 %o% 1
# fine-scale inputs
X_itp = rep(1,fit$spatial_list$n_i) %o% CP %o% c(1)
X_gtp = rep(1,fit$spatial_list$n_g) %o% CP %o% c(1)
}
gc()
fit = fit_model( "settings" = settings,
"Lat_i" = Data_Geostat[,'Lat'],
"Lon_i" = Data_Geostat[,'Lon'],
"t_i" = Data_Geostat[,'Year'],
"c_i" = rep(0,nrow(Data_Geostat)),
"b_i" = Data_Geostat[,'Catch_KG'],
"a_i" = Data_Geostat[,'AreaSwept_km2'],
"v_i" = Data_Geostat[,'Vessel']-1,
"Q_ik" = Q_ik,
"X_itp"=X_itp,
"X_gtp"=X_gtp,
"Xconfig_zcp"=Xconfig_zcp,
"extrapolation_list" = Extrapolation_List,
"optimize_args" =list("lower"=-Inf,"upper"=Inf),
"newtonsteps" = 0,
working_dir = paste0(getwd(), "/",fish_name,"_CPE_results/"))
# save the VAST model
saveRDS(fit,file = paste0(getwd(),"/",fish_name,"_CPE_results/",Species,"VASTfit_oceanography.RDS"))