-
Notifications
You must be signed in to change notification settings - Fork 1
/
Test_figures_tables.Rmd
252 lines (205 loc) · 9.91 KB
/
Test_figures_tables.Rmd
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
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
---
title: ''
author: ''
date: ''
output:
pdf_document:
fig_caption: yes
highlight: haddock
includes:
in_header: header.tex
keep_tex: yes
latex_engine: xelatex
template: Default_template_modified.tex
number_sections: yes
toc: yes
toc_depth: 4
html_document:
toc: yes
word_document: default
documentclass: article
fontsize: 12pt
geometry: margin=1in
csl: CJFAS.csl
bibliography: BibFile.bib
---
```{r global_options, include=FALSE}
# set global options for R code chunks: echo=FALSE (don't include source code);
# warning=FALSE (suppress R warnings); message=FALSE (suppress R messages)
# eval = TRUE is default
knitr::opts_chunk$set(echo = FALSE, warning = FALSE, message = FALSE)
```
```{r}
# Read in preamble R code - including required libraries and the SS file(s)
source('./Rcode/Preamble.R')
# Read in data/manipulations for executive summary tables and figures
# It may take time to edit this file and get it ready for your assessment
# Make small changes in this file and then try to compile the document
# Commit when you have a success!
source('./Rcode/Exec_summary_figs_tables.R')
```
USE THIS .Rmd TO TEST R CODE CHUNKS, FIGURES AND PLOTS BEFORE INSERTING INTO THE MAIN TEXT OR TO DEBUG
\begin{landscape}
```{r}
# test the model parameters table for SS3.30 compatibility
# If you use this for more than one model = change mod1 to mod2 or mod3
mod_params = mod1$parameters[-c(grep('Recr',mod1$parameters$Label),
grep('Impl',mod1$parameters$Label)),
(names(mod1$parameters) %in%
c("Num","Label","Value","Phase","Min",
"Max","Status","Parm_StDev","Gradient",
"PR_type","Prior","Pr_SD","Prior_Like"))]
# So that decimals won't show up in these columns and control digits
# mod_params$Num = as.factor(mod_params$Num)
# mod_params$Phase = as.factor(mod_params$Phase)
# mod_params$Prior = round(mod_params$Prior, digits = 3)
# mod_params$Pr_SD = round(mod_params$Pr_SD, digits = 3)
# Combine bounds into one column
mod_params$Min = paste('(', mod_params$Min, ', ', mod_params$Max, ')', sep='')
# Combine prior info to one column
#mod_params$PR_type = gsub('No_prior', 'None', mod_params$PR_type)
mod_params$PR_type = ifelse(mod_params$PR_type == 'No_prior' , 'None',
paste(mod_params$PR_type,' (', mod_params$Prior,
', ', mod_params$Pr_SD, ')', sep = ''))
# Remove the max, prior and prior sd columns
drops = c('Max', 'Prior', 'Pr_SD')
mod_params = mod_params[, !(names(mod_params) %in% drops)]
# Add column names
colnames(mod_params) = c('No.',
'Parameter',
'Value',
'Phase',
'Bounds',
'Status',
'SD',
'Gradient',
'Prior',
'Prior like')
# Model 1 model parameters
mod_params.table = xtable(mod_params,
caption=c(paste('List of parameters used in
the base model, including estimated
values and standard deviations (SD),
bounds (minimum and maximum),
estimation phase (negative values indicate
not estimated), status (indicates if
parameters are near bounds, and prior type
information (mean, SD).'
, sep='')),
label='tab:model_params')
# Add alignment
align(mod_params.table) = c('lrlrrcrclll')
# Add "continued on next page"" footnote
addtorow <- list()
addtorow$pos <- list()
addtorow$pos[[1]] <- c(0)
addtorow$command <- c(paste("\\hline \n",
"\\endhead \n",
"\\hline \n",
"\\multicolumn{3}{l}",
"{\\footnotesize Continued on next page} \n",
"\\endfoot \n",
"\\endlastfoot \n",sep=""))
# Print Model 1 parameters
print(mod_params.table,
include.rownames = FALSE,
caption.placement = 'top',
tabular.environment = 'longtable',
floating = FALSE,
add.to.row = addtorow)
```
\end{landscape}
```{r}
#Test the time-series table
# Total biomass, extract and subset
Bio_all = mod1$sprseries[, c('Year', 'Bio_all')]
Bio_allyrs = subset(Bio_all, Year > (Dat_start - 1) & Year < (Dat_end + 1))
Bio_allyrs$Bio_all = round(Bio_allyrs$Bio_all, 0)
# Spawning biomass, extract and subset, and turn into scientific notation
SpawningB = mod1$derived_quants[grep('SPB', mod1$derived_quants$LABEL), ]
SpawningB = SpawningB[c(-1, -2), ]
SpawningByrs = SpawningB[SpawningB$LABEL >= paste('SPB_', Dat_start, sep='') &
SpawningB$LABEL <= paste('SPB_', Dat_end, sep=''), ]
SpawningB_units = ''
if(mean(SpawningByrs$Value) > 1000000){
SpawningB_units <- "million"
SpawningByrs$Value <- SpawningByrs$Value/1000000
}
# Depletion, extract, rename and subset
Depl_years = as.data.frame(seq(Dat_start, Dat_end, 1))
colnames(Depl_years) = 'Year'
Depl_years$Depl = 0
Depletion = mod1$derived_quants[grep('Bratio', mod1$derived_quants$LABEL), ]
Depletionyrs = Depletion[Depletion$LABEL >= paste('Bratio_', Dat_start, sep = '') &
Depletion$LABEL <= paste('Bratio_', Dat_end, sep = ''), ]
Depletionyrs$Year = Depletionyrs$Label1 = substr(Depletionyrs$LABEL,
(nchar(Depletionyrs$LABEL) + 1) - 4,
nchar(Depletionyrs$LABEL))
# Make sure depletion is numeric and merge ...
Depletionyrs$Year = as.numeric(Depletionyrs$Year)
Depleteyrs = merge(Depl_years, Depletionyrs, all.x=T, all.y=T, by.x='Year', by.y='Year')
Depleteyrs[is.na(Depleteyrs)] <- 0
Depleteyrs$total = Depleteyrs$Depl + Depleteyrs$Value
# Recruits, extract and subset
Recruit = mod1$derived_quants[grep('Recr', mod1$derived_quants$LABEL), ]
Recruit = Recruit[c(-1, -2), ]
Recruityrs = Recruit[Recruit$LABEL >= paste('Recr_', Dat_start, sep='') &
Recruit$LABEL <= paste('Recr_', Dat_end, sep=''), ]
# Landings/total catch, extract and subset years
Landings = mod1$sprseries[ , c('Year','Dead_Catch_B')]
Landingsyrs = subset(Landings, Year > (Dat_start - 1) & Year < (Dat_end + 1))
# Relatvie exploitation rate, extract, subset and merge
Exploit = mod1$derived_quants[grep('F', mod1$derived_quants$LABEL), ]
Exploit = Exploit[c(-1, -2), ]
Exploityrs = Exploit[Exploit$LABEL >= paste('F_', Dat_start, sep = '') &
Exploit$LABEL <= paste('F_', Dat_end, sep = ''), ]
Exploityrs$Year = Exploityrs$Label1 = substr(Exploityrs$LABEL,
(nchar(Exploityrs$LABEL) + 1) - 4,
nchar(Exploityrs$LABEL))
Exploityrs$Year = as.numeric(Exploityrs$Year)
Exploited = merge(Depl_years, Exploityrs, all.x=T, all.y=T, by.x='Year', by.y='Year')
Exploited[is.na(Exploited)] <- 0
Exploited$total = Exploited$Depl + Exploited$Value
# SPR, extract and subset years
SPR = mod1$sprseries[, c('Year', 'SPR')]
SPRyrs = subset(SPR, Year > (Dat_start - 1) & Year < (Dat_end + 1))
# Bind all the columns together for the table
Timeseries = as.data.frame(cbind(seq(Dat_start, Dat_end, 1),
Bio_allyrs$Bio_all,
round(SpawningByrs$Value, 0),
round(Depleteyrs$total, 2),
Recruityrs$Value,
Landingsyrs$Dead_Catch_B,
round(Exploited$total,2),
round(SPRyrs$SPR, 2)))
# Add colulmn names
colnames(Timeseries)=c('Year',
'Total biomass (mt)',
paste0('Spawning biomass', SpawningB_units, '(mt)'),
'Depletion',
'Age-0 recruits',
'Total catch (mt)',
'Relative exploitation rate',
'SPR')
# Make year a factor so you don't have a decimal
Timeseries$Year = as.factor(Timeseries$Year)
# Remove 2015 values for last three columns since year isn't complete
Timeseries[nrow(Timeseries), c((ncol(Timeseries) - 2):ncol(Timeseries))] <- NA
```
```{r]
OFL_mod1 = mod1$derived_quants[grep('OFL',mod1$derived_quants$LABEL),]
OFL_mod1 = OFL_mod1[c(-1,-2),2]
#Turn into a dataframe and get the total
OFL = as.data.frame(OFL_mod1)
OFL$Year=seq(Project_firstyr+2,Project_lastyr,1)
OFL$Year = as.factor(OFL$Year)
OFL = OFL[,c(2,1)]
colnames(OFL) = c('Year','OFL')
# Create the table
OFL.table = xtable(OFL, caption=c('Projections of potential OFL (mt) for each model, using the base model forecast.'),
label = 'tab:OFL_projection')}
```
```{r, results='asis'}
# Print OFL table
print(OFL.table, include.rownames = FALSE, caption.placement = 'top')
```