In this module we will cover editing and running stock assessments using a script based workflow.
Why?
Containing all steps of an analysis within a script is a key component of an Openscience workflow
Makes analytic steps transparent
Improves reproducibility
Allows for integration with version control software (e.g., Git)
Facilitates documentation
Application
Modify a Stock Synthesis input model file and running the modified model through R using the r4ss package.
Specifically, in our example we will run models that make alternative assumptions for the steepness (\(h\)) parameter which controls the shape of the stock-recruitment relationship (SRR).
#V3.30.23.00;_safe;_compile_date:_Oct 31 2024;_Stock_Synthesis_by_Richard_Methot_(NOAA)_using_ADMB_13.2#_Stock_Synthesis_is_a_work_of_the_U.S._Government_and_is_not_subject_to_copyright_protection_in_the_United_States.#_Foreign_copyrights_may_apply._See_copyright.txt_for_more_information.#_User_support_available_at:NMFS.Stock.Synthesis@noaa.gov#_User_info_available_at:https://vlab.noaa.gov/group/stock-synthesis#_Source_code_at:_https://github.com/nmfs-ost/ss3-source-code#C growth parameters are estimated#C spawner-recruitment bias adjustment Not tuned For optimality#_data_and_control_files: data.ss // control.ss0# 0 means do not read wtatage.ss; 1 means read and use wtatage.ss and also read and use growth parameters1#_N_Growth_Patterns (Growth Patterns, Morphs, Bio Patterns, GP are terms used interchangeably in SS3)1#_N_platoons_Within_GrowthPattern #_Cond 1 #_Platoon_within/between_stdev_ratio (no read if N_platoons=1)#_Cond sd_ratio_rd < 0: platoon_sd_ratio parameter required after movement params.#_Cond 1 #vector_platoon_dist_(-1_in_first_val_gives_normal_approx)#4# recr_dist_method for parameters: 2=main effects for GP, Area, Settle timing; 3=each Settle entity; 4=none (only when N_GP*Nsettle*pop==1)1# not yet implemented; Future usage: Spawner-Recruitment: 1=global; 2=by area1# number of recruitment settlement assignments 0# unused option#GPattern month area age (for each settlement assignment)1110##_Cond 0 # N_movement_definitions goes here if Nareas > 1#_Cond 1.0 # first age that moves (real age at begin of season, not integer) also cond on do_migration>0#_Cond 1 1 1 2 4 10 # example move definition for seas=1, morph=1, source=1 dest=2, age1=4, age2=10#1#_Nblock_Patterns1#_blocks_per_pattern # begin and end years of blocks19701970## controls for all timevary parameters 1#_time-vary parm bound check (1=warn relative to base parm bounds; 3=no bound check); Also see env (3) and dev (5) options to constrain with base bounds## AUTOGEN00000# autogen: 1st element for biology, 2nd for SR, 3rd for Q, 4th reserved, 5th for selex# where: 0 = autogen time-varying parms of this category; 1 = read each time-varying parm line; 2 = read then autogen if parm min==-12345##_Available timevary codes#_Block types: 0: P_block=P_base*exp(TVP); 1: P_block=P_base+TVP; 2: P_block=TVP; 3: P_block=P_block(-1) + TVP#_Block_trends: -1: trend bounded by base parm min-max and parms in transformed units (beware); -2: endtrend and infl_year direct values; -3: end and infl as fraction of base range#_EnvLinks: 1: P(y)=P_base*exp(TVP*env(y)); 2: P(y)=P_base+TVP*env(y); 3: P(y)=f(TVP,env_Zscore) w/ logit to stay in min-max; 4: P(y)=2.0/(1.0+exp(-TVP1*env(y) - TVP2))#_DevLinks: 1: P(y)*=exp(dev(y)*dev_se; 2: P(y)+=dev(y)*dev_se; 3: random walk; 4: zero-reverting random walk with rho; 5: like 4 with logit transform to stay in base min-max#_DevLinks(more): 21-25 keep last dev for rest of years##_Prior_codes: 0=none; 6=normal; 1=symmetric beta; 2=CASAL's beta; 3=lognormal; 4=lognormal with biascorr; 5=gamma## setup for M, growth, wt-len, maturity, fecundity, (hermaphro), recr_distr, cohort_grow, (movement), (age error), (catch_mult), sex ratio #_NATMORT0#_natM_type:_0=1Parm; 1=N_breakpoints;_2=Lorenzen;_3=agespecific;_4=agespec_withseasinterpolate;_5=BETA:_Maunder_link_to_maturity;_6=Lorenzen_range#_no additional input for selected M option; read 1P per morph#1# GrowthModel: 1=vonBert with L1&L2; 2=Richards with L1&L2; 3=age_specific_K_incr; 4=age_specific_K_decr; 5=age_specific_K_each; 6=NA; 7=NA; 8=growth cessation0#_Age(post-settlement) for L1 (aka Amin); first growth parameter is size at this age; linear growth below this25#_Age(post-settlement) for L2 (aka Amax); 999 to treat as Linf-999#_exponential decay for growth above maxage (value should approx initial Z; -999 replicates 3.24; -998 to not allow growth above maxage)0#_placeholder for future growth feature#0#_SD_add_to_LAA (set to 0.1 for SS2 V1.x compatibility)0#_CV_Growth_Pattern: 0 CV=f(LAA); 1 CV=F(A); 2 SD=F(LAA); 3 SD=F(A); 4 logSD=F(A)#1#_maturity_option: 1=length logistic; 2=age logistic; 3=read age-maturity matrix by growth_pattern; 4=read age-fecundity; 5=disabled; 6=read length-maturity1#_First_Mature_Age1#_fecundity_at_length option:(1)eggs=Wt*(a+b*Wt);(2)eggs=a*L^b;(3)eggs=a*Wt^b; (4)eggs=a+b*L; (5)eggs=a+b*W0#_hermaphroditism option: 0=none; 1=female-to-male age-specific fxn; -1=male-to-female age-specific fxn1#_parameter_offset_approach for M, G, CV_G: 1- direct, no offset**; 2- male=fem_parm*exp(male_parm); 3: male=female*exp(parm) then old=young*exp(parm)#_** in option 1, any male parameter with value = 0.0 and phase <0 is set equal to female parameter##_growth_parms#_ LO HI INIT PRIOR PR_SD PR_type PHASE env_var&link dev_link dev_minyr dev_maxyr dev_PH Block Block_Fxn# Sex: 1 BioPattern: 1 NatMort0.050.150.10.10.80-30000000# NatM_uniform_Fem_GP_1# Sex: 1 BioPattern: 1 Growth-104520.937236106-20000000# L_at_Amin_Fem_GP_1409071.556670106-40000000# L_at_Amax_Fem_GP_10.050.250.1636780.150.86-40000000# VonBert_K_Fem_GP_10.050.250.10.10.80-30000000# CV_young_Fem_GP_10.050.250.10.10.80-30000000# CV_old_Fem_GP_1# Sex: 1 BioPattern: 1 WtLen-332.44e-062.44e-060.80-30000000# Wtlen_1_Fem_GP_1-343.346943.346940.80-30000000# Wtlen_2_Fem_GP_1# Sex: 1 BioPattern: 1 Maturity&Fecundity506055550.80-30000000# Mat50%_Fem_GP_1-33-0.25-0.250.80-30000000# Mat_slope_Fem_GP_1-33110.80-30000000# Eggs/kg_inter_Fem_GP_1-33000.80-30000000# Eggs/kg_slope_wt_Fem_GP_1# Sex: 2 BioPattern: 1 NatMort0.050.150.10.10.80-30000000# NatM_uniform_Mal_GP_1# Sex: 2 BioPattern: 1 Growth045036100-30000000# L_at_Amin_Mal_GP_1409069.74840070106-40000000# L_at_Amax_Mal_GP_10.050.250.1735160.150.86-40000000# VonBert_K_Mal_GP_10.050.250.10.10.80-30000000# CV_young_Mal_GP_10.050.250.10.10.80-30000000# CV_old_Mal_GP_1# Sex: 2 BioPattern: 1 WtLen-332.44e-062.44e-060.80-30000000# Wtlen_1_Mal_GP_1-343.346943.346940.80-30000000# Wtlen_2_Mal_GP_1# Hermaphroditism# Recruitment Distribution # Cohort growth dev base0.1101110-10000000# CohortGrowDev# Movement# Platoon StDev Ratio # Age Error from parameters# catch multiplier# fraction female, by GP1e-060.9999990.50.50.50-990000000# FracFemale_GP_1# M2 parameter for each predator fleet##_no timevary MG parameters##_seasonal_effects_on_biology_parms0000000000#_femwtlen1,femwtlen2,mat1,mat2,fec1,fec2,Malewtlen1,malewtlen2,L1,K#_ LO HI INIT PRIOR PR_SD PR_type PHASE#_Cond -2 2 0 0 -1 99 -2 #_placeholder when no seasonal MG parameters#3#_Spawner-Recruitment; Options: 1=NA; 2=Ricker; 3=std_B-H; 4=SCAA; 5=Hockey; 6=B-H_flattop; 7=survival_3Parm; 8=Shepherd_3Parm; 9=RickerPower_3parm0# 0/1 to use steepness in initial equ recruitment calculation0# future feature: 0/1 to make realized sigmaR a function of SR curvature#_ LO HI INIT PRIOR PR_SD PR_type PHASE env-var use_dev dev_mnyr dev_mxyr dev_PH Block Blk_Fxn # parm_name3311010.310010000000# SR_LN(R0)0.210.5738350.70.051-40000000# SR_BH_steep020.60.80.80-40000000# SR_sigmaR-550010-40000000# SR_regime000000-990000000# SR_autocorr#_no timevary SR parameters2#do_recdev: 0=none; 1=devvector (R=F(SSB)+dev); 2=deviations (R=F(SSB)+dev); 3=deviations (R=R0*dev; dev2=R-f(SSB)); 4=like 3 with sum(dev2) adding penalty1971# first year of main recr_devs; early devs can precede this era2001# last year of main recr_devs; forecast devs start in following year2#_recdev phase 1# (0/1) to read 13 advanced options0#_recdev_early_start (0=none; neg value makes relative to recdev_start)-4#_recdev_early_phase0#_forecast_recruitment phase (incl. late recr) (0 value resets to maxphase+1)1#_lambda for Fcast_recr_like occurring before endyr+11949.0#_last_early_yr_nobias_adj_in_MPD 1984.0#_first_yr_fullbias_adj_in_MPD 1997.0#_last_yr_fullbias_adj_in_MPD 2001.9#_first_recent_yr_nobias_adj_in_MPD 0.528#_max_bias_adj_in_MPD (1.0 to mimic pre-2009 models) 0#_period of cycles in recruitment (N parms read below)-5#min rec_dev5#max rec_dev0#_read_recdevs#_end of advanced SR options##_placeholder for full parameter lines for recruitment cycles# read specified recr devs#_year Input_value## all recruitment deviations# 1971R 1972R 1973R 1974R 1975R 1976R 1977R 1978R 1979R 1980R 1981R 1982R 1983R 1984R 1985R 1986R 1987R 1988R 1989R 1990R 1991R 1992R 1993R 1994R 1995R 1996R 1997R 1998R 1999R 2000R 2001R 2002F 2003F 2004F 2005F 2006F 2007F 2008F 2009F 2010F 2011F# 0.134968 -0.0507584 0.0972615 -0.164709 0.0443745 0.710789 -0.00439022 0.0113464 0.266822 0.186805 0.100984 -0.208826 -0.418054 -0.284269 0.420545 0.586619 0.262575 0.204289 -0.321237 0.664807 -0.607033 -0.193662 -0.7506 0.446329 -0.524044 0.539563 1.20007 -0.459983 -0.567356 0.263829 -0.216963 0 0 0 0 0 0 0 0 0 0##Fishing Mortality info 0.3# F ballpark value in units of annual_F-2001# F ballpark year (neg value to disable)3# F_Method: 1=Pope midseason rate; 2=F as parameter; 3=F as hybrid; 4=fleet-specific parm/hybrid (#4 is superset of #2 and #3 and is recommended)2.95# max F (methods 2-4) or harvest fraction (method 1)4# N iterations for tuning in hybrid mode; recommend 3 (faster) to 5 (more precise if many fleets)##_initial_F_parms; for each fleet x season that has init_catch; nest season in fleet; count = 0#_for unconstrained init_F, use an arbitrary initial catch and set lambda=0 for its logL#_ LO HI INIT PRIOR PR_SD PR_type PHASE## F rates by fleet x season#_year: 1971 1972 1973 1974 1975 1976 1977 1978 1979 1980 1981 1982 1983 1984 1985 1986 1987 1988 1989 1990 1991 1992 1993 1994 1995 1996 1997 1998 1999 2000 2001 2002 2003 2004 2005 2006 2007 2008 2009 2010 2011# seas: 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1# FISHERY 0 0.00211813 0.0106457 0.0107404 0.0217792 0.0334421 0.046095 0.0601249 0.0759325 0.10803 0.147259 0.162937 0.181299 0.203352 0.230858 0.266731 0.315247 0.338831 0.3551 0.356617 0.339428 0.2384 0.243238 0.251023 0.263884 0.283723 0.227427 0.23848 0.247863 0.25268 0.253549 0.0132306 0.0281356 0.0381987 0.0448672 0.0493659 0.0525756 0.0550896 0.0572622 0.0592695 0.0611758##_Q_setup for fleets with cpue or survey or deviation data#_1: fleet number#_2: link type: 1=simple q; 2=mirror; 3=power (+1 parm); 4=mirror with scale (+1p); 5=offset (+1p); 6=offset & power (+2p)#_ where power is applied as y = q * x ^ (1 + power); so a power value of 0 has null effect#_ and with the offset included it is y = q * (x + offset) ^ (1 + power)#_3: extra input for link, i.e. mirror fleet# or dev index number#_4: 0/1 to select extra sd parameter#_5: 0/1 for biasadj or not#_6: 0/1 to float#_ fleet link link_info extra_se biasadj float # fleetname210100# SURVEY1310000# SURVEY2-999900000##_Q_parameters#_ LO HI INIT PRIOR PR_SD PR_type PHASE env-var use_dev dev_mnyr dev_mxyr dev_PH Block Blk_Fxn # parm_name-750.51691601010000000# LnQ_base_SURVEY1(2)00.500.0510-40000000# Q_extraSD_SURVEY1(2)-105-6.6254701010000000# LnQ_base_SURVEY2(3)#_no timevary Q parameters##_size_selex_patterns#Pattern:_0; parm=0; selex=1.0 for all sizes#Pattern:_1; parm=2; logistic; with 95% width specification#Pattern:_5; parm=2; mirror another size selex; PARMS pick the min-max bin to mirror#Pattern:_11; parm=2; selex=1.0 for specified min-max population length bin range#Pattern:_15; parm=0; mirror another age or length selex#Pattern:_6; parm=2+special; non-parm len selex#Pattern:_43; parm=2+special+2; like 6, with 2 additional param for scaling (mean over bin range)#Pattern:_8; parm=8; double_logistic with smooth transitions and constant above Linf option#Pattern:_9; parm=6; simple 4-parm double logistic with starting length; parm 5 is first length; parm 6=1 does desc as offset#Pattern:_21; parm=2*special; non-parm len selex, read as N break points, then N selex parameters#Pattern:_22; parm=4; double_normal as in CASAL#Pattern:_23; parm=6; double_normal where final value is directly equal to sp(6) so can be >1.0#Pattern:_24; parm=6; double_normal with sel(minL) and sel(maxL), using joiners#Pattern:_2; parm=6; double_normal with sel(minL) and sel(maxL), using joiners, back compatibile version of 24 with 3.30.18 and older#Pattern:_25; parm=3; exponential-logistic in length#Pattern:_27; parm=special+3; cubic spline in length; parm1==1 resets knots; parm1==2 resets all #Pattern:_42; parm=special+3+2; cubic spline; like 27, with 2 additional param for scaling (mean over bin range)#_discard_options:_0=none;_1=define_retention;_2=retention&mortality;_3=all_discarded_dead;_4=define_dome-shaped_retention#_Pattern Discard Male Special1000# 1 FISHERY1000# 2 SURVEY10000# 3 SURVEY2##_age_selex_patterns#Pattern:_0; parm=0; selex=1.0 for ages 0 to maxage#Pattern:_10; parm=0; selex=1.0 for ages 1 to maxage#Pattern:_11; parm=2; selex=1.0 for specified min-max age#Pattern:_12; parm=2; age logistic#Pattern:_13; parm=8; age double logistic. Recommend using pattern 18 instead.#Pattern:_14; parm=nages+1; age empirical#Pattern:_15; parm=0; mirror another age or length selex#Pattern:_16; parm=2; Coleraine - Gaussian#Pattern:_17; parm=nages+1; empirical as random walk N parameters to read can be overridden by setting special to non-zero#Pattern:_41; parm=2+nages+1; // like 17, with 2 additional param for scaling (mean over bin range)#Pattern:_18; parm=8; double logistic - smooth transition#Pattern:_19; parm=6; simple 4-parm double logistic with starting age#Pattern:_20; parm=6; double_normal,using joiners#Pattern:_26; parm=3; exponential-logistic in age#Pattern:_27; parm=3+special; cubic spline in age; parm1==1 resets knots; parm1==2 resets all #Pattern:_42; parm=2+special+3; // cubic spline; with 2 additional param for scaling (mean over bin range)#Age patterns entered with value >100 create Min_selage from first digit and pattern from remainder#_Pattern Discard Male Special000000000000##_ LO HI INIT PRIOR PR_SD PR_type PHASE env-var use_dev dev_mnyr dev_mxyr dev_PH Block Blk_Fxn # parm_name# 1 FISHERY LenSelex198053.6629500.01120000000# Size_inflection_FISHERY(1)0.016018.9332150.01130000000# Size_95%width_FISHERY(1)# 2 SURVEY1 LenSelex197036.6585300.01120000000# Size_inflection_SURVEY1(2)0.01606.59966100.01130000000# Size_95%width_SURVEY1(2)# 3 SURVEY2 LenSelex#_No_Dirichlet parameters#_no timevary selex parameters#0# use 2D_AR1 selectivity? (0/1)#_no 2D_AR1 selex offset used#_specs: fleet, ymin, ymax, amin, amax, sigma_amax, use_rho, len1/age2, devphase, before_range, after_range#_sigma_amax>amin means create sigma parm for each bin from min to sigma_amax; sigma_amax<0 means just one sigma parm is read and used for all bins#_needed parameters follow each fleet's specifications# -9999 0 0 0 0 0 0 0 0 0 0 # terminator## Tag loss and Tag reporting parameters go next0# TG_custom: 0=no read and autogen if tag data exist; 1=read#_Cond -6 6 1 1 2 0.01 -4 0 0 0 0 0 0 0 #_placeholder if no parameters## no timevary parameters### Input variance adjustments factors: #_1=add_to_survey_CV#_2=add_to_discard_stddev#_3=add_to_bodywt_CV#_4=mult_by_lencomp_N#_5=mult_by_agecomp_N#_6=mult_by_size-at-age_N#_7=mult_by_generalized_sizecomp#_factor fleet value-999910# terminator#4#_maxlambdaphase1#_sd_offset; must be 1 if any growthCV, sigmaR, or survey extraSD is an estimated parameter# read 3 changes to default Lambdas (default value is 1.0)# Like_comp codes: 1=surv; 2=disc; 3=mnwt; 4=length; 5=age; 6=SizeFreq; 7=sizeage; 8=catch; 9=init_equ_catch; # 10=recrdev; 11=parm_prior; 12=parm_dev; 13=CrashPen; 14=Morphcomp; 15=Tag-comp; 16=Tag-negbin; 17=F_ballpark; 18=initEQregime#like_comp fleet phase value sizefreq_method122114221142311-99991111# terminator## lambdas (for info only; columns are phases)# 0 0 0 0 #_CPUE/survey:_1# 1 1 1 1 #_CPUE/survey:_2# 1 1 1 1 #_CPUE/survey:_3# 1 1 1 1 #_lencomp:_1# 1 1 1 1 #_lencomp:_2# 0 0 0 0 #_lencomp:_3# 1 1 1 1 #_agecomp:_1# 1 1 1 1 #_agecomp:_2# 0 0 0 0 #_agecomp:_3# 1 1 1 1 #_size-age:_1# 1 1 1 1 #_size-age:_2# 0 0 0 0 #_size-age:_3# 1 1 1 1 #_init_equ_catch1# 1 1 1 1 #_init_equ_catch2# 1 1 1 1 #_init_equ_catch3# 1 1 1 1 #_recruitments# 1 1 1 1 #_parameter-priors# 1 1 1 1 #_parameter-dev-vectors# 1 1 1 1 #_crashPenLambda# 0 0 0 0 # F_ballpark_lambda0# (0/1/2) read specs for more stddev reporting: 0 = skip, 1 = read specs for reporting stdev for selectivity, size, and numbers, 2 = add options for M,Dyn. Bzero, SmryBio# 1 1 -1 5 # Selectivity: (1) 0 to skip or fleet, (2) 1=len/2=age/3=combined, (3) year, (4) N selex bins; NOTE: combined reports in age bins# 1 5 # Growth: (1) 0 to skip or growth pattern, (2) growth ages; NOTE: does each sex# 1 -1 5 # Numbers-at-age: (1) 0 or area(-1 for all), (2) year, (3) N ages; NOTE: sums across morphs# 5 15 25 35 43 # vector with selex std bins (-1 in first bin to self-generate)# 1 2 14 26 40 # vector with growth std ages picks (-1 in first bin to self-generate)# 1 2 14 26 40 # vector with NatAge std ages (-1 in first bin to self-generate)999
In this case we also copy the executable into our model directory.
ss3-example-calcs.r
# copy over executable dir_exec =paste0(proj_dir,"/executables/stock-synthesis/3.30.23.1/") ss3_exec ="ss3_linux"file.copy(from=paste0(dir_exec,ss3_exec),to=dir_base_stock_synthesis)# run the modelrun(dir=dir_base_stock_synthesis,exe=ss3_exec,show_in_console=TRUE,skipfinished=FALSE)
Testing multiple alternatives
What if we wanted to test models with several different alternatives for \(h\)?
Say also testing \(h = 0.5\) and \(h = 0.6\)?
Put it in a for loop!
ss3-example-calcs.r
# modify input file & run new stock synthesis models steepness_vec =c(0.5,0.6,0.7)for(i inseq_along(steepness_vec)){# make new directory tmp_dir =paste0(dir_model,"steepness-",steepness_vec[i],"/")dir.create(tmp_dir,recursive=TRUE)# modify control file# read in baseline stock synthesis files using r4ss functions tmp_ctl =SS_readctl(file=paste0(dir_base_stock_synthesis,"control.ss_new"),datlist =paste0(dir_base_stock_synthesis,"data_echo.ss_new")) tmp_ctl$SR_parms["SR_BH_steep","INIT"] = steepness_vec[i]# write out file using r4ss functionsSS_writectl(tmp_ctl,outfile=paste0(tmp_dir,"control.ss"),overwrite=TRUE)# copy remaining input files and executablefile.copy(from=paste0(dir_base_stock_synthesis,c("starter.ss","data.ss","forecast.ss")),to=paste0(tmp_dir))file.copy(from=paste0(dir_exec,ss3_exec),to=tmp_dir)# run the modelrun(dir=tmp_dir,exe=ss3_exec,show_in_console=TRUE,skipfinished=FALSE)# clean-up workspacerm(list=c("tmp_dir","tmp_ctl")) }
Group activity!
Using the toy Stock Synthesis model in the stock-synthesis-models/base-model directory make some modifications to the model, and run the different configurations.