-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathfit_model_X_GTP.R
150 lines (146 loc) · 9.36 KB
/
fit_model_X_GTP.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
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
fit_model <- function (settings, Lat_i, Lon_i, t_i, b_i, a_i,
c_iz = rep(0, length(b_i)), v_i = rep(0, length(b_i)),
working_dir = paste0(getwd(), "/"), X1config_cp = NULL,
X2config_cp = NULL, covariate_data, X1_formula = ~0,
X2_formula = ~0, Q1config_k = NULL, Q2config_k = NULL,
catchability_data, Q1_formula = ~0, Q2_formula = ~0,
newtonsteps = 1, silent = TRUE, build_model = TRUE,
run_model = TRUE, test_fit = TRUE,
#####################
##### Zack's Addition
X_gtp = NULL,
#####################
...)
{
extra_args = list(...)
extra_args = c(extra_args, extra_args$extrapolation_args,
extra_args$spatial_args, extra_args$optimize_args, extra_args$model_args)
data_frame = data.frame(Lat_i = Lat_i, Lon_i = Lon_i, a_i = a_i,
v_i = v_i, b_i = b_i, t_i = t_i, c_iz = c_iz)
year_labels = seq(min(t_i), max(t_i))
years_to_plot = which(year_labels %in% t_i)
message("\n### Writing output from `fit_model` in directory: ",
working_dir)
dir.create(working_dir, showWarnings = FALSE, recursive = TRUE)
capture.output(settings, file = file.path(working_dir, "settings.txt"))
message("\n### Making extrapolation-grid")
extrapolation_args_default = list(Region = settings$Region,
strata.limits = settings$strata.limits, zone = settings$zone,
max_cells = settings$max_cells)
extrapolation_args_input = combine_lists(input = extra_args,
default = extrapolation_args_default, args_to_use = formalArgs(make_extrapolation_info))
extrapolation_list = do.call(what = make_extrapolation_info,
args = extrapolation_args_input)
message("\n### Making spatial information")
spatial_args_default = list(grid_size_km = settings$grid_size_km,
n_x = settings$n_x, Method = settings$Method, Lon_i = Lon_i,
Lat_i = Lat_i, Extrapolation_List = extrapolation_list,
DirPath = working_dir, Save_Results = TRUE, fine_scale = settings$fine_scale,
knot_method = settings$knot_method)
spatial_args_input = combine_lists(input = extra_args, default = spatial_args_default,
args_to_use = c(formalArgs(make_spatial_info), formalArgs(INLA::inla.mesh.create)))
spatial_list = do.call(what = make_spatial_info, args = spatial_args_input)
message("\n### Making data object")
if (missing(covariate_data))
covariate_data = NULL
if (missing(catchability_data))
catchability_data = NULL
data_args_default = list(Version = settings$Version, FieldConfig = settings$FieldConfig,
OverdispersionConfig = settings$OverdispersionConfig,
RhoConfig = settings$RhoConfig, VamConfig = settings$VamConfig,
ObsModel = settings$ObsModel, c_iz = c_iz, b_i = b_i,
a_i = a_i, v_i = v_i, s_i = spatial_list$knot_i - 1,
t_i = t_i, spatial_list = spatial_list, Options = settings$Options,
Aniso = settings$use_anisotropy, X1config_cp = X1config_cp,
X2config_cp = X2config_cp, covariate_data = covariate_data,
X1_formula = X1_formula, X2_formula = X2_formula, Q1config_k = Q1config_k,
Q2config_k = Q2config_k, catchability_data = catchability_data,
Q1_formula = Q1_formula, Q2_formula = Q2_formula)
data_args_input = combine_lists(input = extra_args, default = data_args_default)
data_list = do.call(what = make_data, args = data_args_input)
#############Zack Addition####################
data_list$X1_gctp = data_list$X2_gctp = X_gtp
#############Zack Addition####################
message("\n### Making TMB object")
model_args_default = list(TmbData = data_list, RunDir = working_dir,
Version = settings$Version, RhoConfig = settings$RhoConfig,
loc_x = spatial_list$loc_x, Method = spatial_list$Method,
build_model = build_model)
model_args_input = combine_lists(input = extra_args, default = model_args_default,
args_to_use = formalArgs(make_model))
tmb_list = do.call(what = make_model, args = model_args_input)
if (run_model == FALSE | build_model == FALSE) {
input_args = list(extra_args = extra_args, extrapolation_args_input = extrapolation_args_input,
model_args_input = model_args_input, spatial_args_input = spatial_args_input,
data_args_input = data_args_input)
Return = list(data_frame = data_frame, extrapolation_list = extrapolation_list,
spatial_list = spatial_list, data_list = data_list,
tmb_list = tmb_list, year_labels = year_labels, years_to_plot = years_to_plot,
settings = settings, input_args = input_args)
class(Return) = "fit_model"
return(Return)
}
if (silent == TRUE)
tmb_list$Obj$env$beSilent()
if (test_fit == TRUE) {
message("\n### Testing model at initial values")
LogLike0 = tmb_list$Obj$fn(tmb_list$Obj$par)
Gradient0 = tmb_list$Obj$gr(tmb_list$Obj$par)
if (any(Gradient0 == 0)) {
message("\n")
stop("Please check model structure; some parameter has a gradient of zero at starting values\n",
call. = FALSE)
}
else {
message("Looks good: All fixed effects have a nonzero gradient")
}
}
message("\n### Estimating parameters")
optimize_args_default1 = combine_lists(default = list(lower = tmb_list$Lower,
upper = tmb_list$Upper, loopnum = 2), input = extra_args,
args_to_use = formalArgs(TMBhelper::fit_tmb))
optimize_args_input1 = list(obj = tmb_list$Obj, savedir = NULL,
newtonsteps = 0, bias.correct = FALSE, control = list(eval.max = 10000,
iter.max = 10000, trace = 1), quiet = TRUE, getsd = FALSE)
optimize_args_input1 = combine_lists(default = optimize_args_default1,
input = optimize_args_input1, args_to_use = formalArgs(TMBhelper::fit_tmb))
parameter_estimates = do.call(what = TMBhelper::fit_tmb,
args = optimize_args_input1)
if (exists("check_fit") & test_fit == TRUE) {
problem_found = VAST::check_fit(parameter_estimates)
if (problem_found == TRUE) {
message("\n")
stop("Please change model structure to avoid problems with parameter estimates and then re-try; see details in `?check_fit`\n",
call. = FALSE)
}
}
optimize_args_default2 = list(obj = tmb_list$Obj, lower = tmb_list$Lower,
upper = tmb_list$Upper, savedir = working_dir, bias.correct = settings$bias.correct,
newtonsteps = newtonsteps, bias.correct.control = list(sd = FALSE,
split = NULL, nsplit = 1, vars_to_correct = settings$vars_to_correct),
control = list(eval.max = 10000, iter.max = 10000, trace = 1),
loopnum = 1)
optimize_args_input2 = combine_lists(input = extra_args,
default = optimize_args_default2, args_to_use = formalArgs(TMBhelper::fit_tmb))
optimize_args_input2 = combine_lists(input = list(startpar = parameter_estimates$par),
default = optimize_args_input2)
parameter_estimates = do.call(what = TMBhelper::fit_tmb,
args = optimize_args_input2)
Report = tmb_list$Obj$report()
ParHat = tmb_list$Obj$env$parList(parameter_estimates$par)
input_args = list(extra_args = extra_args, extrapolation_args_input = extrapolation_args_input,
model_args_input = model_args_input, spatial_args_input = spatial_args_input,
optimize_args_input1 = optimize_args_input1, optimize_args_input2 = optimize_args_input2,
data_args_input = data_args_input)
Return = list(data_frame = data_frame, extrapolation_list = extrapolation_list,
spatial_list = spatial_list, data_list = data_list, tmb_list = tmb_list,
parameter_estimates = parameter_estimates, Report = Report,
ParHat = ParHat, year_labels = year_labels, years_to_plot = years_to_plot,
settings = settings, input_args = input_args, X1config_cp = X1config_cp,
X2config_cp = X2config_cp, covariate_data = covariate_data,
X1_formula = X1_formula, X2_formula = X2_formula, Q1config_k = Q1config_k,
Q2config_k = Q1config_k, catchability_data = catchability_data,
Q1_formula = Q1_formula, Q2_formula = Q2_formula)
class(Return) = "fit_model"
return(Return)
}