sample_index: Sample the biomass with observation error

Description

This function creates an index of abundance sampled from the expected available biomass for given fleets in given years. Let B_y be the biomass from the operating model for year y. Then the sampled value is calculated as: B_y*exp(rnorm(1, 0, sds_obs)-sds_obs^2/2). The second term adjusts the random samples so that their expected value is B_y (i.e. the log-normal bias correction). If used with run_ss3sim the case file should be named index. A suggested (default) case letter is D for data.

Usage

sample_index(dat_list, outfile, fleets, years, sds_obs, make_plot = FALSE,
  write_file = TRUE)

Value

A modified .dat file if write_file=TRUE. A list object containing the modified .dat file is returned invisibly.

Which arguments to specifiy in case files

All function argument descriptions that start with an asterisk (*) will be passed through the case files to run_ss3sim. If one of these arguments is not specified in a case file, then a value of NULL will be passed, which may or may not be an appropriate value. Other arguments will be ignored if specified.

Examples

Run this code

# Find the example data location:
d <- system.file("extdata", package = "ss3sim")
f_in <- paste0(d, "/example-om/data.ss_new")
dat_list <- r4ss::SS_readdat(f_in, section = 2, verbose = FALSE)
dat_list <- change_fltname(dat_list)
outfile <- "test.dat"
ex1 <- sample_index(dat_list, outfile, fleets=c(2,3),
                    years=list(1938:2012, 1938:2012) ,
                    sds_obs=list(1e-6, 1e-6), write_file=FALSE,
                    make_plot = TRUE)
ex2 <- sample_index(dat_list, outfile, fleets=c(2,3),
                    years=list(1938:2012, 1938:2012) ,
                    sds_obs=list(.05, .05), write_file=FALSE,
                    make_plot = TRUE)
library(ggplot2)
ggplot(ex1, aes(x=year, y=obs, group=index, ymin=0,
                colour=as.factor(index)))+geom_line() + geom_point(data=ex2,
                aes(x=year, y=obs, colour=as.factor(index), group=index))
## Exclude a fleet and have varying sds_obs by year
ex3 <- sample_index(dat_list, outfile, fleets=c(2,NA),
                    years=list(1938:2012, 1950),
                    sds_obs=list(seq(.001, .1, len=75), .1),
                    write_file=FALSE)
ggplot(ex3, aes(x=year, y=obs, group=index, ymin=0,
                colour=as.factor(index)))+geom_point()

Run the code above in your browser using DataLab

Last chance! 50% off unlimited learning