-
Notifications
You must be signed in to change notification settings - Fork 35
/
SSMethod.TA1.8.R
375 lines (367 loc) · 15.3 KB
/
SSMethod.TA1.8.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
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
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
#' Apply Francis composition weighting method TA1.8
#'
#' Uses method TA1.8 (described in Appendix A of Francis 2011) to do
#' stage-2 weighting of composition data from a Stock Synthesis model.
#' Outputs a multiplier, *w* (with bootstrap 95% confidence interval),
#' so that *N2y* = *w* x *N1y*, where *N1y* and
#' *N2y* are the stage-1 and stage-2
#' multinomial sample sizes for the data set in year y. Optionally
#' makes a plot of observed (with confidence limits, based on *N1y*)
#' and expected mean lengths (or ages).
#' \cr\cr
#' CAUTIONARY/EXPLANATORY NOTE.
#' The large number of options available in SS makes it very
#' difficult to be sure that what this function does is
#' appropriate for all combinations of options. The following
#' notes might help anyone wanting to check or correct the code.
#' \enumerate{
#' \item The code first takes the appropriate database (lendbase, sizedbase,
#' agedbase, or condbase) and removes unneeded rows.
#' \item The remaining rows of the database are grouped into individual
#' comps (indexed by vector indx) and relevant statistics (e.g.,
#' observed and expected mean length or age), and ancillary data,
#' are calculated for each comp (these are stored in pldat - one row
#' per comp).
#' If the data are to be plotted, the comps are grouped, with each
#' group corresponding to a panel in the plot, and groups are indexed
#' by plindx.
#' \item A single multiplier is calculated to apply to all the comps.
#' }
#'
#' @param fit Stock Synthesis output as read by r4SS function SS_output
#' @param type either 'len' (for length composition data), 'size' (for
#' generalized size composition data), 'age' (for age composition data),
#' or 'con' (for conditional age at length data)
#' @param fleet vector of one or more fleet numbers whose data are to
#' be analysed simultaneously (the output N multiplier applies
#' to all fleets combined)
#' @template fleetnames
#' @param part vector of one or more partition values; analysis is restricted
#' to composition data with one of these partition values.
#' Default is to include all partition values (0, 1, 2).
#' @param label.part Include labels indicating which partitions are included?
#' @param sexes vector of one or more values for Sexes; analysis is
#' restricted to composition data with one of these
#' Sexes values. Ignored if type=='con'.
#' @param label.sex Include labels indicating which sexes are included?
#' @param seas string indicating how to treat data from multiple seasons
#' 'comb' - combine seasonal data for each year and plot against Yr
#' 'sep' - treat seasons separately, plotting against Yr.S
#' If is.null(seas) it is assumed that there is only one season in
#' the selected data (a warning is output if this is not true) and
#' option 'comb' is used.
#' @param method a vector of one or more size-frequency method numbers
#' (ignored unless type = 'size').
#' If !is.null(method), analysis is restricted to size-frequency
#' methods in this vector. NB comps are separated by method
#' @param plotit if TRUE, make an illustrative plot like one or more
#' panels of Fig. 4 in Francis (2011).
#' @param printit if TRUE, print results to R console.
#' @param datonly if TRUE, don't show the model expectations
#' @param plotadj if TRUE, plot the confidence intervals associated with
#' the adjusted sample sizes (TRUE by default unless datonly = TRUE)
#' @param maxpanel maximum number of panels within a plot
#' @param set.pars Set the graphical parameters such as mar and mfrow.
#' Can be set to FALSE in order to add plots form multiple calls to
#' this function as separate panels in one larger figure.
#' @param add add to existing plot
#' @author R.I.C Chris Francis, Andre E. Punt, Ian G. Taylor
#' @export
#' @family tuning functions
#' @references Francis, R.I.C.C. (2011). Data weighting in statistical
#' fisheries stock assessment models. Canadian Journal of
#' Fisheries and Aquatic Sciences 68: 1124-1138.
#' @examples
#' \dontrun{
#' Nfleet <- length(myreplist[["FleetNames"]])
#' for (Ifleet in 1:Nfleet) {
#' SSMethod.TA1.8(myreplist, "len", fleet = Ifleet, maxpanel = maxpanel)
#' }
#' for (Ifleet in 1:Nfleet) {
#' SSMethod.TA1.8(myreplist, "age", fleet = Ifleet, maxpanel = maxpanel)
#' }
#' for (Ifleet in 1:Nfleet) {
#' SSMethod.TA1.8(myreplist, "size", fleet = Ifleet, maxpanel = maxpanel)
#' }
#' for (Ifleet in 1:Nfleet) {
#' SSMethod.TA1.8(myreplist, "con", fleet = Ifleet, maxpanel = maxpanel)
#' }
#' for (Ifleet in 1:Nfleet) {
#' SSMethod.Cond.TA1.8(myreplist, fleet = Ifleet, maxpanel = maxpanel)
#' }
#' }
#'
SSMethod.TA1.8 <-
function(fit, type, fleet, part = 0:2, sexes = 0:3, seas = NULL,
method = NULL, plotit = TRUE, printit = FALSE,
datonly = FALSE, plotadj = !datonly, maxpanel = 1000,
fleetnames = NULL, label.part = TRUE, label.sex = TRUE,
set.pars = TRUE, add = FALSE) {
# Check the type is correct and the sexes is correct
is.in <- function(x, y) !is.na(match(x, y))
if (!is.in(type, c("age", "len", "size", "con"))) {
stop('Illegal value for type (should be "age", "len", "size", or "con")')
} else {
if (sum(!is.in(sexes, c(0:3))) > 0) {
stop("Unrecognised value for sexes")
}
}
# replace default fleetnames with user input if requested
if (is.null(fleetnames)) {
# use fleetnames in the model
fleetnames <- fit[["FleetNames"]]
} else {
# if custom names input, check length
if (length(fleetnames) != fit[["nfleets"]]) {
stop("fleetnames needs to be NULL or have length = nfleets = ", fit[["nfleets"]])
}
}
# Select the type of datbase
dbase <- fit[[paste(type, "dbase", sep = "")]]
# sel is vector of row indices selected for the plot/calculations
# select row indices matching fleet and partition
sel <- is.in(dbase[["Fleet"]], fleet) & is.in(dbase[["Part"]], part)
# select row indices matching Sexes column
if (type != "con") {
# change column nanme on earlier SS versions to match change from
# Pick_sex to Sexes in 3.30.12 (July 2018)
names(dbase)[names(dbase) == "Pick_sex"] <- "Sexes"
sel <- sel & is.in(dbase$"Sexes", sexes)
}
# for generalized size frequency comps, select chosen size method
if (type == "size" & !is.null(method)) {
sel <- sel & is.in(dbase[["method"]], method)
}
# if there are no rows selected, return empty result
if (sum(sel) == 0) {
return()
}
# subset comp database for selected rows
dbase <- dbase[sel, ]
if (is.null(seas)) {
seas <- "comb"
if (length(unique(dbase[["Seas"]])) > 1) {
message("Combining data from multiple seasons")
}
}
# if generalized size comp is used, check for mix of units
if (type == "size") {
if (length(unique(dbase[["units"]])) > 1) {
warning("Mix of units being compared:", unique(dbase[["units"]]))
}
}
# create label for partitions
partitions <- sort(unique(dbase[["Part"]])) # values are 0, 1, or 2
partition.labels <- c("whole", "discarded", "retained")[partitions + 1]
partition.labels <- paste0("(", paste(partition.labels, collapse = "&"), " catch)")
# indx is string combining fleet, year, and potentially conditional bin
indx <- paste(dbase[["Fleet"]], dbase[["Yr"]], if (type == "con") {
dbase$"Lbin_lo"
} else {
""
}, if (seas == "sep") dbase[["Seas"]] else "")
# if subsetting by sex, add Sexes value to the indx strings
sex.flag <- type != "con" & max(tapply(
dbase$"Sexes",
dbase[["Fleet"]], function(x) length(unique(x))
)) > 1
if (sex.flag) {
indx <- paste(indx, dbase$"Sexes")
}
# if subsetting by generalized size-method, add that value to indx strings
method.flag <- type == "size" && length(unique(dbase[["method"]])) > 1
if (method.flag) {
indx <- paste(indx, dbase[["method"]])
}
# unique strings in indx vector
uindx <- unique(indx)
# test for length 1 results
if (length(uindx) == 1) {
# presumably the method is meaningless of there's only 1 point,
# but it's good to be able to have the function play through
message("Only one point to plot")
return()
}
# create empty data.frame to store information on each observation
pldat <- matrix(0, length(uindx), 10,
dimnames = list(
uindx,
c(
"Obsmn", "Obslo", "Obshi", "semn", "Expmn", "Std.res",
"ObsloAdj", "ObshiAdj", "Fleet", "Yr"
)
)
)
# add columns of zeros to fill with values necessary for subsetting
if (type == "con") pldat <- cbind(pldat, Lbin = 0)
if (sex.flag) pldat <- cbind(pldat, sexes = 0)
if (type == "size") {
pldat <- cbind(pldat, method = 0)
# vector to store units (which are strings and don't fit in pldat matrix)
plunits <- rep(NA, nrow(pldat))
}
# Find the weighting factor for this combination of factors
for (i in seq_along(uindx)) { # each row of pldat is an individual comp
subdbase <- dbase[indx == uindx[i], ]
xvar <- subdbase[["Bin"]]
# observed mean
pldat[i, "Obsmn"] <- sum(subdbase[["Obs"]] * xvar) / sum(subdbase[["Obs"]])
# expected mean
pldat[i, "Expmn"] <- sum(subdbase[["Exp"]] * xvar) / sum(subdbase[["Exp"]])
# use adjusted input sample size for Francis or MI weighting options
Nsamp <- subdbase[["Nsamp_adj"]]
if ("Nsamp_DM" %in% names(subdbase) || any(is.na(subdbase[["Nsamp_DM"]]))) {
# dirichlet multinomial newer format
Nsamp <- subdbase[["Nsamp_DM"]]
}
# standard error of the mean
pldat[i, "semn"] <- sqrt((sum(subdbase[["Exp"]] * xvar^2) / sum(subdbase[["Exp"]]) -
pldat[i, "Expmn"]^2) / mean(Nsamp))
# calculate confidence intervals and other stuff
pldat[i, "Obslo"] <- pldat[i, "Obsmn"] - 2 * pldat[i, "semn"]
pldat[i, "Obshi"] <- pldat[i, "Obsmn"] + 2 * pldat[i, "semn"]
pldat[i, "Std.res"] <- (pldat[i, "Obsmn"] - pldat[i, "Expmn"]) / pldat[i, "semn"]
pldat[i, "Fleet"] <- mean(subdbase[["Fleet"]])
pldat[i, "Yr"] <- mean(if (seas == "comb") subdbase[["Yr"]] else subdbase[["Yr.S"]])
if (type == "con") pldat[i, "Lbin"] <- mean(subdbase$"Lbin_lo")
if (sex.flag) {
pldat[i, "sexes"] <- mean(subdbase$"Sexes")
}
if (type == "size") {
pldat[i, "method"] <- mean(subdbase[["method"]])
plunits[i] <- subdbase[["units"]][1] # units of size comps
}
}
Nmult <- 1 / var(pldat[, "Std.res"], na.rm = TRUE)
# Find the adjusted confidence intervals
for (i in seq_along(uindx)) {
pldat[i, "ObsloAdj"] <- pldat[i, "Obsmn"] - 2 * pldat[i, "semn"] / sqrt(Nmult)
pldat[i, "ObshiAdj"] <- pldat[i, "Obsmn"] + 2 * pldat[i, "semn"] / sqrt(Nmult)
}
Nfleet <- length(unique(pldat[, "Fleet"]))
# make plot if requested
if (plotit) {
plindx <- if (type == "con") {
paste(pldat[, "Fleet"], pldat[, "Yr"])
} else {
pldat[, "Fleet"]
}
if (sex.flag) plindx <- paste(plindx, pldat[, "sexes"])
if (method.flag) plindx <- paste(plindx, pldat[, "method"])
uplindx <- unique(plindx)
# Select number of panels
Npanel <- length(uplindx)
## Ian T. 9/25/14: changing from having at least 4 panels to no minimum
# NpanelSet <- max(4,min(length(uplindx),maxpanel))
NpanelSet <- min(length(uplindx), maxpanel)
Nr <- ceiling(sqrt(NpanelSet))
Nc <- ceiling(NpanelSet / Nr)
if (set.pars) {
# save current graphical parameters
par_current <- par()
# set new parameters
par(
mfrow = c(Nr, Nc), mar = c(2, 2, 1, 1) + 0.1, mgp = c(0, 0.5, 0), oma = c(1.2, 1.2, 0, 0),
las = 1
)
par(cex = 1)
}
for (i in 1:Npanel) {
# loop over panels
subpldat <- pldat[plindx == uplindx[i], , drop = FALSE]
x <- subpldat[, ifelse(type == "con", "Lbin", "Yr")]
# make empty plot (unless adding to existing plot)
if (!add) {
# calculate ylim, including removing Inf values
plot(x, subpldat[, "Obsmn"],
pch = "-",
xlim = if (length(x) > 1) range(x) else c(x - 0.5, x + 0.5),
ylim = range(subpldat[, c("Obslo", "Obshi", "ObsloAdj", "ObshiAdj", "Expmn")],
finite = TRUE, na.rm = TRUE
),
xlab = "", ylab = ""
)
}
segments(x, subpldat[, "Obslo"], x, subpldat[, "Obshi"], lwd = 3, lend = 3)
if (plotadj) {
arrows(x, subpldat[, "ObsloAdj"], x, subpldat[, "ObshiAdj"],
lwd = 1,
length = 0.04, angle = 90, code = 3
)
}
points(x, subpldat[, "Obsmn"], pch = 21, bg = "grey80")
ord <- order(x)
if (!datonly) {
if (length(x) > 1) {
lines(x[ord], subpldat[ord, "Expmn"], lwd = 3, col = 4)
} else {
lines(c(x - 0.5, x + 0.5), rep(subpldat[, "Expmn"], 2), lwd = 3, col = 4)
}
}
# Lines
fl <- fleetnames[subpldat[1, "Fleet"]]
yr <- paste(subpldat[1, "Yr"])
lab <- if (type == "con") ifelse(Nfleet > 1, paste(yr, fl), yr) else fl
if (sex.flag & label.sex) {
lab <- paste(lab, ifelse(subpldat[1, "sexes"] == 0, "comb", "sex"))
}
if (method.flag) {
lab <- paste(lab, "meth", subpldat[1, "method"])
}
if (label.part) {
lab <- paste(lab, partition.labels)
}
if (!add) {
mtext(lab, side = 3, at = mean(x))
}
}
# define y-axis label
ylab <- "Mean age" # default as age unless replaced below
if (type == "len") {
ylab <- "Mean length"
}
if (type == "size") {
# probably more efficient ways to sort out these labels,
# but lots of if-statements make logic easier to follow
units <- unique(plunits[plindx %in% uplindx])
if (length(units) == 1) { # not sure if this will always be true or not
if (units %in% c("kg", "lb")) {
ylab <- paste0("Mean weight (", units, ")")
}
if (units %in% c("cm", "in")) {
ylab <- paste0("Mean length (", units, ")")
}
} else {
# just in case it's possible to have multiple units in one panel
ylab <- paste0("Mean value (", paste(units, collapse = " or "), ")")
}
}
if (!add) {
mtext(ylab, side = 2, las = 0, outer = TRUE)
mtext(ifelse(type == "con", "Length", "Year"), side = 1, outer = TRUE)
}
# restore previous graphics parameters (if changed to begin with
if (set.pars) {
par(
mfrow = par_current[["mfrow"]], mar = par_current[["mar"]], mgp = par_current[["mgp"]],
oma = par_current[["oma"]], las = par_current[["las"]]
)
}
}
if (!datonly) {
tmp <- matrix(sample(pldat[, "Std.res"], 1000 * nrow(pldat), replace = TRUE), nrow(pldat))
confint <- as.vector(quantile(apply(tmp, 2, function(x) 1 / var(x, na.rm = TRUE)),
c(0.025, 0.975),
na.rm = TRUE
))
Output <- c(w = Nmult, lo = confint[1], hi = confint[2])
Outs <- paste("Francis Weights - ", type, ": ", fleetnames[fleet], ": ",
round(Nmult, 4), " (", round(confint[1], 4), "-", round(confint[2], 4), ")",
sep = ""
)
if (printit) {
print(Outs)
}
return(Output)
}
}