##################################################################################### ### PROGRAM FOR SIMULATING PREDATOR PREY POPULATIONS ### by: Andre Buchheister ### First Created: 12-15-11 ### ##################################################################################### # DESCRIPTION # For Chapter 4 of the dissertation, a simulation will be conducted looking at the # influence of prey abundances on predator production. See prospectus for more details, # but in general, this simulation model is comprised of a functional response and growth # component (where consumption and growth are dictated by prey abundances), and a population # component (where the predator population is modeled with recruitment, growth, mortality, # and migraiton). Scenarios are set up to investigate two main objectives: # 1. How do simulated changes in the abundances of modeled prey groups influence # the consumption and productivity of summer flounder at an individual- and population-scale. # 2. How do patterns in predatory migration and connectivity among ecosystems mitigate the # effects on summer flounder population dynamics? #Notes: #1/3/11 - Working on expanding the model to 2 regions. Note, initially, not making immigration (ie phi terms) age-specific #5/30/12 - Added stochastic recruitment variability, and brought in the Functional Response component. # - Consumption of the predator are driven by forced prey abundance, #6/11/12 - Changed ages of simulation to match up with assessment (0 to 7+) #6/12/12 - Trying to add in an optimization routine to estimate reasonable attack rate parameters to get reasonable weight-at-age #6/22/12 - Modifying growth function to be consistent with Mike Wilberg's model (Wt = Wt-1 + Gmax(Ct/Cmax)) #6/25/12 - Started to make attack rate parameters variable with fish ages. First attempt is a simple linear relationship # where atck rate on prey 2 (anchovy) increases with age. #7/8/12 - Making attack rate a logistic function of age, with ages being fractional through the seasons. #7/11/12 - Making a giant function to run the simulation. # - 11pm - Automating the simfxn() to calculate the percent differences for 2 scenarios simultaneously (Prey 1 pulse; prey 2 pulse) #7/12/12 - Made prey pulses only occur in Region 1, and set up "migratory" and "fully mixed" migration scenarios # Immigration made age dependent. #1/29/13 - Modified Cmax to reflect the higher specific rates of consumption by smaller fish # (following Wisconsin Manual and Hartman and Brandt 1995 (striped bass) parameters) # - Modified predator weight array to account for the mixing of fishes among regions. The mean weight of the fish in a region at # a time step will be a weighted mean of the fish that migrated into the area and of those that stayed. # - IE: w.pred = (N1*w1 + N2*w2) / (N1+N2) where N is the number of fish coming (or staying from a region), and # w is the weight-at-age of fish growing in each region. #1/29/13 - Modifying code to add a 3rd prey species... acting as a large fish prey. # - The 3 prey model treats the prey as sequentially preferred prey by predators. The trends in ChesMMAP diets # for crustacean, small fish prey, and large fish prey were used to guide those relationships. (see Flounder_diets by age_1-29-13.xlsx) # - Recoding the prey size and abundances to be relative to one another. (see Flounder_diets by age_1-29-13.xlsx). # To facilitate this, i'm treating all prey as having an individual wt of 1g. That way the relative differences # in prey number are equivalent to the relative biomass, and the biomass is ultimately what dictates the FxnR calculation. #1/30/13 - Added the attackrates() function to optimize the functional response parameters to best fit the weight at age and the # observed diet patterns from ChesMMAP # - To optimize to diet patterns, i made functions to represent the diet patterns (see Flounder_diets by age_1-29-13.xlsx) #2/1/13 - Calculating annual fishery yields #2/2/13 - Updated the stock-recruitment function (using Rothschild et al. 2012 data and formulation, which included a shape parameter beta) # - Also fixed some plotting issues that existed for the S-R plot. #2/4/13 - Modifying this code to be able to run sensitivity analyses #2/5/13 - Worked on doing a better job optimizing the functional response parameters. Had difficulty with the optim() function to work # correctly (don't know why). Ended up making some changes to the growthfxn() and attackrates(): # - Instead of RSS as the metric to minimize, i switched to a Mean Square Error (MSE) measure that was the sum of the MSE for the # weight at age relationship AND the %diet. # - After getting the optimization functions to work, i ended up tweaking some parameters prior to optimization: # - Gmax.a now = rep(0.6,ages) # - atck.infl.3 = 3 (instead of 2) # - Decided on a new set of BASE, default parameter values for the simulation that better matched the w@age and diet patterns. # - Sensitivity analyses conducted and values stored in the Output folder #2/7/13 - Adding fishing scenarios: High, Low, and Declining #2/8/13 - Calculated time to reaching a target SSB level # - Building tihs Time to target output into the simulation function. #2/11/13 - uncovered a problem in which the PredPreySim() function does not gerenate the same SSB estimates if years=60 vs. 80. # - The values should be the same for the first 60 years in common... # - I found the error: the N.prey array was being filled in by columns, not rows! This problem has been fixed as of PredPreySim_2-11-13.r; # the solution was to use aperm() to transpose the N.prey array, fill it in, then transpose it back. # - Also, in this file i've now standardized the Pulse period to be for 25 years. So, this means there is a 25 year burn in period, # then a 25 year "pulse" period (if dictated by P.scen), then the remaining years are simply to reach equilibrium, and no pulses occur. #3/13/13 - Looking into Mike Wilberg's comment regarding whether there is a benefit to the "Pulse" situation vs. a constant increase. # Created a "Step" scenario whereby prey biomass increases by a constant amount during years 26:50. # This increase is an even redistribution of the added productivity created in the pulsed scenarios. # This redistribution or step increase equates to ~1.5-2.5x mean increase in production (see ProdStepIncr ) #3/18/13 - Reworking the simfxn output to behave better with the new step scenarios. # - Changed w.diff.mean to "w.diff.p.mean" to be consistent with the other outputs (S.diff.p.mean, Y.diff.p.mean) where "p" denotes # % diff values calculated for pulse years only. This is important so that i can also calculate % diff across all years for weight at age # for the Step Scenarios. (NOTE: THESE CHANGES NEED TO BE DUPLICATED IN THE GRAPHING CODE) # - TO FIGURE OUT - why are the pulses.all years identical? - decide if output files should be calculated for years 26:50 or 26:years #8/26/13 - After further exploring the implementation of a bioenergetics model to calculate growth within the simulation, (see PredPreySim_BEM-8-21-13.r and BioEnModel_Lit_params_8-15-13.xlsx) # it became evident that the BEM is VERY sensitive to parameter inputs, particularly P. This easily led to explosions in growth. # So, when P (C/Cmax) is changed with the functional response, the model is likely to blow up. # - Ultimately, i decided to back away from the BEM approach and return to the G/Gmax=C/Cmax (ie conversion # efficiency) approach. See GrowthModel_ConvEfficiency_8-25-13.xlsx for more detail on obtaining the parameter values # that will be used within the simulation. # - I got the PredPreySim() function working properly after some minor de-bugging #8/27/13 - Working on optimizing the attack rate parameters to match the diets and growth. #8/29/13 - Updating the Sensitivity function. #8/30/13 - Created Fmax, Fmin, and decrF.time parameters to control the fishing mortality rates. Also, i wanted to be more # in line with the F trends from StockRecruitFxn_8-29-13.xlsx, so i made Fmax 1.5 (instead of 2), and decrF.time = 12 years #9/3/13 - Changed Cmax to be dependent on the previous season's weight. (how much of an effect will this have? Probably note much. #9/4/13 - Re-evaluating the migration (phi.in) parameters and simplified the values make season specific only (ie removing the age-dependence) # - Also, changed recruitment to be 50-50 by region in the mixed scenario. # - Figure out why some sensitivity results are yielding 0% difference from the BASE model!!!! # - Is it due to minute differences that disapear with rounding? Doesn't appear to be the case. (Took the rounding out to 10, and it showed results to 5 decimal figures that were identical. # - Turns out it was due to a coding issue in the simfxn() where i didn't add in the "new" parameters (CA, KR, atck.max, etc.) to the embedded PredPreySim() function. #9/5/13 - Updated offshore temps based on NOAA survey data (See NOAA_Temp Data_offshore_9-5-13.xlsx). Note that the seasonal # trend may be a little off due to aliasing and missing monthly values. # - Decided NOT to make Cmax a function of actual weight (instead of mean empirical weight). # - this would alter the model structure as I have it, and i don't want to make such a large change at this stage. # The changes to the results would be highly negligible. #9/6/13 - Changed s.spawn to be in fall (season 4) once i realized i had the seasons wrongly defined before (correct: 1-winter, 2-spring, 3-summer, 4-fall). #9/7/13 - Deleted "Time To Target" code section that was after the Sim Runs, since it was old code. # - Worked on building in a "Monte-carlo-type" sensitivity analysis #9/9/13 - Updated sensitivity plotting code to plot time to target sensitivities using method 2 (ind paramater perturbation) #9/13/13 - Explored an alternative way to determine Fmsy and Smsy, by running PredPreySim() at different Fs to find the one that generates the biggest yield # - In doing so, i also found that the S.target should really be calculated for the new spawning season (s=4). # This was fixed by makeing the reference to s.spawn directly. #9/14/13 - Built different F scenarios that represent moratorium or that replicates the F trajectory given by the stock assessment. #9/19/13 - Built additional F scenarios that represent a knife-edge change from Fmax to Fmin "knifF"; # - Also added an exponential decrease in F scenario to Fmin ("exp.F"). # - Fmin changed to 0.31 (from 0.255) because this equates to FMSY, based on simulation runs without stochasticity. # - Realized that the SR parameters were not set as function parameters, but instead hard-coded. Fixed this. It will have some # impact on the sensitivity analyses... #9/20/13 - Added a decrF scenario where the linear decrease occurs in half the time of decrF.time # - Realized that by have pulse.range = c(2,10) that this led to instances where the prey biomass was 30X greater # than the mean prey biomass. This was because the multiplier could amplify the natural stochastic variability # created by the prey.dev.sd = 0.4. # - To alleviate this problem i explored calculating the prey pulses as an additive effect (N(t) = N(t) + N.mean*Multiplier. # - I decided that this additive approach was not the way to go as it completely muted the random seasonal pattern # and is likely not how the prey dynamics operate. # - I decided to reduce pulse.range so that the prey biomass pulses weren't so drastic. Pulse.range = c(2:6) seemed # to generate more reasonable values where the mean of the max biomass increases across sims was ~10X. # See "PulseRange effect on biomass pulses_Exploration_9-24-13.xlsx" for values, or look at code section following the PredPreySim() code (before simfxn() section) #9/24/13 - Re-set the default pulse.range to c(2:6) and re-ran the simulation... again. Lord, help me. # - Realized that the S, Y, and w outputs from simfxn() for all years were being calculated from 26:years, instead of 26:50 (ie the years # when a pulse COULD have occurred). Fixed this but need to re-run simulation. Did this at 1240am 9-25-13. #9/30/13 - Decided to add RECRUITMENT and %Change in R (in the year following the pulse) as an output of the simulation function. # The simulation was re-run to generate this new data for plotting and possible inclusion in the manuscript. #10/5/13 - Updated the code for the sensitivity plots so that the model parameters were expressed with the official greek letters and subscripts # - Rerunning the MC Sensitivity analysis to get the Recruitment output. #7/23/14 - Finally getting back around to this Chapter to get it published. Based on Tom Miller's suggestion, i'm going to re-run the simulation with # more than 100 simulations... aiming for 1000. Also, now i'm doing this from my CBL laptop, so some code may change due to that as well. # - This file builds off of version PredPreySim_ce_9-24-13.r which was the version used for the dissertation. # INITIAL CONDITIONS/DESCRIPTION # - 1 predator, 8 age-classes, with plus group (0-7+) # - 3 prey (with forced abundances) # - prey 1 = crustaceans (mysids & shrimps) # - prey 2 = small fish (e.g. bay anchovy) # - prey 3 = larger fish (e.g. spot, scup) # - 2 Regions # - seasonal time step over 50 years ### ORDER OF OPERATIONS # Within each timestep: # [ RECRUITMENT - MORTALITY - GROWTH - MOVEMENT - CENSUS ] # Notes: - Thus, the abundance at time t is the abundance at the END of the time step # - also, note that the CENSUS is for both the P (or abundance) matrix as well as for SSB library(reshape2) #for melt() library(ggplot2) #for ggplot() library(plyr) #for ddply() library(Hmisc) #for minor.tick() ### IMPORT MISCELLANEOUS DATA #Importing LW and age data for flounder from ChesMMAP LWdata=read.csv("C:\\Users\\andre\\VIMS\\PhD\\Research\\PredPrey_Simulations\\Data\\Flounder_LW_Age_6-11-12.csv") #Importing SR and F data for flounder from assessment SRdata=read.csv("C:\\Users\\andre\\VIMS\\PhD\\Research\\PredPrey_Simulations\\Data\\SR_data_from_assessment_9-14-13.csv") SRdata = data.frame(SimYear = c(11:39),SRdata) SRdata ## SIMULATION BOUNDS/CONDITIONS #Numbers of predators, ages, prey, seasons, years, and regions in simulation predators <<- 1 ages <<- 8 #represents ages 0 to 7+ to be consistent with stock assessment preys <<- 3 #prey 1 = Crustaceans (mysid), prey 2 = Small fishes (anchovy), 3 = larger fishes (spot) years <<- 80 years.2 = years #This is for a coding nuance... necessary when estimating the target SSB within the sim function seasons <<- 4 regions <<- 2 #Region 1 is the Chesapeake Bay ############################################################################################################ ############################################################################################################ ### PREDATOR / PREY SIMULATION FUNCTION ############################################################################################################ ############################################################################################################ PredPreySim = function(seed=500, #set the random seed value M=0.25, #Instantaneous natural mortality rate (per year) Fmax = 1.5, #Maximum instantaneous fishing mortality. Previously 2.0; exploring how 1.5 looks since that was more the average. (See StockRecruitFxn_8-29-13.xlsx) Fmin = 0.31, #Minimum instantaneous fishing mortality (equal to Fmsy for my simulation; see code following the PredPreySim() definition. decrF.time = 12, #Number of years to go from Fmax to Fmin (F=0.255 is the target F as of Jan 2013). F.sel = c(0,0.1,0.5,rep(1.0,ages-3)), #Fishing selectivity by age R.dev.sd=0.4, #Standard deviation for the lognormal random recruitment error (SD=0.4 see Beddington and Cook, 1983 from Pat Lynch; note Maunder 2012 uses 0.6) alpha.SR = 3.4, #alpha parameter for Shepherd S-R function (Rothschild et al. 2012) K.SR = 27.4, #K parameter for Shepherd S-R function (Rothschild et al. 2012) beta.SR = 1, #beta shape parameter for Shepherd S-R function. beta.SR=1 is for Beverton Holt. #OLD Gmax.a = rep(0.6,ages), #Age-specific maximum growth (kg/yr) possible for species CTO= 22 , #Temperature for optimum consumption; assumed value (see eq in Hanson et al. 1997) CTM= 35 , #Maximum temp parameter for temperature equation; assumed value (see eq in Hanson et al. 1997) CQ = 2.5 , #Rate parameter for temperature equation; assumed value (see eq in Hanson et al. 1997) CA=0.3 , #Intercept Parameter for allometric relationship of Cmax (Hanson et al. 1997; GrowthModel_ConvEfficiency_8-25-13.xlsx) Cmax = CA*W^CB CB=-0.2 , #Decay Parameter for allometric relationship of Cmax (Hanson et al. 1997; GrowthModel_ConvEfficiency_8-25-13.xlsx) Cmax = CA*W^CB KL = 0.5 , #Lambda parameter for maximum conversion effiency (aka Kmax) sigmoidal model(GrowthModel_ConvEfficiency_8-25-13.xlsx) KR = -0.0014 , #Rate of decline for sigmoidal K model (GrowthModel_ConvEfficiency_8-25-13.xlsx) KW = -600 , #Weight at which KL/2 is attained in Sigmoidal K model (GrowthModel_ConvEfficiency_8-25-13.xlsx) n.pulse=5, #Number of prey pulses occuring from years 26-50 pulse.range = c(2:6), #Range of multiplier values for when a pulse in production occurs. E.g. 2-10x production in pulse years N.prey.2.mean = 3e+10 , #Number of prey 2 (calc from Jung and Houde 2004 (MEPS) based on assumption of 1g/fish). Note that this number gets multiplied by the 1g to convert to biomass in the FxnR model. N.prey.dev = 0.4, #Standard deviation for the lognormal random error of prey populations relB.1 = 10 , #Number of prey 1 relative to prey 2; because treating each prey as 1g (for all prey), this is equivalent to relative biomass relB.3 = .1, #Number of prey 3 relative to prey 2; because treating each prey as 1g (for all prey), this is equivalent to relative biomass atck.max.1 =exp(-19.5), atck.rho.1= -1.1, atck.infl.1= 2, #Parameters for sigmoidal Attack rate for prey 1 (max value, rate value rho (neg for decr fxn), and inflection point (age at 50%max)) atck.max.2a=exp(-16.3), #Parameters for double-logistic attack rate for prey 2 (max.2a = scaling parameter for max value) atck.rho.2a=0.46, atck.rho.2b = 0.46, #Parameters for double-logistic attack rate for prey 2 (rho = rate parameter for ascending (2a) and descending (2b) limbs) atck.infl.2a=1, atck.infl.2b=4, #Parameters for double-logistic attack rate for prey 2 (infl = age of inflection point for ascending (2a) and descending (2b) limbs) (see Buchheister disseration Chapter 3) #atck.max.2 =exp(-16.7), atck.rho.2= 0.0, atck.infl.2= 2, #OLD Parameters for sigmoidal Attack rate for prey 2 (changed to double logistic) atck.max.3 =exp(-13.5), atck.rho.3= 1.9, atck.infl.3= 3.5, #Parameters for sigmoidal Attack rate for prey 3 (max value, rate value rho (neg for decr fxn), and inflection point (age at 50%max)) F.scenario = "highF", #F scenarios = "highF","lowF","decrF". P.scenario = 1, #1-no prey trends, "Pulse1", "Pulse2", "Pulse3", "Step.1", "Step.2", "Step.3" I.scenario="mig", #Migration scenario, either "mig", or "mix" graph=T, #Switch to turn output plots on or off. fxn=F #Leave this as FALSE. It is needed to allow the growthfxn() to be imbedded in the PredPreySim() code. It acts as a switch to turn on section of the code when running the growthfxn() ) { #For Troubleshooting: if(FALSE) { #Comment out trouble shooting parameters seed=500; M=0.25; Fmax=1.5; Fmin=0.31; decrF.time=12; F.sel = c(0,0.1,0.5,rep(1.0,ages-3)); R.dev.sd=0.4; alpha.SR = 3.4; K.SR = 27.4; beta.SR = 1; Gmax.a = rep(0.6,ages); CTO=22; CTM=35; CQ=2.5; CA=0.3; CB=-0.2; KL = 0.5; KR = -0.0014; KW = -600; n.pulse=5; pulse.range = c(2:6);N.prey.2.mean = 3e+10 ; N.prey.dev = 0.4; relB.1 = 10 ; relB.3 = .1; atck.max.1 =exp(-19.5); atck.rho.1= -1.1; atck.infl.1= 2; atck.max.2a=exp(-16.3); atck.rho.2a=0.46; atck.rho.2b = 0.46; atck.infl.2a=1; atck.infl.2b=4; atck.max.3 =exp(-13.5); atck.rho.3= 1.9; atck.infl.3= 3.5; F.scenario = "user.specified"; P.scenario = 1; I.scenario="mig"; graph=T; fxn=F } #################################### ### DEFINE PARAMETERS ### #################################### seed = seed #for setting the random seed number ## RECRUITMENT PARAMETERS (Can i get S-R parameters from Terceiro?... not clear in his Stock Assessment report) #Set seasonal life history parameters s.spawn <<- 4 #season of spawning (1=spring, 2=summer, 3=fall, 4=winter) (REALLY THIS SHOULD BE: (1=winter, 2=spring, 3=summer, 4=fall !! 9-3-13) s.recruit <<- 1 #season in which new individuals recruit to the population #Stock-recruitment parameters (see StockRecruitFxn_12-20-11.xlsx, Crecco 2008, and Rothschild et al. 2012) #Where R = alpha.SR * S / (1 + S/K.SR)^beta.SR (Rothschild et al. 2012) (UNITS: SSB in kmt; recrtuitment in 10^6) #alpha.P = 23.65; beta.P = 1/1488.3 #OLD values. Taken from Table 10 in Crecco 2008. (UNITS: SSB in mt, age-0 recruitment in thousands(?) of fish) alpha.SR = alpha.SR #Rothschild(UNITS: SSB in kmt (1E6 kg), age-0 recruitment in millions (1E6) of fish) K.SR = K.SR beta.SR = beta.SR #This parameter determines the shape of the SR function. beta.SR=1 is for Beverton Holt. R.dev.sd = R.dev.sd #the standard deviation for random error on recruitment. SD=0.4 is temporary (see Beddington and Cook, 1983 from Pat Lynch) set.seed(seed*1) R.deviates = exp(rnorm(years,mean=0,sd=R.dev.sd)) #time series of random deviates for recruitment. #Proportion of recruits that recruit to region r (theta). Note, sum of values must be 1.0. theta = matrix(nrow=regions, ncol=1, dimnames=list(paste("r",1:regions,sep=""),"theta")) #regions as rows if (I.scenario=="mig") theta[,] = c(0.9,0.1) #for migration scen. assume 90% recruitment to region 1. if (I.scenario=="mix") theta[,] = c(0.5,0.5) #for mixed scen. assume 50% recruitment to region 1 ## MORTALITY PARAMETERS #Mortality parameters for predator (Natural mortality=M; Fishing mortality = F) (UNITS = year**(-1)) #M = 0.25 #Value specified in function call. See Terceiro 2011 Assessment, p. 21 #Selectivity by Age sel.a = matrix(nrow=ages, ncol=1, dimnames=list(paste("age",0:(ages-1),sep=""),"sel.a")) #Fishing selectivity by age sel.a[,1] = F.sel #sel.a[,1] = c(0,0.1,0.5,rep(1.0,ages-3)) #See Terceiro 2011 Stock assessment, p 21); assume fully recruited from age 3 on. #F scenarios if(F.scenario == "highF") F.rate = seq(Fmax,Fmax,length=years) #F=2 was the maximum F obtained from 1982 to 2010 (Terceiro 2011 Table 56) if(F.scenario == "lowF") F.rate = seq(Fmin,Fmin,length=years) #F=0.255 is the target F as of Jan 2013. if(F.scenario == "decrF") F.rate = c(seq(Fmax,Fmax,length=24), seq(Fmax,Fmin,length=decrF.time), rep(Fmin,years-(24+decrF.time))) #Fmax for 25 years, F drops linearly to target F (Fmin) for x years, then holds at target for remaining years. this mimics the observed trend (Terceiro 2011) if(F.scenario == "dec2F") F.rate = c(seq(Fmax,Fmax,length=24), seq(Fmax,Fmin,length=decrF.time/2), rep(Fmin,years-(24+decrF.time/2))) #Fmax for 25 years, F drops linearly to target F (Fmin) for x/2 years, then holds at target for remaining years. this mimics the observed trend (Terceiro 2011) if(F.scenario == "exp.F") {exp.decl = Fmax*exp(1/(decrF.time-1)*log(Fmin/Fmax)*c(0:(decrF.time-1))) #This equation is for an exponential decline in F from Fmax to Fmin over the course of decrF.time years. Note the "-1" is F.rate = c(seq(Fmax,Fmax,length=24), exp.decl, rep(Fmin,years-(24+decrF.time))) } #Fmax for 25 years, F drops exponentially to target F (Fmin) for 20 years, then holds at target for remaining years. if(F.scenario == "ass.F") F.rate = c(seq(Fmax,Fmax,length=10), SRdata$F, rep(Fmin,years-39)) #This uses the estimated annual F's from the stock assessment; with earlier years forced to be Fmax, and later years to be Fmin (Terceiro 2011) if(F.scenario == "knifF") F.rate = c(seq(Fmax,Fmax,length=25), rep(Fmin,years-25)) #Knife-edge: Fmax for 25 years, then holds at target for remaining years. if(F.scenario == "mor.F") F.rate = c(seq(Fmax,Fmax,length=24), seq(0,0,length=decrF.time), rep(0,years-(24+decrF.time))) #Fmax for 25 years, then moratorium imposed #Sub scenarios for "decrF", where the decline in F occurs over a shorter or longer period if(F.scenario == "decrF.24") F.rate = c(seq(Fmax,Fmax,length=24), seq(Fmax,Fmin,length=24), rep(Fmin,years-48)) #Fmax for 25 years, F drops linearly to target F (Fmin) for 24 years, then holds at target for remaining years. if(F.scenario == "decrF.22") F.rate = c(seq(Fmax,Fmax,length=24), seq(Fmax,Fmin,length=22), rep(Fmin,years-46)) #Fmax for 25 years, F drops linearly to target F (Fmin) for 22 years, then holds at target for remaining years. if(F.scenario == "decrF.21") F.rate = c(seq(Fmax,Fmax,length=24), seq(Fmax,Fmin,length=21), rep(Fmin,years-45)) #Fmax for 25 years, F drops linearly to target F (Fmin) for 21 years, then holds at target for remaining years. if(F.scenario == "decrF.19") F.rate = c(seq(Fmax,Fmax,length=24), seq(Fmax,Fmin,length=19), rep(Fmin,years-43)) #Fmax for 25 years, F drops linearly to target F (Fmin) for 19 years, then holds at target for remaining years. if(F.scenario == "decrF.18") F.rate = c(seq(Fmax,Fmax,length=24), seq(Fmax,Fmin,length=18), rep(Fmin,years-42)) #Fmax for 25 years, F drops linearly to target F (Fmin) for 18 years, then holds at target for remaining years. if(F.scenario == "decrF.16") F.rate = c(seq(Fmax,Fmax,length=24), seq(Fmax,Fmin,length=16), rep(Fmin,years-40)) #Fmax for 25 years, F drops linearly to target F (Fmin) for 16 years, then holds at target for remaining years. if(F.scenario == "decrF.14") F.rate = c(seq(Fmax,Fmax,length=24), seq(Fmax,Fmin,length=14), rep(Fmin,years-38)) #Fmax for 25 years, F drops linearly to target F (Fmin) for 14 years, then holds at target for remaining years. if(F.scenario == "decrF.12") F.rate = c(seq(Fmax,Fmax,length=24), seq(Fmax,Fmin,length=12), rep(Fmin,years-36)) #Fmax for 25 years, F drops linearly to target F (Fmin) for 12 years, then holds at target for remaining years. if(F.scenario == "decrF.10") F.rate = c(seq(Fmax,Fmax,length=24), seq(Fmax,Fmin,length=10), rep(Fmin,years-34)) #Fmax for 25 years, F drops linearly to target F (Fmin) for 10 years, then holds at target for remaining years. if(F.scenario == "decrF.08") F.rate = c(seq(Fmax,Fmax,length=24), seq(Fmax,Fmin,length=08), rep(Fmin,years-32)) #Fmax for 25 years, F drops linearly to target F (Fmin) for 08 years, then holds at target for remaining years. #Calculate F at each age F.t = matrix(nrow=years, ncol=1, dimnames=list(paste("yr",1:years,sep=""),"F.t")) #F trend over time F.t[,1] = F.rate F.a = matrix(nrow=years, ncol=ages, dimnames=list(paste("yr",1:years,sep=""),paste("age",0:(ages-1),sep=""))) #age-specific F over time for (y in 1:years) (F.a[y,] = sel.a * F.t[y] ) ## IMMIGRATION PARAMETERS (TEMPORARY VALUES!!! Perhaps make this AGE DEPENDENT?) #Immigration coefficients. [phi.in = Proportion of fish from region n moving into region r at the end of season s] #Movement occurs at the end of the seasonal time step (after the fish die and grow) #Sources: Packer et al. 1999, figure 3, p45 for seasonal maps; Terceiro 2002; Bonzek et al. 2010 (ChesMMAP Reports). ## SPECIFY MIGRATION SCENARIOS (NOTE: IMMIGRATION (PHI.IN) IS AGE DEPENDENT. phi.in = array(NA,dim=c(ages,regions,regions,seasons),dimnames=list(paste("age",0:(ages-1),sep=""),paste("r",1:regions,sep=""),paste("n",1:regions,sep=""),paste("s",1:seasons,sep=""))) if (I.scenario=="mig_old"){ #OLD/ORIGINAL MIGRATION SCENARIO phi.in[,,,1] = c(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, #values for n1 into r1 for each age (season 1) 0.1, 0.1, 0.1, 0.2, 0.2, 0.2, 0.2, 0.2, #values for n1 into r2 for each age 0.0, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, #values for n2 into r1 for each age 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0) #values for n2 into r2 for each age phi.in[,,,2] = c(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, #values for n1 into r1 for each age (season 2) 0.0, 0.0, 0.0, 0.1, 0.1, 0.1, 0.1, 0.1, #values for n1 into r2 for each age 0.7, 0.7, 0.7, 0.7, 0.7, 0.7, 0.7, 0.7, #values for n2 into r1 for each age 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0) #values for n2 into r2 for each age phi.in[,,,3] = c(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, #values for n1 into r1 for each age (season 3) 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, #values for n1 into r2 for each age 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, #values for n2 into r1 for each age 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0) #values for n2 into r2 for each age phi.in[,,,4] = c(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, #values for n1 into r1 for each age (season 4) 0.4, 0.8, 0.8, 0.9, 0.9, 0.9, 0.95, 0.95, #values for n1 into r2 for each age #see Packer et al. 1999, p6&7; YOY fish (and adults can overwinter in the bay 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, #values for n2 into r1 for each age 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0) } #values for n2 into r2 for each age if (I.scenario=="mig"){ #NEW MIGRATION SCENARIO - 9-4-13 phi.in[,,,1] = c(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, #values for n1 into r1 for each age (season 1) 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, #values for n1 into r2 for each age 0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3, #values for n2 into r1 for each age #A little bit of a discrepancy between Packer et al. 1999 and Terceiro 2002 (Terceiro showed lower catches offshore); so, i split the diff a bit. 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0) #values for n2 into r2 for each age phi.in[,,,2] = c(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, #values for n1 into r1 for each age (season 2) 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, #values for n1 into r2 for each age 0.95, 0.95, 0.95, 0.95, 0.95, 0.95, 0.95, 0.95, #values for n2 into r1 for each age #Very few flounder caught offshore (Packer 1999; Terceiro 2002) 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0) #values for n2 into r2 for each age phi.in[,,,3] = c(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, #values for n1 into r1 for each age (season 3) 0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4, #values for n1 into r2 for each age #Emmigration to offshore begins, as supported with the higher catch rates in the NMFS surveys at this time. 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, #values for n2 into r1 for each age #Assume no more fish coming into estuaries/nearshore at this time... instead they are starting to leave. 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0) #values for n2 into r2 for each age phi.in[,,,4] = c(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, #values for n1 into r1 for each age (season 4) 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, #values for n1 into r2 for each age #March min trawlable number in ChesBay in March ~10% of peak value (Bonzek et al. 2010). 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, #values for n2 into r1 for each age 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0) } #values for n2 into r2 for each age if (I.scenario=="mix"){ #FULLY MIXED SCENARIO phi.in[,,,1] = c(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, #values for n1 into r1 for each age (season 1) 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, #values for n1 into r2 for each age 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, #values for n2 into r1 for each age 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0) #values for n2 into r2 for each age phi.in[,,,2] = c(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, #values for n1 into r1 for each age (season 2) 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, #values for n1 into r2 for each age 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, #values for n2 into r1 for each age 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0) #values for n2 into r2 for each age phi.in[,,,3] = c(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, #values for n1 into r1 for each age (season 3) 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, #values for n1 into r2 for each age 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, #values for n2 into r1 for each age 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0) #values for n2 into r2 for each age phi.in[,,,4] = c(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, #values for n1 into r1 for each age (season 4) 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, #values for n1 into r2 for each age 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, #values for n2 into r1 for each age 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0) } #values for n2 into r2 for each age ## INITIAL POPULATION SIZE #Initial population abundance for year 1 (See Terceiro 2011, Table 57 for population size estimates) #Note: P[years, seasons, regions, ages] P0 = array(NA, dim=c(years,seasons,regions,ages),dimnames=list(1:years,paste("s",1:seasons,sep=""),paste("reg",1:regions,sep=""),paste("age",0:(ages-1),sep=""))) for (s in 1:seasons) { #Total population abundance estimate for 1989, divided equally among the regions, and constant across seasons for first year (from Terceiro 2011, Table 57) for (r in 1:regions) { P0[1,s,r,] = c(28000,9000,11000,2000,300,40,10,rep(0,ages-7))*1000/regions }} ## TEMPERATURE INPUTS AND EFFECTS ON CONSUMPTION (by season and region) #Mean seasonal temperatures by region to modify prey consumption in the growth submodel Temp.seas = array(NA,dim=c(regions, seasons),dimnames=list(paste("reg",1:regions,sep=""), paste("s",1:4,sep="") )) #UNITS: degrees C Temp.seas[1,] = c(4.4, 13.6, 24.4, 14.8) #Mean seasonal bottom temps from CB4.2C (see Temp_ForcingFunction_ChesBay_xxx.xlsx Temp.seas[2,] = c(7.6, 10.8, 13.5, 15.2) #Based on bottom Temp data from NOAA cruises (see NOAA_Temp Data_offshore_9-5-13.xlsx) #Temp modifier function (f.Temp) Z=log(CQ)*(CTM-CTO) Y=log(CQ)*(CTM-CTO+2) V=(CTM-Temp.seas)/(CTM-CTO) X=(Z^2*(1+(1+40/Y)^0.5)^2)/400 f.Temp=V^X*(exp(X*(1-V))) ## GROWTH PARAMETERS (See also Table 27 of Terceiro 2011 for mean weight at age of summer flounder catch) #Weight at age (w) (UNITS = KG) (Is this the weight at the beginning or end of the season???!) #VonBert relationship: L = 60.79094*(1-exp(-0.29949*age-0.27934) (see vonBert_ab.sas; sexes combined, L in cm, age is fractional (jan1 birthday)) L.a = array(NA,dim=c(seasons,ages),dimnames=list(paste("s",1:seasons,sep=""),paste("age",0:(ages-1),sep=""))) for (s in 1:seasons) { L.a[s,] = 60.79094*(1-exp(-0.29949*(c(0:(ages-1))+s/4)-0.27934) ) } #Von Bertalanffy equation. (note values are calc for end of each season.) #L-W flounder relationship: W = 0.004324*L^3.232744 (n=4981,R2=0.99; ChesMMAP data; W in g, L in cm. See "LW_Parameters_10-25-11.xls") w = array(NA,dim=c(seasons,ages),dimnames=list(paste("s",1:seasons,sep=""),paste("age",0:(ages-1),sep=""))) for (s in 1:seasons) { w[s,] = (0.004324*L.a[s,]^3.232744)/1000 } #L-W equation. (note values are calc for end of each season.) UNITS: KG #Starting weight of age-0 fish in season 1 (w0) w0 = w[1,1] #Maturity at age (m) - see Terceiro 2011 Assessment, p.18. Note, age-0 omitted! m = array(NA,dim=c(seasons,ages),dimnames=list(paste("s",1:seasons,sep=""),paste("age",0:(ages-1),sep=""))) m[,1] = 0.38 #age-0 % mature in Fall. m[,2] = 0.72 #age-1 % mature in Fall. m[,3] = 0.90 #age-2 % mature in Fall. m[,4:ages] = 1 #age-3+ % mature in Fall. ## INDIVIDUAL GROWTH SUBMODEL (based on conversion efficiency, K - proportion of consumed biomass converted to predator biomass; size-dependent) #Growth function ("G.func") defined as G = K*C, where G and C (UNITS = g/season). G.func = function(WT, C, KL.x = KL, KR.x = KR, KW.x = KW) { K = KL.x / (1 + exp(-KR.x*(WT*1000 - KW.x))) #K modeled as sigmoidal function. Need to convert weight to g for K equation. See GrowthModel_ConvEfficiency_8-25-13.xlsx. G = K*C return(G) #(UNITS = g/seas) } #End G.func ## FUNCTIONAL RESPONSE (ATTACK RATE) PARAMETERS #These values have been optimized to have results match weight-at-age and %diet data (see MS Functional Response_8-27-13.xlsx and Flounder_Diets by age_1-29-13.xlsx) #Values optimized using N.prey.2.mean=3e+10, relB.1=10, relB.3=0.1 #Attack rate parameters for logistic function (for prey 1 and 3) is defined as atck = atck.max / (1+exp(-atck.rho(age - atck.infl))); atck.max = maximum attack rate, atck.rho = rate to reach max, atck.infl = age at 50% attack rate. #Attack rate for double logistic function (prey 2) is defined as: atck = atck.max.2a*(1/(1+exp(-atck.rho.2a*(age-atck.infl.2a))))*(1-(1/(1+exp(-atck.rho.2b*(age-atck.infl.2b))))); see definitions in function call. #To optimize the values, i needed to use this embedded growth function: ################################################## ### GROWTH FUNCTION ################################################## #This function is necessary for optimizing the attack rate parameters... #To run the Full PredPreySim, comment out the growth fxn definition, and wrap this all up with the full PredPreySim function #Otherwise, use this function and the attackrates() function later in the code to test, evaluate and optimize the functional response parameters if (FALSE) { growthfxn <<- function( #Parameters as above, except unnecessary parms commented out. seed = 500, #M=0.25, #Instantaneous natural mortality rate (per year) #F.rate = seq(2,2,length=years), #annual fishing mortality rate over the time series #F.sel = c(0,0.1,0.5,rep(1.0,ages-3)), #Fishing selectivity by age #R.dev.sd=0.4, #Standard deviation for the lognormal random recruitment error #alpha.SR = 3.4, #alpha parameter for Shepherd S-R function (Rothschild et al. 2012) #K.SR = 27.4, #K parameter for Shepherd S-R function (Rothschild et al. 2012) #beta.SR = 1, #beta shape parameter for Shepherd S-R function. beta.SR=1 is for Beverton Holt. Gmax.a = rep(0.6,ages), #Age-specific maximum growth (kg/yr) possible for species n.pulse=5, #Number of prey pulses pulse.range = c(2:6), #Range of multiplier values for when a pulse in production occurs. E.g. 2-10x production in pulse years N.prey.2.mean = 3e+10 , #Number of prey 2 N.prey.dev = 0.4, #Standard deviation for the lognormal random error of prey populations relB.1 = 10 , #Number of prey 1 relative to prey 2; because treating each prey as 1g (for all prey), this is equivalent to relative biomass relB.3 = .1, #Number of prey 3 relative to prey 2; because treating each prey as 1g (for all prey), this is equivalent to relative biomass atck.max.1 =exp(-19.5), atck.rho.1= -1.1, atck.infl.1= 2, #Parameters for sigmoidal Attack rate for prey 1 (max value, rate value rho (neg for decr fxn), and inflection point (age at 50%max)) atck.max.2a=exp(-16.3), #Parameters for double-logistic attack rate for prey 2 (max.2a = scaling parameter for max value) atck.rho.2a=0.46, atck.rho.2b = 0.46, #Parameters for double-logistic attack rate for prey 2 (rho = rate parameter for ascending (2a) and descending (2b) limbs) atck.infl.2a=1, atck.infl.2b=4, #Parameters for double-logistic attack rate for prey 2 (infl = age of inflection point for ascending (2a) and descending (2b) limbs) #atck.max.2 =exp(-16.7), atck.rho.2= 0.0, atck.infl.2= 2, #OLD Parameters for sigmoidal Attack rate for prey 2 (changed to double logistic) atck.max.3 =exp(-13.5), atck.rho.3= 1.9, atck.infl.3= 3.5, #Parameters for sigmoidal Attack rate for prey 3 (max value, rate value rho (neg for decr fxn), and inflection point (age at 50%max)) F.scenario = "user.specified", #F scenarios = "highF","lowF","decrF". Use "user.specified" to use the F.rate value identified in the function call. P.scenario = 1, #1-no prey trends, "Pulse1", "Pulse2", "Pulse3" I.scenario="mig", #Migration scenario, either "mig", or "mix" graph=T, #Switch to turn output plots on or off. #Parameters specific or crucial to growthfxn fxn=T, #Switch indicating that the growthfxn() is being run, to calculate w.pred (without accounting for migration) optim.w=T, #Use this to output the RSS for the weight at age relationship (needed if optimizing FxnR parameters for this) optim.diet=T, #Use this to output the RSS for the percent diet relationships (needed if optimizing FxnR parameters for this) optim.both=T, #Use this to output the RSS for the w@age AND percent diet relationships summed(needed if optimizing FxnR parameters for this) wt4optim = 3, #this is a weighting factor indicating how much more to weight the weight at age RSS over the diet % optim.weighted = T #This outputs a weighted sum of both RSS according to the wt4optim parameter ) { } #Use this if fxn=F, to close the growthfxn() code } #COMMENTING OUT THE GROWTH FUNCTION DEFINITION ## FUNCTIONAL RESPONSE PARAMETERS (CONT'D) #Attack rate parameters for logistic function (where atck = atck.max / (1+exp(-atck.rho(age - atck.infl))); atck.max = maximum attack rate, atck.rho = rate to reach max, atck.infl = age at 50% attack rate. #NEW (8/27/13) Attack rates (atck) (UNITS: numbers of prey / ind. predator / time) Assume atck changes logistically with age (but double logistic for prey 2), age is fractional over the year atck = array(NA,dim=c(seasons,ages,preys),dimnames=list(paste("s",1:seasons,sep=""),paste("age",0:(ages-1),sep=""),paste("prey",1:preys,sep=""))) atck[,,1] = atck.max.1 /(1 + exp(-atck.rho.1*(seq(from=0,by=0.25,length=ages*seasons)-atck.infl.1))) #Assume attack rate for prey 1 is decreasing sigmoidally with age. atck[,,2] = atck.max.2a*(1/(1+exp(-atck.rho.2a*(seq(from=0,by=0.25,length=ages*seasons)-atck.infl.2a)))) * #Assume attack rate follows a double logistic (dome-shaped) (1-(1/(1+exp(-atck.rho.2b*(seq(from=0,by=0.25,length=ages*seasons)-atck.infl.2b))))) atck[,,3] = atck.max.3 /(1 + exp(-atck.rho.3*(seq(from=0,by=0.25,length=ages*seasons)-atck.infl.3))) #Assume attack rate for prey 3 is increasing sigmoidally with age. if(FALSE) {#Comment out plotting code and OLD code #Plot attack rates ylim.atck = c(0,max(atck[,,3])) #ylim.atck = NULL plot(seq(from=0,by=0.25,length=ages*seasons),atck[,,1], ylim=ylim.atck, type="l", col="black") par(new=T) plot(seq(from=0,by=0.25,length=ages*seasons),atck[,,2], ylim=ylim.atck, type="l", col="red") par(new=T) plot(seq(from=0,by=0.25,length=ages*seasons),atck[,,3], ylim=ylim.atck, type="l", col="blue") } #OLD (1/29/13) Attack Rates (atck) (UNITS: numbers of prey / ind. predator / time). Assume atck changes logistically with age, age is fractional over the year #atck = array(NA,dim=c(seasons,ages,preys),dimnames=list(paste("s",1:seasons,sep=""),paste("age",0:(ages-1),sep=""),paste("prey",1:preys,sep=""))) #atck[,,1] = atck.max.1 /(1 + exp(-atck.rho.1*(seq(from=0,by=0.25,length=ages*seasons)-atck.infl.1))) #Assume attack rate for prey 1 is decreasing with age. #atck[,,2] = atck.max.2 /(1 + exp(-atck.rho.2*(seq(from=0,by=0.25,length=ages*seasons)-atck.infl.2))) #Assume attack rate for prey 2 is increasing with age. #if(preys==3) atck[,,3] = atck.max.3 /(1 + exp(-atck.rho.3*(seq(from=0,by=0.25,length=ages*seasons)-atck.infl.3))) #Assume attack rate for prey 3 is increasing with age. #Max Consumption (Cmax) (UNITS: prey mass(KG) / ind. predator / season) #Cmax (in g/g/d) = CA * weight ^ CB (following allometric model from Wisconsin model manual - Hanson et al. 1997) #CA=CA #CA and CB values mimic observed trends max stomach fullness and are comparable to values for other temperate species (Hartman and Brandt 1995a) #CB=CB Cmax = array(as.vector(c(NA,w[,] * (CA*(w[,]*1000)^CB) *91)), #Need to shift Cmax values by 1 season, because Cmax is calculated from season s-1 (because weights are measured at the end of a season). 9-3-13 dim=c(seasons,ages),dimnames=list(paste("s",1:seasons,sep=""),paste("age",0:(ages-1),sep=""))) #OLD: Cmax = array(NA,dim=c(seasons,ages),dimnames=list(paste("s",1:seasons,sep=""),paste("age",0:(ages-1),sep=""))) #OLD: Cmax[,] = w[,] * (CA*(w[,]*1000)^CB) *91 #1000 is for converting kg to g, and 91 is for expanding consumption from day to season. #OLDer: Cmax[,] = w[,] * Cmax.perc #Assume Cmax = 5% BW/d * 91 days for all fish (see Wisconsin model manual) if(FALSE){ #Comment out because used with old method of calculating growth #OLD: Max Growth (Gmax) (UNITS: kg/season) (See "MS Functional Response_11-29-11.xlsx"; tuned to match weight at age data) Gmax = array(NA,dim=c(seasons,ages),dimnames=list(paste("s",1:seasons,sep=""),paste("age",0:(ages-1),sep=""))) for (s in 1:4) { Gmax[s,] = Gmax.a/4 } #Gmax.a (age-specific max growth) is input into the function #Gmax[s,] = c(0.6,0.6,0.5,0.5,0.4,0.4,0.3,0.2)/4 } #New values (1/29/13). See "MS Functional Response_11-29-11.xlsx". Note that the fastest growth of my age-0 fish for the masters project was 300g in 1/2 year, so 0.6 kg/yr is good. } ## STORE PARAMETERS AND ARRAYS AS NEEDED (for use outside the function) L.a<<-L.a w <<- w w0 <<- w0 m <<- m phi.in <<- phi.in atck<<-atck Cmax <<- Cmax #################################### ### DEFINE PREY SCENARIOS ### #################################### ## SPECIFY PREY PARAMETERS (FOR MEAN AND VARIATION), CREATE EMPTY ARRAY #Define mean prey abundance levels (UNITS = number of individuals; or Biomass/1000 because treating all individual prey weights as 1g): N.prey.2.mean = N.prey.2.mean #3e+10 individuals... calc from ~30*10^3 metric tons (see Jung and Houde 2004 (MEPS)), assuming each fish is 1g N.prey.1.mean = N.prey.2.mean * relB.1 #biomass of prey 1 is a scalar of prey 2 biomass. N.prey.3.mean = N.prey.2.mean * relB.3 #biomass of prey 3 is a scalar of prey 2 biomass. N.prey.1.dev = N.prey.dev #Arbitrary parameter governing the lognormal deviates of prey abundance N.prey.2.dev = N.prey.dev N.prey.3.dev = N.prey.dev N.prey.mean.all = c(N.prey.1.mean, N.prey.2.mean, N.prey.3.mean) #This is being used to test the alternative method for calculating the prey pulses. #Avg. Prey Weights (w) (UNITS: mass-kg) w.prey = matrix(nrow=preys, ncol=1, dimnames=list(paste("prey",1:preys,sep=""),"prey.wt")) w.prey[,1] = rep(0.001,preys) #Functional response ultimately dependent on biomasses of prey. So, standardizing prey to 1g and then have biomasses controlled above by prey numbers #OLD: w.prey[,1] = c(0.0002,0.003) #Mysid weight (.2g) estimated from masters isotope analyses (see Mysid Weights_6-11-12.xls); Anchovy weight (3g) is mean from ChesMMAP catches; #Create empty prey abundance table (N.prey) for years, seasons, regions, preys N.prey = array(NA, dim=c(years,seasons, regions, preys),dimnames=list(1:years,paste("s",1:seasons,sep=""),paste("reg",1:regions,sep=""),paste("prey",1:preys,sep=""))) #Matrix for prey abundance time series ##DEFINE PREY SCENARIOS (N.prey[years,seas,regs,preys]) #Prey Scenario 1 (constant prey abundances); This array will be re-written or modified based on the other scenarios #if (P.scenario == 1) { #Transpose/permute the array to allow N.prey values to be entered by rows instead of columns. Otherwise, runs with years=60 or 80 will result in different output for the first 60 years, when they should be identical. N.prey.perm = aperm(N.prey, c(2:4,1)); dim(N.prey.perm) #Fill in array set.seed(seed*2) N.prey.perm[,,1,] = N.prey.1.mean * exp(rnorm(years*seasons*regions,mean=0,sd=N.prey.1.dev)) #temporarily assume that the trend is identical for regions and same value across seasons. set.seed(seed*3) N.prey.perm[,,2,] = N.prey.2.mean * exp(rnorm(years*seasons*regions,mean=0,sd=N.prey.2.dev)) #temporarily assume that the trend is identical for regions and same value across seasons. set.seed(seed*3.5) if(preys==3) N.prey.perm[,,3,] = N.prey.3.mean * exp(rnorm(years*seasons*regions,mean=0,sd=N.prey.3.dev)) #temporarily assume that the trend is identical for regions and same value across seasons. #Transpose/permute the array back to the original orientation N.prey = aperm(N.prey.perm, c(4,1:3)) N.prey.base = N.prey #NEW 3/13/13. Store the base prey biomasses #plot(seq(1,by=0.25,length=years*seasons),t(N.prey.base[,,1,preyX])/1E10, type="l") #abline(h=c(N.prey.1.mean, N.prey.2.mean, N.prey.3.mean)/1E10,col=c("black","red","blue"),lty=2) #} #Prey Scenario 2 (decline in prey 1, increase in prey 2) #temporarily assume trend is identical for regions if (P.scenario == 2) { prey.multplyr = 2 #starting prey value will be this value times the prey mean, and ending value is prey mean divided by this multiplier set.seed(seed*4) N.prey[,,,1] = matrix(data=seq(N.prey.1.mean*prey.multplyr,N.prey.1.mean/prey.multplyr,length=years*seasons)* exp(rnorm(years*seasons,mean=0,sd=N.prey.1.dev)),nrow=years,byrow=T) #also incorporating error into the prey trajectories set.seed(seed*5) N.prey[,,,2] = matrix(data=seq(N.prey.2.mean*prey.multplyr,N.prey.2.mean/prey.multplyr,length=years*seasons)* exp(rnorm(years*seasons,mean=0,sd=N.prey.2.dev)),nrow=years,byrow=T) } #plot(seq(1,by=1,length=years*seasons),t(N.prey[,,1,1])) #Prey Pulse Scenarios (pulses of prey X production, anywhere from 2 to 10x normal production ) if (P.scenario %in% c("Pulse1","Pulse2","Pulse3")) { preyX = as.numeric(substr(P.scenario,6,6)) #Store which prey # is experiencing the pulse set.seed(seed*6) pulse.years = sample(26:50,n.pulse) #Pulses can only occur from year 26-50 set.seed(seed*6) N.prey[pulse.years,,1,preyX] = N.prey[pulse.years,,1,preyX] * runif(n.pulse, min(pulse.range),max(pulse.range)) #Prey pulses are anywhere from 2 to 10x normal production #N.prey[pulse.years,,1,preyX] = N.prey[pulse.years,,1,preyX] * sample(pulse.range,n.pulse) #Prey pulses are anywhere from 2 to 10x normal production (integers only; OLD WAY) #plot(seq(1,by=0.25,length=years*seasons),t(N.prey[,,1,preyX])/1E11, type="l") #abline(h=c(N.prey.1.mean, N.prey.2.mean, N.prey.3.mean)/1E11,col=c("black","red","blue"),lty=2) } #Alternative calculation for prey pulses (Pulses f prey X production, adding 2-10x mean production) if (P.scenario %in% c("PulsB1","PulsB2","PulsB3")) { preyX = as.numeric(substr(P.scenario,6,6)) #Store which prey # is experiencing the pulse set.seed(seed*6) pulse.years = sample(26:50,n.pulse) #Pulses can only occur from year 26-50 set.seed(seed*6) N.prey[pulse.years,,1,preyX] = N.prey[pulse.years,,1,preyX] + sample(pulse.range,n.pulse)*N.prey.mean.all[preyX] #Prey pulses are anywhere from 2 to 10x normal production #plot(seq(1,by=0.25,length=years*seasons),t(N.prey[,,1,preyX])/1E11, type="l") #abline(h=c(N.prey.1.mean, N.prey.2.mean, N.prey.3.mean)/1E11,col=c("black","red","blue"),lty=2) } #Prey "Step" Scenarios (step increase in prey X production from years 26-50 (in Region 1 only) to match the total production in the pulse scenario) NEW: 3/13/13 if (P.scenario %in% c("Step.1","Step.2","Step.3")) { #note: use of "." is to maintain same number of characters as "Pulse" preyX = as.numeric(substr(P.scenario,6,6)) #Store which prey # is experiencing the pulse #Calc total production added by pulses (analogous to area under the curve) set.seed(seed*6) pulse.years = sample(26:50,n.pulse) #Pulses can only occur from year 26-50 N.prey.pulse = N.prey.base #Create a copy of N.prey.base to store the new pulsed production values set.seed(seed*6) N.prey.pulse[pulse.years,,1,preyX] = N.prey.base[pulse.years,,1,preyX] * runif(n.pulse,min(pulse.range),max(pulse.range)) #Prey pulses are anywhere from 2 to 10x normal production #N.prey.pulse[pulse.years,,1,preyX] = N.prey.base[pulse.years,,1,preyX] * sample(pulse.range,n.pulse) #Prey pulses are anywhere from 2 to 10x normal production (integer only; OLD WAY) TotProdIncr = sum(N.prey.pulse[pulse.years,,1,preyX] - N.prey.base[pulse.years,,1,preyX]) #Re-allocate the production added by pulses across all 25 years so that there is a step increase in the base prey productivity that matches the productivity added in the pulse scenario. N.prey[26:50,,1,preyX] = N.prey.base[26:50,,1,preyX] + TotProdIncr/(25*4) #MODIFY THE TotProdIncr DENOMINATOR IF NEEDED!! (currently 25 yrs * 4 seasons) ProdStepIncr = (mean(N.prey[26:50,,1,preyX]) - mean(N.prey.base[26:50,,1,preyX]) ) / mean(N.prey.base[26:50,,1,preyX]) * 100 # %increase in production during "step" phase. (Value of 100 means a doubling of the base productivity (100% increase)) #Plot the prey time series #plot(seq(1,by=1,length=years*seasons),t(N.prey[,,1,1]), type="l", ylim=c(0,5E12)) #par(new=T) #plot(seq(1,by=1,length=years*seasons),t(N.prey.step[,,1,1]), type="l", col="red", ylim=c(0,5E12)) #plot(seq(1,by=1,length=years*seasons),t(N.prey[,,1,3]), type="l") #, ylim=c(0,5E12) #par(new=T) #plot(seq(1,by=1,length=years*seasons),t(N.prey[,,1,3]), type="l", col="red") #, ylim=c(0,5E12) } #Store data for use outside function N.prey <<- N.prey n.pulse<<- n.pulse N.prey.1.mean <<- N.prey.1.mean N.prey.2.mean <<- N.prey.2.mean N.prey.3.mean <<- N.prey.3.mean if (P.scenario %in% c("Pulse1","Pulse2","Pulse3","Step.1","Step.2","Step.3")) { pulse.years<<-pulse.years } if (P.scenario %in% c("Step.1","Step.2","Step.3")) { ProdStepIncr<<-ProdStepIncr } ################################################## ### RUN FUNCTIONAL RESPONSE AND GROWTH LOOPS ### ################################################## ## CREATE EMPTY MATRICES FOR DATA STORAGE #Matrices for functional response calculations sum.a.N.w = array(NA,dim=c(years,seasons, regions, ages),dimnames=list(1:years,paste("s",1:seasons,sep=""),paste("reg",1:regions,sep=""),paste("age",0:(ages-1),sep=""))) #Table of sum.a.N.w (across prey for all combinations of ages and years). This is in the denominator of the functional respones equation temp=matrix(NA, nrow=1, ncol=preys, dimnames=list("a.N.w",paste("prey",1:preys,sep=""))) #Temp table necessary for calculation FxnR = array(NA,dim=c(years,seasons,regions,ages,preys),dimnames=list(1:years,paste("s",1:seasons,sep=""),paste("reg",1:regions,sep=""),paste("age",0:(ages-1),sep=""),paste("prey",1:preys,sep=""))) #Table of sum.a.N.w (across prey for all combinations of ages and years). This is in the denominator of the functional respones equation ## CALCULATE PER CAPITA CONSUMPTION OF PREY BASED ON THEIR FORCED ABUNDANCES #Loop for calculating Functional response for (y in 1:years) { for (s in 1:seasons) { for (ai in 1:ages) { for (r in 1:regions) { for (j in 1:preys) { for (p in 1:preys) { #calculate a.N.w for each prey, and then sum it for use in the FxnR calculation denominator temp[1,p] = atck[s,ai,p]*N.prey[y,s,r,p]*w.prey[p] } sum.a.N.w[y,s,r,ai] = sum(temp) #Store sum(a.N.w) across prey and easier reference in the code #Populate the FxnR matrix. This gives the biomass of prey consumed per individ. predator by year, season, region, and age. FxnR[y,s,r,ai,j] = (Cmax[s,ai]*atck[s,ai,j]*N.prey[y,s,r,j]*w.prey[j]) / (Cmax[s,ai]+sum.a.N.w[y,s,r,ai]) }}}}} #FxnR[1,1,1,,]; N.prey[1,1,1,]; sum.a.N.w[1,1,1,] #Checking the matrices... #Store arrays for use outside of a function FxnR <<- FxnR sum.a.N.w <<- sum.a.N.w ## 1/29/13 - COMMENTING OUT THIS SECTION. Need to move it to population section below to have pred weights account for different weights of migrating fish. (Ie calculate weighted means of fish) if(fxn==T) { #ONLY INCLUDE THIS w.pred SECTION IF RUNNING THE growthfxn() CODE, NOT IF RUNNING PredPreySim() #Calculate weight-at-age based on the estimated consumption w.pred <<- array(NA,dim=c(years,seasons,regions,ages),dimnames=list(1:years,paste("s",1:seasons,sep=""),paste("reg",1:regions,sep=""),paste("age",0:(ages-1),sep=""))) for (r in 1:regions) (w.pred[1,,r,] = w ) #Set first year values to previously defined start values. for (y in 2:years) { for (s in 1:seasons) { for (ai in 1:ages) { for (r in 1:regions) { #Calculate weight at age for age-0 fish: if(ss.recruit && ai==1) ( w.pred[y,s,r,ai] = w.pred[y,s-1,r,ai] + G.func(WT=w.pred[y,s-1,r,ai],C=sum(FxnR[y,s,r,ai,]))*f.Temp[r,s] ) #weight of fish increases by Consumption * Conversion efficiency (K) * f.Temp #OLD: if(s>s.recruit && ai==1) ( w.pred[y,s,r,ai] = w.pred[y,s-1,r,ai] + Gmax[s,ai]*sum(FxnR[y,s,r,ai,])/Cmax[s,ai] ) #weight of fish increases by Gmax * (Fraction of max consumption) (following Wilberg's code) #Calculate weight at age for age>1 fish: (NOTE: i am not doing anything special for the plus group at this point. if(ai>1) { if(s==1) { w.pred[y,s,r,ai] = w.pred[y-1,seasons,r,ai-1] + G.func(WT=w.pred[y-1,seasons,r,ai-1],C=sum(FxnR[y,s,r,ai,]))*f.Temp[r,s] } #updated on 8/27/13 #{ w.pred[y,s,r,ai] = w.pred[y-1,seasons,r,ai-1] + Gmax[s,ai]*sum(FxnR[y,s,r,ai,])/Cmax[s,ai] } #OLD VERSION (1/29/13) #{ w.pred[y,s,r,ai] = w.pred[y-1,seasons,r,ai-1] + sum(FxnR[y,s,r,ai,]*mean(ce)) } #OLD VERSION (6/25/12): weight of fish increases by consumption * conversion efficiency if(s>1) { w.pred[y,s,r,ai] = w.pred[y,s-1,r,ai] + G.func(WT=w.pred[y,s-1,r,ai],C=sum(FxnR[y,s,r,ai,]))*f.Temp[r,s] } } #updated on 8/27/13 #{ w.pred[y,s,r,ai] = w.pred[y,s-1,r,ai] + Gmax[s,ai]*sum(FxnR[y,s,r,ai,])/Cmax[s,ai] } } #OLD VERSION (1/29/13) #{ w.pred[y,s,r,ai] = w.pred[y,s-1,r,ai] + sum(FxnR[y,s,r,ai,]*mean(ce)) } } }}}} #Store arrays for use outside of a function w.pred <<- w.pred } ## CALCULATE OBSERVED AND EXPECTED DIETS - for graphing and calculating RSS #"Expected" values of diets, based on empirical estimates from ChesMMAP if(preys==2) { diet.1.exp = 1 / (1 + exp(-(-0.5)*(seq(0.25,ages,by=0.25) - 1.5))) # see Flounder_Diets by age_1-29-13.xlsx for equation diet.2.exp = 1 / (1 + exp(-(0.5)*(seq(0.25,ages,by=0.25) - 1.5))) # see Flounder_Diets by age_1-29-13.xlsx for equation } if(preys==3) { diet.1.exp = 1 / (1 + exp(-(-1)*(seq(0.25,ages,by=0.25) - 1))) # see Flounder_Diets by age_1-29-13.xlsx for equation diet.3.exp = 0.9 / (1 + exp(-(1)*(seq(0.25,ages,by=0.25) - 3.5))) # see Flounder_Diets by age_1-29-13.xlsx for equation diet.2.exp = 1 - (diet.1.exp + diet.3.exp) # see Flounder_Diets by age_1-29-13.xlsx for equation #plot(seq(0.25,ages,by=0.25),diet.1.exp, type="l", ylim=c(0,1)) #lines(seq(0.25,ages,by=0.25),diet.2.exp, col="red") #lines(seq(0.25,ages,by=0.25),diet.3.exp, col="blue") diet.1.exp <<- diet.1.exp; diet.2.exp <<- diet.2.exp; diet.3.exp <<- diet.3.exp #Make accessible outside of function for additional plotting. } #"Observed" values of diets, based on the simulation diet.1.obs = apply(FxnR[,,1,,1]/apply(FxnR[,,1,,],c(1,2,3),sum),c(2,3),mean) # calculate mean diet proportion across years diet.2.obs = apply(FxnR[,,1,,2]/apply(FxnR[,,1,,],c(1,2,3),sum),c(2,3),mean) if(preys==3) diet.3.obs = apply(FxnR[,,1,,3]/apply(FxnR[,,1,,],c(1,2,3),sum),c(2,3),mean) ## GRAPHING GROWTH OUTPUT if(fxn==T && graph==T) { #This section is essentially commented out if i am not running growthfxn(). This code is copied below for when i run PredPreySim() windows() #open new graph window for plotting width=7,height=12 #Temporary plotting to investigate Growth trajectories and Cmax par(mfrow=c(2,2)) x.ages = seq(0.25,ages,by=0.25) #Weight at age plot plot(x.ages,w, type="l",lwd=2, col="red", main="Weight at age",ylab="Weight (kg)",xlab="Age",xlim=c(0,12),ylim=c(0,4)) for (y in (ages+1):years){ lines(x.ages,w.pred[y,,1,], col="blue") } points(LWdata$AGE2,LWdata$WEIGHT_KG) legend("topleft",c("Mean from wild","Simulated"), col=c("red","blue"), lty=1, lwd=c(2,1)) #Consumption relative to Cmax by age plot(x.ages,Cmax, type="l",lwd=3, col="red", ylim=c(0,12), main="Consumption at age", ylab="Consumption (kg/fish/season)",xlab="Age" ) for (y in (ages+1):years){ lines(x.ages,apply(FxnR[y,,1,,],c(1,2),sum)) } #calculates sum of prey consumed by season and age (for each year y, and region 1) legend("topleft",c("Cmax","Consumption"), col=c("red","black"), lty=1, lwd=c(3,1)) #Relative diet by age plot(x.ages,rep(-50,ages*seasons), ylim=c(0,100),ylab="Diet (%W)",xlab="Age", main="Estimated diet") #Create empty plot for (y in (ages+1):years){ lines(x.ages,100*FxnR[y,,1,,1]/apply(FxnR[y,,1,,],c(1,2),sum),col="black") lines(x.ages,100*FxnR[y,,1,,2]/apply(FxnR[y,,1,,],c(1,2),sum),col="red") if(preys==3) lines(x.ages,100*FxnR[y,,1,,3]/apply(FxnR[y,,1,,],c(1,2),sum),col="blue") } lines(x.ages, diet.1.exp*100, col="black", lwd=4) lines(x.ages, diet.2.exp*100, col="red", lwd=4) if(preys==3) lines(x.ages, diet.3.exp*100, col="blue", lwd=4) legend("right",c("Prey 1","Prey 2","Prey 3"), col=c("black","red","blue"), lty=1, lwd=1, bg="white") #Time series of prey abundance plot(seq(1,by=0.25,length=years*seasons),t(N.prey[,,1,1]), type="l", lwd=2, main=paste("Prey Abundances, Reg 1"),xlab="Year",ylab="Prey Abundance (numbers)") par(new=T) plot(seq(1,by=0.25,length=years*seasons),t(N.prey[,,1,2]), type="l", lwd=2, col="red", axes=F, xlab="",ylab="") axis(4, col="red") legend("topleft",c("Prey 1","Prey 2"), col=c("black","red"), lty=1, bg="white") } #End Plotting section ##CALCULATE RSS for weight-at-age and %diet composition - needed to optimize attackrates() below if(fxn==T) { RSS = sum((apply(w.pred[,,1,], c(2,3), mean) - w)^2) MSE.w = RSS/length(w) #Mean square error #if(optim.w==T) print(RSS) #These were commented out temporarily as I tried to troubleshoot why teh optim function wasn't working below. (2-5-13) #RSS.diets = sum((diet.1.obs - diet.1.exp)^2) + sum((diet.2.obs - diet.2.exp)^2) if(preys==3) { RSS.diets = sum((diet.1.obs - diet.1.exp)^2) + sum((diet.2.obs - diet.2.exp)^2) + sum((diet.3.obs - diet.3.exp)^2) MSE.diets = RSS.diets/(length(w)*preys) #the three is to account for the 3 prey groups } #if(optim.diet==T) print(RSS.diets) if(optim.both==T) print(MSE.w + MSE.diets) #print(RSS + RSS.diets) #if(optim.weighted ==T ) print( wt4optim * RSS + RSS.diets) } #end RSS section #} #END growthfxn() #Run the growth function to generate values if(FALSE) { #COMMENT OUT #OLD DEFAULTS: atck.max.1=exp(-19.5),atck.max.2=exp(-15.5), atck.max.3=exp(-18), atck.rho.1 = -0.2, atck.rho.2 = 3, atck.rho.3 = 0, # atck.infl.1= 2, atck.infl.2 = 2, atck.infl.3 = 0, growthfxn() } #COMMENT OUT ################################################# # CODE FOR OPTIMIZING ATTACK RATES ################################################# if(FALSE) { #COMMENT OUT ALL OPTIMIZATION CODE WHEN RUNNING PredPreySim() #NEW 8/27/13 - Function for optimizing attack rates (with prey 2 attack rates modeled with a double logistic) attackrates = function(x) { #x is a vector of parameter values (NOTE the unique and critical order of the inputs!) x.1 = exp(x[1]) #max attackrate for prey 1 NOTE: easier to optimize the log of the parameters because they are so small! x.2 = exp(x[2]) #max attackrate for prey 2 x.3 = exp(x[3]) #max attackrate for prey 3 x.4 = x[4] #rate of reaching max attack rate (atck.rho) for prey 1 NOTE: no need to exp x.5 = x[5] #rate of ascending and descending attack rates (atck.rho) for prey 2 (assume function is symetric) x.6 = x[6] #rate of reaching max attack rate (atck.rho) for prey 3 #Run the growth function to get the MSE for the weight at age and %diet fits, given the attack rate parameters defined by x #Note: to optimize BOTH w@age and %diet, set optim.w AND optim.diet to FALSE. growthfxn can only spit out one RSS value for optim to work. #Set the other two to false depending on what you want to optimize (optim.w , optim.diet , optim.both) growthfxn(atck.max.1=x.1, atck.max.2a=x.2, atck.max.3=x.3, atck.rho.1=x.4, atck.rho.2a=x.5, atck.rho.2b=x.5, atck.rho.3=x.6, #Note: setting atck.rho.2a and 2b equal. atck.infl.3=3.5, optim.w=F, optim.diet=F, optim.both=T, fxn=T, graph=F) } #end attackrates() #OLD Function for optimizing attack rates (2/5/13) attackrates = function(x) { #x is a vector of parameter values (NOTE the unique and critical order of the inputs!) x.1 = exp(x[1]) #max attackrate for prey 1 NOTE: easier to optimize the log of the parameters because they are so small! x.2 = exp(x[2]) #max attackrate for prey 2 x.3 = exp(x[3]) #max attackrate for prey 3 x.4 = x[4] #rate of reaching max attack rate (atck.rho) for prey 1 NOTE: no need to exp x.5 = x[5] #rate of reaching max attack rate (atck.rho) for prey 2 x.6 = x[6] #rate of reaching max attack rate (atck.rho) for prey 3 #Run the growth function to get the MSE for the weight at age and %diet fits, given the attack rate parameters defined by x #Note: to optimize BOTH w@age and %diet, set optim.w AND optim.diet to FALSE. growthfxn can only spit out one RSS value for optim to work. #Set the other two to false depending on what you want to optimize (optim.w , optim.diet , optim.both) growthfxn(atck.max.1=x.1, atck.max.2=x.2, atck.max.3=x.3, atck.rho.1=x.4, atck.rho.2=x.5, atck.rho.3=x.6, atck.infl.3=3, optim.w=F, optim.diet=F, optim.both=T, fxn=T, graph=F, Gmax.a=rep(0.6,ages)) } #end attackrates() #Exploration and Troublshooting growthfxn() growthfxn(atck.max.1=exp(-19),atck.max.2=exp(-16), atck.rho.1=-0.2, atck.rho.2=1, fxn=T, graph=T) attackrates(c(-19,-16, -0.2, 1 , -18, 1)) ##OPTIMIZE FXN RESPONSE PARMS - FOR 3 PREY #Optimize attack rates and rho for w@age AND diets (preys=3) WITH Double-logistic prey 2 attack rate (8-27-13) test = optim(c(-19,-16, -14, -1.2, 1 , 1.7), attackrates,lower=c(-Inf,-Inf,-Inf,-3,-Inf,0), upper=c(0,0,0,0,Inf,3), method="L-BFGS-B",hessian=T) test #Best-fit parameters; 8/27/13 (Using double-logistic prey 2 attack rates; atck.infl=3.5) growthfxn(atck.max.1=exp(-19.5), atck.max.2a=exp(-16.3), atck.max.3=exp(-13.5), atck.rho.1=-1.1, atck.rho.2a=0.46, atck.rho.2b=0.46,atck.rho.3=1.9, atck.infl.3=3.5, graph=T ) #OLD results (2-5-13) #Optimize attack rates and rho for w@age AND diets (preys=3) (2-5-13) # test = optim(c(-19,-16, -14, -1.2, -0.4 , 1.7), attackrates,lower=c(-Inf,-Inf,-Inf,-3,-3,0), upper=c(0,0,0,0,0,3), method="L-BFGS-B",hessian=T) # test #growthfxn(atck.max.1=exp(-19.5), atck.max.2=exp(-16.7), atck.max.3=exp(-13.7), # atck.rho.1=-0.7, atck.rho.2=0, atck.rho.3=2.4, # atck.infl.3=3, # Gmax.a = rep(0.6,ages), graph=T ) } #END COMMENT OUT FOR OPTIMIZATION SECTION ################################### # POPULATION SIMULATION ################################### ## CREATE EMPTY MATRICES TO POPULATE WITH DATA. #Array with seasonal and regional Spawning stock biomasses (S.r) S.r = array(NA, dim=c(years,seasons,regions),dimnames=list(1:years,paste("s",1:seasons,sep=""),paste("reg",1:regions,sep=""))) #Array with seasonal Spawning stock biomass (summed across regions) (S) S = array(NA, dim=c(years,seasons),dimnames=list(1:years,paste("s",1:seasons,sep=""))) #Array with predator abundance (P) (UNITS: no. of fish) P = array(NA, dim=c(years,seasons,regions,ages),dimnames=list(1:years,paste("s",1:seasons,sep=""),paste("reg",1:regions,sep=""),paste("age",0:(ages-1),sep=""))) #Array with fishery yields (UNITS: kg) NEW (2/1/13) Y = array(NA, dim=c(years,seasons,regions,ages),dimnames=list(1:years,paste("s",1:seasons,sep=""),paste("reg",1:regions,sep=""),paste("age",0:(ages-1),sep=""))) #NEW (1/29/13) - Array with predator weights-at-age over time and space. (UNITS = kg) w.pred <<- array(NA,dim=c(years,seasons,regions,ages),dimnames=list(1:years,paste("s",1:seasons,sep=""),paste("reg",1:regions,sep=""),paste("age",0:(ages-1),sep=""))) for (r in 1:regions) (w.pred[1,,r,] = w ) #Set first year values to previously defined start values. ## INITIATE POPULATION AND SSB IN YEAR 1, based on previously defined P0 #Initiate Population and SSB in year 1 based on previously defined P0. P[1,,,] = P0[1,,,] for (s in 1:seasons) { for (r in 1:regions) { S.r[1,s,r] = sum(P[1,s,r,] * w.pred[1,s,r,] * m[s,]) } #SSB for each year, season and region (summed across ages) S[1,s] = sum(S.r[1,s,]) } #SSB for each year and season (summed across regions and ages) #S.r[1,,] #S[1,] #Arrays for use in the loop #P [years, seasons, regions, ages] #S.r[years, seasons, regions] #S [years, seasons] #phi.in[ages,r,n,seasons] #JUNK y=2; s=1; r=1; nr=2; ai=2; M=0.02; F.a[,1] = c(0.5,1.0,4) ### RUN POPULATION DYNAMICS LOOP ### # Order of operations [ RECRUITMENT - MORTALITY - GROWTH - MOVEMENT - CENSUS ] for (y in 2:years){ for (s in 1:seasons) { for (ai in 1:ages) { for (r in 1:regions) { #Calculate recruits AND Weight-at-age for age-0 fish if (s1 && ais.recruit) || (ai>1 && s>1)) { #Region1 #Number of fish P[y,s,1,ai] = ( P[y,s-1,1,ai] * exp(-0.25*(M+F.a[y,ai])) #account for mortality * (1-sum(phi.in[ai,,1,s])) #proportion of survivors staying in the region (ie account for losses to other regions) + P[y,s-1,2,ai] * exp(-0.25*(M+F.a[y,ai])) * phi.in[ai,1,2,s] #survivors of region 2 that move into region 1 + 0 ) #ADD ADDITIONAL LINES IF ADDING REGIONS!! #NEW 8/27/13 - Weight-at-age (weighted mean of fish weights: w.pred[y,s,1,ai] = (N1*w1 + N2*w2) / (N1+N2) w.pred[y,s,1,ai] = ((( P[y,s-1,1,ai] * exp(-0.25*(M+F.a[y,ai])) * (1-sum(phi.in[ai,,1,s])) #Number of surviving fish staying in region 1 * (w.pred[y,s-1,1,ai] + G.func(WT=w.pred[y,s-1,1,ai],C=sum(FxnR[y,s,1,ai,]))*f.Temp[1,s]) ) #Weight-at-age of fish after growing in region 1 + ( P[y,s-1,2,ai] * exp(-0.25*(M+F.a[y,ai])) * phi.in[ai,1,2,s] #Number of surviving fish coming from region 2 * (w.pred[y,s-1,2,ai] + G.func(WT=w.pred[y,s-1,2,ai],C=sum(FxnR[y,s,2,ai,]))*f.Temp[2,s]) )) #Weight-at-age of fish after growing in region 2 / P[y,s,1,ai] ) #Total number of fish in region 1 if(FALSE) { #Comment out OLD Gmax Method #OLD 1/29/13 - Weight-at-age (weighted mean of fish weights: w.pred[y,s,1,ai] = (N1*w1 + N2*w2) / (N1+N2) w.pred[y,s,1,ai] = ((( P[y,s-1,1,ai] * exp(-0.25*(M+F.a[y,ai])) * (1-sum(phi.in[ai,,1,s])) #Number of surviving fish staying in region 1 * (w.pred[y,s-1,1,ai] + Gmax[s,ai]*sum(FxnR[y,s,1,ai,])/Cmax[s,ai]) ) #Weight-at-age of fish after growing in region 1 + ( P[y,s-1,2,ai] * exp(-0.25*(M+F.a[y,ai])) * phi.in[ai,1,2,s] #Number of surviving fish coming from region 2 * (w.pred[y,s-1,2,ai] + Gmax[s,ai]*sum(FxnR[y,s,2,ai,])/Cmax[s,ai]) )) #Weight-at-age of fish after growing in region 2 / P[y,s,1,ai] ) #Total number of fish in region 1 } #Region2 #Number of fish P[y,s,2,ai] = ( P[y,s-1,2,ai] * exp(-0.25*(M+F.a[y,ai])) #account for mortality * (1-sum(phi.in[ai,,2,s])) #proportion of survivors staying in the region (ie account for losses to other regions) + P[y,s-1,1,ai] * exp(-0.25*(M+F.a[y,ai])) * phi.in[ai,2,1,s] #survivors of region 1 that move into region 2 + 0 ) #ADD ADDITIONAL LINES IF ADDING REGIONS!! #NEW 8/27/13 - Weight-at-age (weighted mean of fish weights: w.pred[y,s,1,ai] = (N1*w1 + N2*w2) / (N1+N2) w.pred[y,s,2,ai] = ((( P[y,s-1,2,ai] * exp(-0.25*(M+F.a[y,ai])) * (1-sum(phi.in[ai,,2,s])) #Number of surviving fish staying in region 2 * (w.pred[y,s-1,2,ai] + G.func(WT=w.pred[y,s-1,2,ai],C=sum(FxnR[y,s,2,ai,]))*f.Temp[2,s]) ) #Weight-at-age of fish after growing in region 2 + ( P[y,s-1,1,ai] * exp(-0.25*(M+F.a[y,ai])) * phi.in[ai,2,1,s] #Number of surviving fish coming from region 1 * (w.pred[y,s-1,1,ai] + G.func(WT=w.pred[y,s-1,1,ai],C=sum(FxnR[y,s,1,ai,]))*f.Temp[1,s]) )) #Weight-at-age of fish after growing in region 1 / P[y,s,2,ai] ) #Total number of fish in region 2 if(FALSE) { #Comment out OLD Gmax Method #OLD 1/29/13 - Weight-at-age (weighted mean of fish weights: w.pred[y,s,1,ai] = (N1*w1 + N2*w2) / (N1+N2) w.pred[y,s,2,ai] = ((( P[y,s-1,2,ai] * exp(-0.25*(M+F.a[y,ai])) * (1-sum(phi.in[ai,,2,s])) #Number of surviving fish staying in region 2 * (w.pred[y,s-1,2,ai] + Gmax[s,ai]*sum(FxnR[y,s,2,ai,])/Cmax[s,ai]) ) #Weight-at-age of fish after growing in region 2 + ( P[y,s-1,1,ai] * exp(-0.25*(M+F.a[y,ai])) * phi.in[ai,2,1,s] #Number of surviving fish coming from region 1 * (w.pred[y,s-1,1,ai] + Gmax[s,ai]*sum(FxnR[y,s,1,ai,])/Cmax[s,ai]) )) #Weight-at-age of fish after growing in region 1 / P[y,s,2,ai] ) #Total number of fish in region 2 } } #Calculate seasonal SSB (to calculate recruits for following year) S.r[y,s,r] = sum(P[y,s,r,] * w.pred[y,s,r,] * m[s,]) #SSB for each year, season and region S[y,s] = sum(S.r[y,s,]) #SSB for each year, season, (summed across regions) #Calculate fishery yield (kg) - using Baranov's Catch equation Y[y,s,r,ai] = F.a[y,ai]/(F.a[y,ai] + M) * (1 - exp(-0.25*(M+F.a[y,ai]))) * P[y,s,r,ai] * w.pred[y,s,r,ai] }}}} #END POP DY LOOP #P[,3,1,]; P[,,,1]; S; N.prey; w.pred[1:10,,1,] #Checking matrices ## Calculate annual fishery landings Y.yr = apply(Y,c(1),sum) ## Calculate annual recruitment (R) (numbers of fish from both regions in recruitment season) (9-30-13) R = apply(P[,s.recruit,,1],1,sum) #Store matrices for access outside of the function S.r <<- S.r S <<- S P <<- P w.pred <<- w.pred Y <<- Y Y.yr<<- Y.yr R <<- R ## GRAPHING GROWTH OUTPUT if(graph==T) { windows() #open new graph window for plotting width=7,height=12 #Temporary plotting to investigate Growth trajectories and Cmax par(mfrow=c(2,2)) x.ages = seq(0.25,ages,by=0.25) #Weight at age plot plot(x.ages,w, type="l",lwd=2, col="red", main="Weight at age",ylab="Weight (kg)",xlab="Age",xlim=c(0,12),ylim=c(0,4)) for (y in (ages+1):years){ lines(x.ages,w.pred[y,,1,], col="blue") } points(LWdata$AGE2,LWdata$WEIGHT_KG) legend("topleft",c("Mean from wild","Simulated"), col=c("red","black"), lty=1, lwd=c(2,1)) #Consumption relative to Cmax by age plot(x.ages,Cmax, type="l",lwd=3, col="red", main="Consumption at age", ylab="Consumption (kg/fish/season)",xlab="Age" ) for (y in (ages+1):years){ lines(x.ages,apply(FxnR[y,,1,,],c(1,2),sum)) } #calculates sum of prey consumed by season and age (for each year y, and region 1) legend("topleft",c("Cmax","Consumption"), col=c("red","black"), lty=1, lwd=c(3,1)) #Relative diet by age plot(x.ages,rep(-50,ages*seasons), ylim=c(0,100),ylab="Diet (%W)",xlab="Age", main="Estimated diet") #Create empty plot for (y in (ages+1):years){ lines(x.ages,100*FxnR[y,,1,,1]/apply(FxnR[y,,1,,],c(1,2),sum),col="black") lines(x.ages,100*FxnR[y,,1,,2]/apply(FxnR[y,,1,,],c(1,2),sum),col="red") if(preys==3) lines(x.ages,100*FxnR[y,,1,,3]/apply(FxnR[y,,1,,],c(1,2),sum),col="blue") } lines(x.ages, diet.1.exp*100, col="gray40", lwd=4) lines(x.ages, diet.2.exp*100, col="red4", lwd=4) if(preys==3) lines(x.ages, diet.3.exp*100, col="midnightblue", lwd=4) legend("right",c("Prey 1","Prey 2","Prey 3"), col=c("black","red","blue"), lty=1, lwd=1, bg="white") #Time series of prey biomass plot(seq(1,by=0.25,length=years*seasons),t(N.prey[,,1,1]), type="l", lwd=2, main=paste("Prey Biomass, Reg 1"),xlab="Year",ylab="Prey Abundance (fish)") par(new=T) plot(seq(1,by=0.25,length=years*seasons),t(N.prey[,,1,2]), type="l", lwd=2, col="red", axes=F, xlab="",ylab="") axis(4, col="red") par(new=T) plot(seq(1,by=0.25,length=years*seasons),t(N.prey[,,1,3]), type="l", lwd=2, col="blue", axes=F, xlab="",ylab="") legend("topleft",c("Prey 1","Prey 2", "Prey 3"), col=c("black","red","blue"), lty=1, bg="white") } #End Plotting section ### GRAPHING OF MODEL OUTPUT #P[years,seasons,regions,ages] #S[years,seasons] #P.gr.1[years,seasons,regs] #P.gr.2[years,seasons,ages] if(graph==T) { yr.min = ages+1 if(FALSE){ #Commenting out less important graphs #Plot of seasonal total abundances by years and regions #Sum abundances across all ages P.gr.1 = apply(P,c(1,2,3),sum) windows() par(mfrow=c(2,2)) for (s in 1:4) { plot(P.gr.1[,s,1],type="l",lwd=3, col="black",ylim=c(min(P.gr.1[,s,]),max(P.gr.1[,s,])), ylab="Total Abundance (all ages)",xlab="Year") lines(P.gr.1[,s,2], lwd=3, col="red") legend("bottomright",c("Reg 1","Reg 2"), col=c("black","red"),lty=1,lwd=3,merge=T,bg="white") title(paste("Season",s)) } #Plot of age-specific abundances by years (summed across regions) #Sum abundances across all regions P.gr.2 = apply(P,c(1,2,4),sum) windows() par(mfrow=c(2,2)) for (s in 1:4) { plot(P.gr.2[,s,1],type="l",lwd=3, col="black", ylim=c(min(P.gr.2[,s,]),max(P.gr.2[,s,])),ylab="Abundance (all regions)",xlab="Year") lines(P.gr.2[,s,2], lwd=3, col="red") lines(P.gr.2[,s,3], lwd=3, col="green") legend("bottomright",c("Age 1","Age 2","Age 3"), col=c("black","red","green"),lty=1,lwd=3,merge=T,bg="white") title(paste("Season",s)) } #Plot of population proportion in region 1 for each season #Calculate proportion for year 50 only P.gr.3 = round(P[50,,1,] / P.gr.2[50,,] * 100 ,1) #Note: using P.gr.2 from above windows() par(mfrow=c(1,1)) plot(0:(ages-1),P.gr.3[1,] ,type="l",lwd=3, col="blue", ylim=c(0,100),main="Population % in Reg 1 by age", ylab="% of population in Reg 1",xlab="Age", cex.lab=1.5) lines(0:(ages-1),P.gr.3[2,], lwd=3, col="green") lines(0:(ages-1),P.gr.3[3,], lwd=3, col="red") lines(0:(ages-1),P.gr.3[4,], lwd=3, col="black") legend("right",c("Seas 1","Seas 2","Seas 3","Seas 4"), col=c("blue","green","red","black"),lty=1,lwd=3,bg="white") } #End of commented out graphs windows() #Set new graph window for plotting par(mfrow=c(2,2)) #Plot of SSB by years during season 3 (summed across regions and ages) #plot(S[(ages+1):years,3],type="l",lwd=3, col="black",ylab="SSB (kg)",xlab="Year", main="SSB") #Removes initial burn-in period of 1 full generation (ages+1) plot(S[,s.spawn]/1E6,type="l",lwd=3, col="black",ylab="SSB (kmt)",xlab="Year", main="SSB") abline(h=91.4, lty=2, col="red") #This is the S.target based on Fmin=Fmsy=0.31 (9/19/13) #Plots of Recruitment #Calculate recruits to each region based on the SSB in matrix S R.gr.1 = array(data=NA,dim=c(years,regions),dimnames=list(1:years,paste("reg",1:regions,sep=""))) for (y in 2:years) { for (r in 1:regions) { R.gr.1[y,r] = (((alpha.SR * S[y-1,s.spawn]/1E6) / (1 + (S[y-1,s.spawn]/1E6)/K.SR)^beta.SR *1E6 * #Calculate age-0 recruits from BH equation if it is the recruitment season. Note: 1E3 and 1E6 are needed because the SR function deals with SSB in mt and R in millions of age-0 fish. theta[r,] * R.deviates[y])) #* exp(-0.25*(M+F.a[y,ai])) ) #allocate recruits to each region, and add stochastic variability (recruitment deviates) #account for mortality }} #Plot Recruits by year and region plot(R.gr.1[,1]/1E6,type="l",lwd=3, col="black", ylim=c(0,max(na.omit(R.gr.1))/1E6),ylab="Number of recruits (millions)",xlab="Year") lines(R.gr.1[,2]/1E6, lwd=3, col="red") legend("bottomright",c("Reg 1","Reg 2"), col=c("black","red"),lty=1,lwd=3,merge=T,bg="white") title("Recruits by year and region") #Plot S-R relationship for observed SSB S.junk = seq(0,max(S[,s.spawn]/1E6),by=1) #Range of SSB values for plotting the predicted recruitment based on the SR function. (units=kmt) #Recruits come from Season 1 of P array. [Note that the "-years" and "-1" are to match up the years of SSB with Year of recruitment (since there is a delay in when Recruits show up!] plot(S[-years, s.spawn]/1E6,rowSums(P[-1,1,,1])/1E6,type="p",lwd=2, col="black", ylim=c(0,max(rowSums(P[,1,,1])/1E6)),xlim=c(0,max(S[,s.spawn])/1E6), ylab="Number of recruits (millions)",xlab="SSB (kmt)") abline(h=alpha.SR*K.SR * exp(-0.25*(M+mean(F.a[,1]))),col="gray") #Maximum production of recruits (after accounting for mortality in the first timestep) lines(S.junk,alpha.SR * (S.junk) / (1 + S.junk/K.SR)^beta.SR * exp(-0.25*(M+mean(F.a[,1]))) , col="red",lwd=2 ) title("S-R relationship") #Plot of annual fishery landings plot(Y.yr/1E6, xlab="Year", ylab="Fishery yield (kmt)", type="b", lwd=2, main="Fishery Landings") if(substr(P.scenario,1,5)=="Pulse") { abline(v=pulse.years, lty=2, col="red") legend("topleft",c("Pulse year"), col="red",lty=2, bg="white") } } #End of Plotting Section } ## END OF PREDPREYSIM FUNCTION PredPreySim() PredPreySim(N.prey.dev=0, R.dev.sd=0, Fmin=0.31, F.scenario="decrF") PredPreySim(N.prey.dev=0, R.dev.sd=0, F.scenario="mor.F") PredPreySim(P.scenario="PulsB2") PredPreySim(P.scenario="Pulse1") PredPreySim(P.scenario="Pulse2", F.scenario="highF") PredPreySim(P.scenario="Step.1", seed=10); ProdStepIncr PredPreySim(P.scenario="Step.2", seed=10); ProdStepIncr PredPreySim(P.scenario="Step.3", seed=10); ProdStepIncr PredPreySim(Fmax=0,Fmin=0, R.dev.sd=0, N.prey.dev=0) S #9/14/13 - Comparing model SSB trends under decrF scenarios to the assessment trends in SSB. plot(S[,s.spawn]/1E6,type="l",lwd=3, col="black",ylab="SSB (kmt)",xlab="Year", main="SSB - model") abline(h=91,col="red", lty=2) plot(SRdata$Year,SRdata$SSB_kmt,type="l",lwd=3, col="red",ylab="SSB (kmt)",xlab="Year", main="SSB - assessment") plot(S[,s.spawn]/1E6,type="l",lwd=3, col="black",ylab="SSB (kmt)",xlab="Year", main="SSB", xlim=c(10,50)) #Remove initial burn-in period of 1 full generation abline(h=106,col="red", lty=2) par(new=T) plot(SRdata$ID+10,SRdata$SSB_kmt,type="l",lwd=3, col="red",ylab="SSB (kmt)",xlab="Year", main="SSB", xlim=c(10,40), ylim=c(10,65)) #Remove initial burn-in period of 1 full generation #This is the code used to establish S.target: (for troubleshooting) years=80 PredPreySim(R.dev.sd=0,N.prey.dev=0, F.scenario="decrF") #Run simulation without stochasticity S.target = max(S[,3]/1E6)*1; S.target #Target SSB defined as 100% of the 80 year equilibrium value. S.target=112.2 kmt (95% of equilibrium is 106.6) #Obtaining MSY from model at different F's to determing the best Fmsy and corresponding Smsy Fseq = c(0,0.1,0.2,0.255,seq(0.26,0.35,by=0.01),0.5,1,1.5) #Generate vector of F's to use Y.F = data.frame(Fmax = Fseq, Yield_kmt = NA) #Create dataframe to store yield values for a given F. for (i in Fseq) { #Run loop to calculate Yield for each F PredPreySim(Fmax=i, N.prey.dev=0, R.dev.sd=0, F.scenario="highF", graph=F) Y.F[Y.F$Fmax==i,"Yield_kmt"] = Y.yr[60]/1E6 print(paste("i =",i)) } plot(Y.F$Fmax, Y.F$Yield_kmt, type="b", xlab="F", ylab="Yield (kmt)", main="MSY by F rate") #Plot Y vs. F curve to find MSY and FMSY Y.F #FMSY = 0.31 (with Smsy = 91.4 kmt) #F used for target (as of 9-13-13): F=0.255 (with Smsy = 106.2 kmt (at season 4); 108.5 kmt at season 3) #The difference in Yield at this Fmsy vs F=0.255 was only 1.00%; but with F=0.255 your S.target is 16.2% greater PredPreySim(Fmax=0.31, N.prey.dev=0, R.dev.sd=0, F.scenario="highF") S[60,s.spawn]/1E6 #Smsy = 91.4 kmt PredPreySim(Fmax=0.255, N.prey.dev=0, R.dev.sd=0, F.scenario="highF") S[60,s.spawn]/1E6 #Smsy = 106.2 kmt #Exploring different F.scenarios (linear decr; exponential, moratorium, assessment) (9-19-13) PredPreySim(R.dev.sd=0.4,N.prey.dev=0, F.scenario="decrF") PredPreySim(R.dev.sd=0.4,N.prey.dev=0, F.scenario="exp.F") PredPreySim(R.dev.sd=0.4,N.prey.dev=0, F.scenario="ass.F") PredPreySim(R.dev.sd=0.4,N.prey.dev=0, F.scenario="mor.F") PredPreySim(R.dev.sd=0.4,N.prey.dev=0, F.scenario="knifF") #Exploring the prey biomass trajectories to see if the % increases in prey biomass are reasonable. (9-23-13) #CONCLUSIONS: # with pulse.range = c(2:10), the prey biomass could increase as much as 30X, but averaged a max increase of 16. # Decided to re-run the simulation with pulse.range = c(2:6) where biomass would increase up to 18X, but the average maximum value across sims was 10.4 # See "PulseRange effect on biomass pulses_Exploration_9-24-13.xlsx" par(mfrow=c(2,2)) for (i in 1:50){ PredPreySim(P.scenario="Pulse1", graph=F, seed=500*i, pulse.range=c(2:6), N.prey.dev=0.4) print(max(N.prey[,,1,1])/N.prey.1.mean) preyX=1 windows() #plot(seq(1,by=0.25,length=years*seasons),t(N.prey[,,1,preyX])/N.prey.1.mean, type="l", xlab="Year", ylab="Increase (X fold)") #Normal scale # abline(h=c(1,10),col=c("red"),lty=2) plot(seq(1,by=0.25,length=years*seasons),log(t(N.prey[,,1,preyX])), type="l", xlab="Year", ylab="Increase (X fold)") #Lognormal scale abline(h=c(log(N.prey.1.mean),log(N.prey.1.mean*10)),col=c("red"),lty=2) } #Checking that the new code for N.all is correct. (9/24/13) plot(seq(1,by=0.25,length=years*seasons),t(N.all[1,,,1,"Pulse"])/N.prey.1.mean, type="l", xlab="Year", ylab="Increase (X fold)") lines(seq(1,by=0.25,length=years*seasons),t(N.all[1,,,1,"NoPulse"])/N.prey.1.mean, col="red", lty=2) abline(h=c(1,10),col=c("red"),lty=2) plot(seq(1,by=0.25,length=years*seasons),t(N.all[1,,,2,"Pulse"])/N.prey.2.mean, type="l", xlab="Year", ylab="Increase (X fold)") lines(seq(1,by=0.25,length=years*seasons),t(N.all[1,,,2,"NoPulse"])/N.prey.2.mean, col="blue", lty=1) abline(h=c(1,10),col=c("red"),lty=2) plot(seq(1,by=0.25,length=years*seasons),log(t(N.all[1,,,3,"Pulse"])/N.prey.3.mean), type="l", xlab="Year", ylab="Increase (X fold)") lines(seq(1,by=0.25,length=years*seasons),log(t(N.all[1,,,3,"NoPulse"])/N.prey.3.mean), col="red", lty=1) abline(h=c(1,10),col=c("red"),lty=2) #CALCULATE PROPORTION OF MAX CONSUMPTION - 10-7-13 #FxnR[years,seas,regs,age,prey] #FxnR[,,1,,] Cons.all = apply(FxnR[,,1,,],c(1,2,3),sum) #sum of all prey consumed in region 1, by year, seas, age Cons.all = aperm(Cons.all, c(2,3,1)) #re-arrange array so years are the third dimension (to facilitate code below) Pmax = Cons.all #Create Pmax array with same dimensions as Cons.all Pmax[,,] = 0 #Set all values to zero to avoid errors for (i in 1:years) { #Calculate the proportion of max consumption (Pmax) consumed for each season, age, and year Pmax[,,i] = Cons.all[,,i]/Cmax } Pmax[,,1] #Check Year 1 Pmax.min = apply(Pmax,c(1,2),min) #Calculate min, max, and mean Pmax values across all years for this single simulation. Pmax.max = apply(Pmax,c(1,2),max) Pmax.mean= apply(Pmax,c(1,2),mean) min(na.omit(Pmax.mean)); max(na.omit(Pmax.mean)) #range of means across ages: 28-73% #JUNK - TROUBLESHOOTING 8-27-13 #FxnR[years,seas,regs,age,prey] FxnR[1,,1,,] Cons.all = apply(FxnR[1,,1,,],c(1,2),sum) #Cons.all Pmax = Cons.all/Cmax #Proportion of max consumption by seas, age Pmax write.table(Pmax, "clipboard", sep="\t") Pmax.seas = Pmax*f.Temp[1,] #Proportion of max consumption after factoring in temp effect. #Checking Growth by seas and age w.pred[10,,1,] G.func(WT=w.pred[10,,1,], C=Cmax) G.func(WT=w, C=Cmax) K = KL / (1 + exp(-KR*(w[4,8] - KW))); K G.func(WT=3030.7, C=16649) w.pred[30:40,,1,] ############################################################################################################ ############################################################################################################ ### SIMULATION LOOP FUNCTION ############################################################################################################ ############################################################################################################ ### DEFINE SCENARIOS #Define the number of global scenarios to be investigated (e.g. prey 1 pulses, prey 2 pulses) ,"PulsB1","PulsB2","PulsB3" ,"Step.1","Step.2","Step.3" "exp.F","ass.F","knifF" scen.names = expand.grid(P.scen=c("Pulse1","Pulse2","Pulse3"),I.scen=c("mig","mix"),F.scen=c("highF","lowF","decrF","dec2F","mor.F")) # scen.names$scenario = c(1:nrow(scen.names)) scen.names$scen.name = paste(scen.names$P.scen, scen.names$I.scen, scen.names$F.scen, sep="-") scen.names = subset(scen.names, select=c("scenario","scen.name","P.scen","I.scen","F.scen")) #Re-order columns to match the previous version scenarios = nrow(scen.names) #Sub Scenarios - focused on differences in the rate of fishing decrease if(FALSE) { #COMMENT OUT UNLESS YOU WANT TO RUN THE SUBSCENARIOS "Step.1","Step.2","Step.3" scen.names = expand.grid(P.scen=c("Pulse1","Pulse2","Pulse3"),I.scen=c("mig"),F.scen=c("decrF.14","decrF.12","decrF.10")) scen.names$scenario = c(1:nrow(scen.names)) scen.names$scen.name = paste(scen.names$P.scen, scen.names$I.scen, scen.names$F.scen, sep="-") scen.names = subset(scen.names, select=c("scenario","scen.name","P.scen","I.scen","F.scen")) #Re-order columns to match the previous version scenarios = nrow(scen.names) } #end Comment out. #Exploring "Step" Scenarios - compare effect of how biomass increases are distributed temporally (pulse years vs. step increase) if(FALSE) { #COMMENT OUT UNLESS YOU WANT TO RUN THE SUBSCENARIOS scen.names = expand.grid(P.scen=c("Pulse1","Pulse2","Pulse3","Step.1","Step.2","Step.3"),I.scen=c("mig"),F.scen=c("decrF")) scen.names$scenario = c(1:nrow(scen.names)) scen.names$scen.name = paste(scen.names$P.scen, scen.names$I.scen, scen.names$F.scen, sep="-") scen.names = subset(scen.names, select=c("scenario","scen.name","P.scen","I.scen","F.scen")) #Re-order columns to match the previous version scenarios = nrow(scen.names) } #end Comment out. #For Troubleshooting: scen.names = expand.grid(P.scen=c("Pulse1","Pulse2","Pulse3"),I.scen=c("mig"),F.scen=c("highF","lowF","decrF")) # "highF", ,"mix" scen.names$scenario = c(1:nrow(scen.names)) scen.names$scen.name = paste(scen.names$P.scen, scen.names$I.scen, scen.names$F.scen, sep="-") scen.names = subset(scen.names, select=c("scenario","scen.name","P.scen","I.scen","F.scen")) #Re-order columns to match the previous version scenarios = nrow(scen.names) #scenarios = 2 #scen.names = data.frame(scenario=c(1:scenarios), # scen.name=c("Pulse1-mig-highF","Pulse2-mix-decrF"), P.scen=c("Pulse1","Pulse2"), I.scen=c("mig","mix"), F.scen=c("highF","decrF")) ### SIMULATION FUNCTION (RUNNING PREDPREYSIM MANY TIMES) ### simfxn = function (seed.init=500, sims=3, #Number of simulation runs to conduct graph=T, #Switch for turning graphs on or off, notch.gr=F, #Indicate whether boxplots have a notch #PredPreySim parameters: (copying them here so i can store the info as "function.call" when saving the simulation runs. M=0.25, #Instantaneous natural mortality rate (per year) Fmax = 1.5, #Maximum instantaneous fishing mortality. Previously 2.0; exploring how 1.5 looks since that was more the average. (See StockRecruitFxn_8-29-13.xlsx) Fmin = 0.31, #Minimum instantaneous fishing mortality. decrF.time = 12, #Number of years to go from Fmax to Fmin (F=0.255 is the target F as of Jan 2013; Changed to F=0.31 on 9/19/13). F.sel = c(0,0.1,0.5,rep(1.0,ages-3)), #Fishing selectivity by age R.dev.sd=0.4, #Standard deviation for the lognormal random recruitment error alpha.SR = 3.4, #alpha parameter for Shepherd S-R function (Rothschild et al. 2012) K.SR = 27.4, #K parameter for Shepherd S-R function (Rothschild et al. 2012) beta.SR = 1, #beta shape parameter for Shepherd S-R function. beta.SR=1 is for Beverton Holt. CTO= 22 , #Temperature for optimum consumption; assumed value (see eq in Hanson et al. 1997) CTM= 35 , #Maximum temp parameter for temperature equation; assumed value (see eq in Hanson et al. 1997) CQ = 2.5 , #Rate parameter for temperature equation; assumed value (see eq in Hanson et al. 1997) CA = 0.3 , #Intercept Parameter for allometric relationship of Cmax (Hanson et al. 1997; GrowthModel_ConvEfficiency_8-25-13.xlsx) Cmax = CA*W^CB CB =-0.2 , #Decay Parameter for allometric relationship of Cmax (Hanson et al. 1997; GrowthModel_ConvEfficiency_8-25-13.xlsx) Cmax = CA*W^CB KL = 0.5 , #Lambda parameter for maximum conversion effiency (K) sigmoidal model(GrowthModel_ConvEfficiency_8-25-13.xlsx) KR = -0.0014 , #Rate of decline for sigmoidal K model (GrowthModel_ConvEfficiency_8-25-13.xlsx) KW = -600 , #Weight at which KL/2 is attained in Sigmoidal K model (GrowthModel_ConvEfficiency_8-25-13.xlsx) n.pulse=5, #Number of prey pulses pulse.range = c(2:6), #Range of multiplier values for when a pulse in production occurs. E.g. 2-10x production in pulse years N.prey.2.mean = 3e+10 , #Number of prey 2 N.prey.dev = 0.4, #Standard deviation for the lognormal random error of prey populations relB.1 = 10 , #Number of prey 1 relative to prey 2; because treating each prey as 1g (for all prey), this is equivalent to relative biomass relB.3 = .1, #Number of prey 3 relative to prey 2; because treating each prey as 1g (for all prey), this is equivalent to relative biomass atck.max.1 =exp(-19.5), atck.rho.1= -1.1, atck.infl.1= 2, #Parameters for sigmoidal Attack rate for prey 1 (max value, rate value rho (neg for decr fxn), and inflection point (age at 50%max)) atck.max.2a=exp(-16.3), #Parameters for double-logistic attack rate for prey 2 (max.2a = scaling parameter for max value) atck.rho.2a=0.46, atck.rho.2b = 0.46, #Parameters for double-logistic attack rate for prey 2 (rho = rate parameter for ascending (2a) and descending (2b) limbs) atck.infl.2a=1, atck.infl.2b=4, #Parameters for double-logistic attack rate for prey 2 (infl = age of inflection point for ascending (2a) and descending (2b) limbs) (see Buchheister disseration Chapter 3) atck.max.3 =exp(-13.5), atck.rho.3= 1.9, atck.infl.3= 3.5 #Parameters for sigmoidal Attack rate for prey 3 (max value, rate value rho (neg for decr fxn), and inflection point (age at 50%max)) #OLD Gmax.a = rep(0.6,ages), #Age-specific maximum growth (kg/yr) possible for species #OLD atck.max.1 =exp(-19.5), atck.max.2=exp(-16.7), atck.max.3=exp(-13.7), #Max attack rate on prey 1 and 2 #OLD atck.rho.1 = -0.7, atck.rho.2 = 0.0, atck.rho.3=2.4, #Rate for reaching max #OLD atck.infl.1= 2, atck.infl.2 = 2, atck.infl.3=3 #Age at 50% max (ie the "inflection point") ) { ## CREATE EMPTY ARRAYS TO FILL WITH SIM RESULTS #Spawning stock biomass results across simulations (summed across regions, spawning season only) (S.all) (UNITS: kg) S.all <<- array(NA, dim=c(scenarios,sims,years,2),dimnames=list(scen.names[,2],1:sims,1:years,c("NoPulse","Pulse"))) S.diff.mean <<- array(NA, dim=c(sims,scenarios),dimnames=list(paste("sim",1:sims,sep=''),scen.names[,"scen.name"])) #SSB %diff for years 26:end S.diff.p.mean <<- array(NA, dim=c(sims,scenarios),dimnames=list(paste("sim",1:sims,sep=''),scen.names[,"scen.name"])) #SSB %diff for pulse years only #Weight at age array across sims (region 1 only) (w.all) (UNITS: kg) w.all <<- array(NA,dim=c(scenarios,sims,years,seasons,ages,2),dimnames=list(scen.names[,2],1:sims,1:years,paste("s",1:seasons,sep=""),paste("age",0:(ages-1),sep=""),c("NoPulse","Pulse"))) w.diff.mean <<- array(NA,dim=c(sims,ages,scenarios),dimnames=list(paste("sim",1:sims,sep=''),paste("age",0:(ages-1),sep=""),scen.names[,2])) #weight@age %diff for years 26:end w.diff.p.mean <<- array(NA,dim=c(sims,ages,scenarios),dimnames=list(paste("sim",1:sims,sep=''),paste("age",0:(ages-1),sep=""),scen.names[,2])) #weight@age %diff for pulse years only #Total Annual fishery yields (Y.all) (UNITS: kg) Y.all <<- array(NA, dim=c(scenarios,sims,years,2),dimnames=list(scen.names[,2],1:sims,1:years,c("NoPulse","Pulse"))) Y.diff.mean <<- array(NA, dim=c(sims,scenarios),dimnames=list(paste("sim",1:sims,sep=''),scen.names[,"scen.name"])) #Yield %diff for years 26:end Y.diff.p.mean <<- array(NA, dim=c(sims,scenarios),dimnames=list(paste("sim",1:sims,sep=''),scen.names[,"scen.name"])) #Yield %diff for pulse years only #Annual Recruitment by year (R.all) (UNITS: millions of fish) (9-30-13) R.all <<- array(NA, dim=c(scenarios,sims,years,2),dimnames=list(scen.names[,2],1:sims,1:years,c("NoPulse","Pulse"))) R.diff.mean <<- array(NA, dim=c(sims,scenarios),dimnames=list(paste("sim",1:sims,sep=''),scen.names[,"scen.name"])) #Recruitment %diff for years 26:50 R.diff.p.mean <<- array(NA, dim=c(sims,scenarios),dimnames=list(paste("sim",1:sims,sep=''),scen.names[,"scen.name"])) #Recruitment %diff for pulse years only #Storing the pulse years for each simulation (i.e. the years in which prey pulses occured) (pulses.all) pulses.all <<- array(NA, dim=c(sims,n.pulse),dimnames=list(1:sims,1:n.pulse)) #Storing the Prey Biomass time series (N.all) N.prey[years,seas,reg,prey] N.all <<- array(NA, dim=c(sims,years,seasons,preys,2),dimnames=list(1:sims,1:years,paste("s",1:seasons,sep=""),paste("prey",1:preys,sep=""),c("NoPulse","Pulse"))) ## RUN SCENARIO and SIMULATION LOOP #Run pulse and non-pulse conditions for each scenario. The Pulse results will ultimately be in relation to the non-pulse results for (sc in 1:scenarios) { for (i in 1:sims) { #NoPulse (P.Scenario 1) #run PredPreySim for NoPulse (with appropriate I and F scenarios) and store pertinent output PredPreySim(P.scenario=1, I.scenario=scen.names[sc,"I.scen"], F.scenario=scen.names[sc,"F.scen"], graph=F, seed=(seed.init*i), #PredPreySim parameters: M = M, #Instantaneous natural mortality rate (per year) Fmax = Fmax, #Maximum instantaneous fishing mortality. Previously 2.0; exploring how 1.5 looks since that was more the average. (See StockRecruitFxn_8-29-13.xlsx) Fmin = Fmin, #Minimum instantaneous fishing mortality. decrF.time = decrF.time, #Number of years to go from Fmax to Fmin F.sel = F.sel, #Fishing selectivity by age R.dev.sd = R.dev.sd, #Standard deviation for the lognormal random recruitment error alpha.SR = alpha.SR, #alpha parameter for Shepherd S-R function (Rothschild et al. 2012) K.SR = K.SR, #K parameter for Shepherd S-R function (Rothschild et al. 2012) beta.SR = beta.SR, #beta shape parameter for Shepherd S-R function. beta.SR=1 is for Beverton Holt. CTO= CTO , #Temperature for optimum consumption; assumed value (see eq in Hanson et al. 1997) CTM= CTM , #Maximum temp parameter for temperature equation; assumed value (see eq in Hanson et al. 1997) CQ = CQ , #Rate parameter for temperature equation; assumed value (see eq in Hanson et al. 1997) CA = CA , #Intercept Parameter for allometric relationship of Cmax (Hanson et al. 1997; GrowthModel_ConvEfficiency_8-25-13.xlsx) Cmax = CA*W^CB CB = CB , #Decay Parameter for allometric relationship of Cmax (Hanson et al. 1997; GrowthModel_ConvEfficiency_8-25-13.xlsx) Cmax = CA*W^CB KL = KL , #Lambda parameter for maximum conversion effiency (K) sigmoidal model(GrowthModel_ConvEfficiency_8-25-13.xlsx) KR = KR , #Rate of decline for sigmoidal K model (GrowthModel_ConvEfficiency_8-25-13.xlsx) KW = KW , #Weight at which KL/2 is attained in Sigmoidal K model (GrowthModel_ConvEfficiency_8-25-13.xlsx) n.pulse=n.pulse, #Number of prey pulses pulse.range = pulse.range, #Range of multiplier values for when a pulse in production occurs. E.g. 2-10x production in pulse years N.prey.2.mean = N.prey.2.mean , #Number of prey 2 N.prey.dev = N.prey.dev, #Standard deviation for the lognormal random error of prey populations relB.1 = relB.1 , #Number of prey 1 relative to prey 2; because treating each prey as 1g (for all prey), this is equivalent to relative biomass relB.3 = relB.3, #Number of prey 3 relative to prey 2; because treating each prey as 1g (for all prey), this is equivalent to relative biomass atck.max.1 =atck.max.1, atck.rho.1= atck.rho.1, atck.infl.1= atck.infl.1, #Parameters for sigmoidal Attack rate for prey 1 (max value, rate value rho (neg for decr fxn), and inflection point (age at 50%max)) atck.max.2a=atck.max.2a, #Parameters for double-logistic attack rate for prey 2 (max.2a = scaling parameter for max value) atck.rho.2a=atck.rho.2a, atck.rho.2b = atck.rho.2b, #Parameters for double-logistic attack rate for prey 2 (rho = rate parameter for ascending (2a) and descending (2b) limbs) atck.infl.2a=atck.infl.2a, atck.infl.2b=atck.infl.2b, #Parameters for double-logistic attack rate for prey 2 (infl = age of inflection point for ascending (2a) and descending (2b) limbs) (see Buchheister disseration Chapter 3) atck.max.3 =atck.max.3, atck.rho.3= atck.rho.3, atck.infl.3= atck.infl.3 #Parameters for sigmoidal Attack rate for prey 3 (max value, rate value rho (neg for decr fxn), and inflection point (age at 50%max)) ) S.all[sc,i,,1] = S[,s.spawn] #SSB Y.all[sc,i,,1] = Y.yr #Total annual fishery yields R.all[sc,i,,1] = R #Total annual recruitment w.all[sc,i,,,,1] = w.pred[,,1,] #Weight at age for region 1. N.all[i,,,,1] = N.prey[,,1,] #Prey Biomass (N) in region 1 for the NoPulse scenario. #Pulse Scenarios (P.scenario of "Pulse1", etc. (pulses in prey 1, 2, or 3) PredPreySim(P.scenario=scen.names[sc,"P.scen"], I.scenario=scen.names[sc,"I.scen"], F.scenario=scen.names[sc,"F.scen"], graph=F, seed=(seed.init*i), #PredPreySim parameters: M = M, #Instantaneous natural mortality rate (per year) Fmax = Fmax, #Maximum instantaneous fishing mortality. Previously 2.0; exploring how 1.5 looks since that was more the average. (See StockRecruitFxn_8-29-13.xlsx) Fmin = Fmin, #Minimum instantaneous fishing mortality. decrF.time = decrF.time, #Number of years to go from Fmax to Fmin F.sel = F.sel, #Fishing selectivity by age R.dev.sd = R.dev.sd, #Standard deviation for the lognormal random recruitment error alpha.SR = alpha.SR, #alpha parameter for Shepherd S-R function (Rothschild et al. 2012) K.SR = K.SR, #K parameter for Shepherd S-R function (Rothschild et al. 2012) beta.SR = beta.SR, #beta shape parameter for Shepherd S-R function. beta.SR=1 is for Beverton Holt. CTO= CTO , #Temperature for optimum consumption; assumed value (see eq in Hanson et al. 1997) CTM= CTM , #Maximum temp parameter for temperature equation; assumed value (see eq in Hanson et al. 1997) CQ = CQ , #Rate parameter for temperature equation; assumed value (see eq in Hanson et al. 1997) CA = CA , #Intercept Parameter for allometric relationship of Cmax (Hanson et al. 1997; GrowthModel_ConvEfficiency_8-25-13.xlsx) Cmax = CA*W^CB CB = CB , #Decay Parameter for allometric relationship of Cmax (Hanson et al. 1997; GrowthModel_ConvEfficiency_8-25-13.xlsx) Cmax = CA*W^CB KL = KL , #Lambda parameter for maximum conversion effiency (K) sigmoidal model(GrowthModel_ConvEfficiency_8-25-13.xlsx) KR = KR , #Rate of decline for sigmoidal K model (GrowthModel_ConvEfficiency_8-25-13.xlsx) KW = KW , #Weight at which KL/2 is attained in Sigmoidal K model (GrowthModel_ConvEfficiency_8-25-13.xlsx) n.pulse=n.pulse, #Number of prey pulses pulse.range = pulse.range, #Range of multiplier values for when a pulse in production occurs. E.g. 2-10x production in pulse years N.prey.2.mean = N.prey.2.mean , #Number of prey 2 N.prey.dev = N.prey.dev, #Standard deviation for the lognormal random error of prey populations relB.1 = relB.1 , #Number of prey 1 relative to prey 2; because treating each prey as 1g (for all prey), this is equivalent to relative biomass relB.3 = relB.3, #Number of prey 3 relative to prey 2; because treating each prey as 1g (for all prey), this is equivalent to relative biomass atck.max.1 =atck.max.1, atck.rho.1= atck.rho.1, atck.infl.1= atck.infl.1, #Parameters for sigmoidal Attack rate for prey 1 (max value, rate value rho (neg for decr fxn), and inflection point (age at 50%max)) atck.max.2a=atck.max.2a, #Parameters for double-logistic attack rate for prey 2 (max.2a = scaling parameter for max value) atck.rho.2a=atck.rho.2a, atck.rho.2b = atck.rho.2b, #Parameters for double-logistic attack rate for prey 2 (rho = rate parameter for ascending (2a) and descending (2b) limbs) atck.infl.2a=atck.infl.2a, atck.infl.2b=atck.infl.2b, #Parameters for double-logistic attack rate for prey 2 (infl = age of inflection point for ascending (2a) and descending (2b) limbs) (see Buchheister disseration Chapter 3) atck.max.3 =atck.max.3, atck.rho.3= atck.rho.3, atck.infl.3= atck.infl.3 #Parameters for sigmoidal Attack rate for prey 3 (max value, rate value rho (neg for decr fxn), and inflection point (age at 50%max)) ) S.all[sc,i,,2] = S[,s.spawn] Y.all[sc,i,,2] = Y.yr R.all[sc,i,,2] = R w.all[sc,i,,,,2] = w.pred[,,1,] pulses.all[i,] = pulse.years if(sc==1 && scen.names[sc,"P.scen"] == "Pulse1") N.all[i,,,1,2] = N.prey[,,1,1] #Store Prey Biomass time series for all sims when it is Pulse1 #Note: Note that the prey biomasses are the same for all F and I scenarios if(sc==2 && scen.names[sc,"P.scen"] == "Pulse2") N.all[i,,,2,2] = N.prey[,,1,2] #Store Prey Biomass time series for all sims when it is Pulse2 if(sc==3 && scen.names[sc,"P.scen"] == "Pulse3") N.all[i,,,3,2] = N.prey[,,1,3] #Store Prey Biomass time series for all sims when it is Pulse3 print(paste("scenario =",sc,"of",scenarios,"; sim =",i,"of",sims)) } #end sims loop #w.all[,3,26:50,2,c(1,7),2] ## CALCULATE OUTPUTS ## SSB OUTPUTS #Calculate mean annual difference in pulse and nonpulse SSB for years 26:50 S.diff = round((S.all[sc,,,2] - S.all[sc,,,1])/S.all[sc,,,1]*100,10) #%diff for each year and sim. 100*(Pulse SSB - nonPulse SSB) / (nonPulse SSB) S.diff.mean[,sc] = round(rowMeans(S.diff[,26:50],na.rm=T),10) #take the average %diff value across years in which a pulse could have occurred #Calculate mean annual difference in pulse and nonpulse SSB for PULSE YEARS ONLY (".p") S.diff.p = S.diff[,1:n.pulse] #Create data set for only the pulse years for (i in 1:sims) { #looping across the simulations to calculate %diff in SSB, but only in the years in which there was a pulse event S.diff.p[i,] = round((S.all[sc,i,pulses.all[i,],2] - S.all[sc,i,pulses.all[i,],1])/S.all[sc,i,pulses.all[i,],1]*100,10) } S.diff.p.mean[,sc] = round(rowMeans(S.diff.p[,],na.rm=T),10) #take the average %diff value across years ONLY for years with a pulse ## YIELD OUTPUTS #Calculate mean annual difference in pulse and nonpulse YIELD for years 26:50 Y.diff = round((Y.all[sc,,,2] - Y.all[sc,,,1])/Y.all[sc,,,1]*100,10) #%diff for each year and sim. 100*(Pulse YIELD - nonPulse YIELD) / (nonPulse YIELD) Y.diff.mean[,sc] = round(rowMeans(Y.diff[,26:50],na.rm=T),10) #take the average %diff value across years in which a pulse could have occurred #Calculate mean annual difference in pulse and nonpulse YIELD for PULSE YEARS ONLY (".p") Y.diff.p = Y.diff[,1:n.pulse] #Create data set for only the pulse years for (i in 1:sims) { #looping across the simulations to calculate %diff in YIELD, but only in the years in which there was a pulse event Y.diff.p[i,] = round((Y.all[sc,i,pulses.all[i,],2] - Y.all[sc,i,pulses.all[i,],1])/Y.all[sc,i,pulses.all[i,],1]*100,10) } Y.diff.p.mean[,sc] = round(rowMeans(Y.diff.p[,],na.rm=T),10) #take the average %diff value across years ONLY for years with a pulse ## RECRUITMENT OUTPUTS #Calculate mean annual difference in pulse and nonpulse RECRUITMENT for years 26:50 R.diff = round((R.all[sc,,,2] - R.all[sc,,,1])/R.all[sc,,,1]*100,10) #%diff for each year and sim. 100*(Pulse R - nonPulse R) / (nonPulse R) R.diff.mean[,sc] = round(rowMeans(R.diff[,26:50],na.rm=T),10) #take the average %diff value across years in which a pulse could have occurred #Calculate mean annual difference in pulse and nonpulse RECRUITMENT for PULSE YEARS ONLY (".p") - But note that the comparisons are actually for the year AFTER each pulse event, because i'm interested in the R boost that follows the pulse! R.diff.p = R.diff[,1:n.pulse] #Create data set for only the pulse years for (i in 1:sims) { #looping across the simulations to calculate %diff in RECRUITMENT, but only in the years in which there was a pulse event R.diff.p[i,] = round((R.all[sc,i,(pulses.all[i,]+1),2] - R.all[sc,i,(pulses.all[i,]+1),1])/R.all[sc,i,(pulses.all[i,]+1),1]*100,10) #Note that 1 is added to pulses.all to look at the following year's recruitment } R.diff.p.mean[,sc] = round(rowMeans(R.diff.p[,],na.rm=T),10) #take the average %diff value across years ONLY for years with a pulse ## WEIGHT-AT-AGE #Calculate mean annual difference in pulse and nonpulse WEIGHT@AGE for years 26:50 (note, using season 4 weights at age!) #w.all[scenarios,sims,years,seasons,ages,NP/P] #w.diff[sims,n.pulse,ages] or w.diff[sims,years,ages] #w.diff.mean[sims,ages,scen] w.diff = array(NA,dim=dim(w.all[sc,,,4,,2]), dimnames=list(1:sims,1:years,paste("age",0:(ages-1),sep=""))) #%diff for scenario sc, all sims, all years, season 4, all ages, Pulse. w.diff[,,]=round((w.all[sc,,,4,,2] - w.all[sc,,,4,,1]) / w.all[sc,,,4,,1]*100,10) #Note use of season 4 values w.diff.mean[,,sc] = round(apply(w.diff[,26:50,],c(1,3),mean),10) #take the average %diff value across years in which a pulse could have occurred #Calculate mean differences in weight at age for pulse and nonpulse pairs of years (note, using season 4 weights at age!) w.diff.p = w.diff[,1:n.pulse,] #Create new array, but with # years limited to pulse years w.diff.p[,,] = 99 #Change all values to 99 to help mark if there are any problems with code below. for (i in 1:sims) { #looping across the simulations to calculate %diff in wt@age, but only in the years in which there was a pulse event w.diff.p[i,,] = w.diff[i,pulses.all[i,],] } w.diff.p.mean[,,sc] = round(apply(w.diff.p,c(1,3),mean),10) #take the average %diff value across years in which a pulse occurred } #end scenarios loop ## TIME TO TARGET SSB #Derive a target SSB value, based on a PredPreySim run without any stochastic noise years=80 #Increase num of years to allow population to equilibrate PredPreySim(R.dev.sd=0,N.prey.dev=0, F.scenario="decrF",graph=F) #Run simulation without stochasticity S.target = max(S[,s.spawn]/1E6)*1; S.target #Target SSB defined as 100% of the 80 year equilibrium value. S.target=112.2 kmt (95% of equilibrium is 106.6) years=years.2 #Re-set years to be the value used for the simulation. #Calculate the number of years needed to reach the target SSB ("time (to) target" = time.target) time.target = S.all[,,1,] #Create array to be filled in time.target[,,] = 0 #Recode values to zero to help avoid errors for (sc in 1:dim(time.target)[1]) { #Loop for calculating the # years to reach target SSB for (p in 1:2) { for (s in 1:dim(time.target)[2]) { time.target[sc,s,p] = min(which(S.all[sc,s,,p]/1E6>S.target)) - 25 #extracting the first (ie lowest) year in which SSB was > than the target and subtract the burn in period }}} #NOTE: the previous line will generate a warning if the target SSB is never reached in the sim. (e.g. if F is not decreasing). time.target[,,] = ifelse(time.target[,,]<=0 | time.target[,,] == Inf, NA, time.target[,,]) #Recode nonsensical or uninformative values as NA. #Convert array to data frame for easier plotting time.target.df = melt(time.target, varnames=c("Scenario","Sim","PvNP"), value.name="Time.target") #melt array to facilitate plotting by group time.target.df = merge(time.target.df,scen.names[,2:5], by.x="Scenario", by.y="scen.name") #Bring in some ancillary information ## MAKING PLOTS if(graph==T) { #par(mfrow=c(2,2) #% Increase in SSB plot plot.mar = c(8,5,4,3) windows() par(mar=plot.mar) boxplot(S.diff.mean, main="SSB change due to prey pulses - all years", cex.lab=1.4,xlab="",ylab="mean SSB %diff", notch=notch.gr, las=2) windows() par(mar=plot.mar) boxplot(S.diff.p.mean, main="SSB change due to prey pulses - Pulse years", cex.lab=1.4,xlab="",ylab="mean SSB %diff", notch=notch.gr, las=2) #%Increase (%) in YIELD plot windows() par(mar=plot.mar) boxplot(Y.diff.mean, main="YIELD change due to prey pulses - all years", cex.lab=1.4,xlab="",ylab="mean %diff in YIELD", notch=notch.gr, las=2) windows() par(mar=plot.mar) boxplot(Y.diff.p.mean, main="YIELD change due to prey pulses - Pulse years", cex.lab=1.4,xlab="",ylab="mean %diff in YIELD", notch=notch.gr, las=2) #%Increase (%) in RECRUITMENT plot windows() par(mar=plot.mar) boxplot(R.diff.mean, main="RECRUITMENT change due to prey pulses - all years", cex.lab=1.4,xlab="",ylab="mean %diff in RECRUITMENT", notch=notch.gr, las=2) windows() par(mar=plot.mar) boxplot(R.diff.p.mean, main="RECRUITMENT change due to prey pulses - Pulse years", cex.lab=1.4,xlab="",ylab="mean %diff in RECRUITMENT", notch=notch.gr, las=2) #Change in WEIGHT-AT-AGE if(preys==2) { windows() boxplot(w.diff.p.mean[,,1], main="Weight-at-age change due to prey pulses",ylab="mean %diff",ylim=c(0,max(w.diff.p.mean))) par(new=T) boxplot(w.diff.p.mean[,,2], main="Weight-at-age change due to prey pulses",ylab="mean %diff",col="#FFB6C150",border="red",ylim=c(0,max(w.diff.p.mean))) legend("topright",c("Prey 1 Pulse","Prey 2 Pulse"), col=c("black","red"),lty=1,lwd=3,merge=T,bg="white") } if(preys==3) { library(reshape2) w.diff.p.mean.df = melt(w.diff.p.mean, varnames=c("Sim","Age","Scenario"), value.name="w.diff.p.mean") #melt array to facilitate plotting by group for (j in unique(substr(scen.names$P.scen,1,5))) { #Loop for the different Prey scenarios (Pulse and step) for (i in unique(paste(scen.names$I.scen,scen.names$F.scen[1],sep="-"))) { #Loop for the different migration scenarios w.diff.p.mean.df.1 = subset(w.diff.p.mean.df, substr(Scenario,1,5)==j) #subset for the Prey scenario w.diff.p.mean.df.2 = subset(w.diff.p.mean.df.1, Scenario %in% paste(unique(scen.names$P.scen),i,sep="-")) #subset for each migration scenario w.diff.p.mean.df.2$Scenario = factor(w.diff.p.mean.df.2$Scenario) #re-"factorize" Scenario to remove any unused factor levels windows() boxplot(w.diff.p.mean~Scenario*Age, data=w.diff.p.mean.df.2, main=paste(j," scenarios (",i,")",sep=""), col=c("gray","red","blue"), ylab="Weight-at-age increase (%) during prey pulse", xlab="Age", xaxt="n", notch=notch.gr ) axis(1,at=c(1:(preys*ages)),labels=rep(0:(ages-1),each=preys)) legend("topright",c("Prey 1","Prey 2","Prey 3"), title="Increase in:", col=c("gray","red","blue"),pch=15,bg="white") } #end i loop } #end j loop } #End preys==3 #Weight at age plot #w.all[scen,sim,years,seas,ages,PulseNoPulse] #Make vector of fractional ages x.ages = seq(0.25,ages,by=0.25) #Create an array with only the weight data for the years in which a pulse occurred w.pulses = w.all[,,c(1:n.pulse),,,2] # make a smaller array with only n.pulse # of years w.pulses[,,,,] =0 #turn all entries into zero to help avoid errors for (s in 1:sims){ #fill in the array with data from PULSE years only. w.pulses[,s,,,] = w.all[,s,pulses.all[s,],,,2] } #Make Plot - PULSE YEARS #w.all[scenarios,sims,years,seasons,ages,NP/P] windows() plot(LWdata$AGE2,LWdata$WEIGHT_KG, col="darkgray", main=paste("Mean weights at age (Pulse years)",scen.names$I.scen[1],sep=" - "), ylab="Weight (kg)",xlab="Age",xlim=c(0,12),ylim=c(0,4)) lines(x.ages, apply(w.all[1,,,,,1],c(3,4),mean), col="black", lwd=3, lty=2) #Scen=1 (NO PULSE); averaging across years and simulations. NOTE: The NO-pulse values are identical across Scenarios as long as the I.scen is the same. lines(x.ages, apply(w.pulses[1,,,,],c(3,4),mean), col="black", lwd=6) #Scen=Pulse1 (PULSE YEARS); averaging across years and simulations lines(x.ages, apply(w.pulses[2,,,,],c(3,4),mean), col="red", lwd=3) #Scen=Pulse2 (PULSE YEARS); averaging across years and simulations lines(x.ages, apply(w.pulses[3,,,,],c(3,4),mean), col="blue", lwd=5) #Scen=Pulse3 (PULSE YEARS); averaging across years and simulations legend("topleft",c("No Pulse",dimnames(w.pulses)[[1]][1:3]), col=c("black","black","red","blue"), lty=c(2,1,1,1), lwd=c(2,6,3,5)) #Make Plot - ALL YEARS (for Step scenarios if present) if("Step." %in% substr(scen.names$P.scen,1,5)) { windows() plot(LWdata$AGE2,LWdata$WEIGHT_KG, col="darkgray", main=paste("Mean weights at age (All years)",scen.names$I.scen[1],sep=" - "), ylab="Weight (kg)",xlab="Age",xlim=c(0,12),ylim=c(0,4)) lines(x.ages, apply(w.all[1,,,,,1],c(3,4),mean), col="black", lwd=3, lty=2) #Scen=1 (NO PULSE); averaging across years and simulations. NOTE: The NO-pulse values are identical across Scenarios as long as the I.scen is the same. lines(x.ages, apply(w.all[4,,26:50,,,2],c(3,4),mean), col="black", lwd=6) #Scen=Step.1 (All years); averaging across years and simulations. lines(x.ages, apply(w.all[5,,26:50,,,2],c(3,4),mean), col="red", lwd=3) #Scen=Step.2 (All years); averaging across years and simulations. lines(x.ages, apply(w.all[6,,26:50,,,2],c(3,4),mean), col="blue", lwd=5) #Scen=Step.3 (All years); averaging across years and simulations. legend("topleft",c("No Pulse",dimnames(w.pulses)[[1]][4:6]), col=c("black","black","red","blue"), lty=c(2,1,1,1), lwd=c(2,6,3,5)) } #Time to SSB target plots #Make boxplot of time to target time.target.df$Scen = paste(time.target.df$Scenario,time.target.df$PvNP,sep="-") #Create new variable to use to subset the data time.target.sub = rbind(subset(time.target.df, Scen == "Pulse1-mig-decrF-NoPulse"), #Subset to include 1 NoPulse scenario (they're all the same) subset(time.target.df, PvNP == "Pulse" & substr(time.target.df$Scenario,8,16) == "mig-decrF")) #And also only the decreasing F scenarios with migration. time.target.sub$Scen = factor(time.target.sub$Scen ) #Re-factorize variable to get rid of unwanted factor levels. windows() if(nrow(time.target.sub)>0){ #only generate time to SSB plot if there is data to plot boxplot(Time.target~Scen, data=time.target.sub, main="Pulse Effects on Management Goals", #col=c("gray","red","blue"), ylab="Time to reach target SSB (yr)", xlab="Prey Scenario", xaxt="n", las=2, notch=notch.gr, cex.lab=1.2 ) axis(1,at=c(1:length(levels(time.target.sub$Scen))),labels=c("None", substr(levels(time.target.sub$Scen),1,6)[-1])) minor.tick(ny=5, tick.ratio=0.5) } }#END graph==T section #Store tables if(TRUE) { #COMMENT THIS OUT, use the list output below instead S.all <<- S.all S.diff.mean <<- S.diff.mean S.diff.p.mean <<- S.diff.p.mean Y.all <<- Y.all Y.diff.mean <<- Y.diff.mean Y.diff.p.mean <<- Y.diff.p.mean R.all <<- R.all R.diff.mean <<- R.diff.mean R.diff.p.mean <<- R.diff.p.mean w.all <<- w.all w.diff.mean <<- w.diff.mean w.diff.p.mean <<- w.diff.p.mean time.target <<- time.target #needed for sensitivity function time.target.df <<- time.target.df N.all <<- N.all pulses.all <<- pulses.all } #Store function call with parameter values: function.call = list( scen.names = scen.names, predators = predators, ages = ages, preys = preys, years = years, years.2 = years.2, seasons = seasons, regions = regions, seed.init=seed.init, sims=sims, #Number of simulation runs to conduct graph=graph, #Switch for turning graphs on or off, notch.gr=notch.gr, #Indicate whether boxplots have a notch #PredPreySim parameters: M=M, #Instantaneous natural mortality rate (per year) Fmax = Fmax, #Maximum instantaneous fishing mortality. Previously 2.0; exploring how 1.5 looks since that was more the average. (See StockRecruitFxn_8-29-13.xlsx) Fmin = Fmin, #Minimum instantaneous fishing mortality. decrF.time = decrF.time, #Number of years to go from Fmax to Fmin F.sel = F.sel, #Fishing selectivity by age R.dev.sd=R.dev.sd, #Standard deviation for the lognormal random recruitment error alpha.SR = alpha.SR, #alpha parameter for Shepherd S-R function (Rothschild et al. 2012) K.SR = K.SR, #K parameter for Shepherd S-R function (Rothschild et al. 2012) beta.SR = beta.SR, #beta shape parameter for Shepherd S-R function. beta.SR=1 is for Beverton Holt. CTO= CTO , #Temperature for optimum consumption; assumed value (see eq in Hanson et al. 1997) CTM= CTM , #Maximum temp parameter for temperature equation; assumed value (see eq in Hanson et al. 1997) CQ = CQ , #Rate parameter for temperature equation; assumed value (see eq in Hanson et al. 1997) CA = CA , #Intercept Parameter for allometric relationship of Cmax (Hanson et al. 1997; GrowthModel_ConvEfficiency_8-25-13.xlsx) Cmax = CA*W^CB CB = CB , #Decay Parameter for allometric relationship of Cmax (Hanson et al. 1997; GrowthModel_ConvEfficiency_8-25-13.xlsx) Cmax = CA*W^CB KL = KL , #Lambda parameter for maximum conversion effiency (K) sigmoidal model(GrowthModel_ConvEfficiency_8-25-13.xlsx) KR = KR , #Rate of decline for sigmoidal K model (GrowthModel_ConvEfficiency_8-25-13.xlsx) KW = KW , #Weight at which KL/2 is attained in Sigmoidal K model (GrowthModel_ConvEfficiency_8-25-13.xlsx) #Gmax.a = Gmax.a, #Age-specific maximum growth (kg/yr) possible for species n.pulse=n.pulse, #Number of prey pulses pulse.range = pulse.range, #Range of multiplier values for when a pulse in production occurs. E.g. 2-10x production in pulse years N.prey.2.mean = N.prey.2.mean , #Number of prey 2 N.prey.dev = N.prey.dev, #Standard deviation for the lognormal random error of prey populations relB.1 = relB.1 , #Number of prey 1 relative to prey 2; because treating each prey as 1g (for all prey), this is equivalent to relative biomass relB.3 = relB.3, #Number of prey 3 relative to prey 2; because treating each prey as 1g (for all prey), this is equivalent to relative biomass atck.max.1 =atck.max.1, atck.rho.1= atck.rho.1, atck.infl.1= atck.infl.1, #Parameters for sigmoidal Attack rate for prey 1 (max value, rate value rho (neg for decr fxn), and inflection point (age at 50%max)) atck.max.2a=atck.max.2a, #Parameters for double-logistic attack rate for prey 2 (max.2a = scaling parameter for max value) atck.rho.2a=atck.rho.2a, atck.rho.2b = atck.rho.2b, #Parameters for double-logistic attack rate for prey 2 (rho = rate parameter for ascending (2a) and descending (2b) limbs) atck.infl.2a=atck.infl.2a, atck.infl.2b=atck.infl.2b, #Parameters for double-logistic attack rate for prey 2 (infl = age of inflection point for ascending (2a) and descending (2b) limbs) (see Buchheister disseration Chapter 3) atck.max.3 =atck.max.3, atck.rho.3=atck.rho.3, atck.infl.3=atck.infl.3 #Parameters for sigmoidal Attack rate for prey 3 atck.max.1 =atck.max.1, atck.max.2=atck.max.2, atck.max.3=atck.max.3, #Max attack rate on prey 1 and 2 ) #Return results in a list sim.out = list(function.call = function.call, S.all = S.all, S.diff.mean = S.diff.mean, S.diff.p.mean=S.diff.p.mean, Y.all = Y.all, Y.diff.mean = Y.diff.mean, Y.diff.p.mean=Y.diff.p.mean, R.all = R.all, R.diff.mean = R.diff.mean, R.diff.p.mean=R.diff.p.mean, w.all = w.all, w.diff.mean=w.diff.mean, w.diff.p.mean=w.diff.p.mean, time.target.df = time.target.df, time.target = time.target, N.all = N.all, pulses.all=pulses.all ) return(sim.out) } #END simfxn() ##NOTES ON ERRORS/WARNINGS AND RUNNING THE simfxn(): # - warnings about the time.target calculation are ok. They just indicate situations where the target SSB was not reached. # - Plots of %incr in W@age will be screwy if you don't have pulses for all three prey (the colors will be off and non-informative. #JUNK - Working on coding for Step scenarios test = simfxn(sims=5, graph=T) #started 1140 pm - test$function.call time.target.df #TROUBLESHOOTING/EXPLORATION (9-9-13) - Figuring out why sim#12,13 "NoPulse" are anomalous for time to target - yielding a large effect of "Pulses" on time to target SSB. #See manuscript outline> results> time to target, for notes. #Run Simulation with 14 runs to see sims 12 and 13 test = simfxn(sims=14, graph=T) test$S.all[1,,,"NoPulse"] test$S.all[1,1,,] #Plot the SSB trajectories, and highlight some specific cases library(Hmisc) plot(1:years, test$S.all[1,12,,"NoPulse"], col="red", type="l", xlab="Year", ylab="SSB (kg)", main = "SSB trajectories for NoPulse, decrF, mig scenario") abline(v=25) abline(v=58) abline(h=1.085e+08, lty=2) #This is the S.target minor.tick(nx=4) for (i in 1:14) { lines(test$S.all[1,i,,"NoPulse"], col="black", type="l") } lines(test$S.all[1,12,,"NoPulse"], col="red", type="l", lwd=3) #Took 33 yrs to reach target (sim=12) lines(test$S.all[1,13,,"NoPulse"], col="blue", type="l", lwd=3) #Took 34 yrs to reach target (sim=13) lines(test$S.all[1,9,,"NoPulse"], col="green", type="l", lwd=3) #Never reached the target (sim=9) #JUNK 2-10-13 - exploring the different decrF subscenarios. #Time to SSB target plots #Make boxplot of time to target time.target.df$Scen = paste(time.target.df$Scenario,time.target.df$PvNP,sep="-") #Create new variable to use to subset the data time.target.sub = time.target.df[-which(time.target.df[,"PvNP"]=="NoPulse" & time.target.df[,"P.scen"] !="Pulse1"),] #time.target.sub$Scen = factor(time.target.sub$Scen, levels = # c("Pulse1-mig-decrF.24-NoPulse","Pulse1-mig-decrF.22-NoPulse", "Pulse1-mig-decrF.21-NoPulse","Pulse1-mig-decrF-NoPulse","Pulse1-mig-decrF.19-NoPulse","Pulse1-mig-decrF.18-NoPulse","Pulse1-mig-decrF.16-NoPulse", # "Pulse1-mig-decrF.24-Pulse","Pulse1-mig-decrF.22-Pulse", "Pulse1-mig-decrF.21-Pulse","Pulse1-mig-decrF-Pulse","Pulse1-mig-decrF.19-Pulse","Pulse1-mig-decrF.18-Pulse","Pulse1-mig-decrF.16-Pulse", # "Pulse2-mig-decrF.24-Pulse","Pulse2-mig-decrF.22-Pulse", "Pulse2-mig-decrF.21-Pulse","Pulse2-mig-decrF-Pulse","Pulse2-mig-decrF.19-Pulse","Pulse2-mig-decrF.18-Pulse","Pulse2-mig-decrF.16-Pulse", # "Pulse3-mig-decrF.24-Pulse","Pulse3-mig-decrF.22-Pulse", "Pulse3-mig-decrF.21-Pulse","Pulse3-mig-decrF-Pulse","Pulse3-mig-decrF.19-Pulse","Pulse3-mig-decrF.18-Pulse","Pulse3-mig-decrF.16-Pulse") # ) time.target.sub$Scen = factor(time.target.sub$Scen, levels = c("Pulse1-mig-decrF.14-NoPulse", "Pulse1-mig-decrF.14-Pulse","Pulse2-mig-decrF.14-Pulse", "Pulse3-mig-decrF.14-Pulse", "Pulse1-mig-decrF.12-NoPulse", "Pulse1-mig-decrF.12-Pulse","Pulse2-mig-decrF.12-Pulse", "Pulse3-mig-decrF.12-Pulse", "Pulse1-mig-decrF.10-NoPulse", "Pulse1-mig-decrF.10-Pulse","Pulse2-mig-decrF.10-Pulse", "Pulse3-mig-decrF.10-Pulse") ) #time.target.sub$Scen = factor(time.target.sub$Scen ) #Re-factorize variable to get rid of unwanted factor levels. windows() par(mfrow=c(2,2)) par(mar=c(14,5,5,2)) boxplot(Time.target~Scen, data=time.target.sub, main="Pulse Effects on Management Goals", #col=c("gray","red","blue"), ylab="Time to reach target SSB (yr)", xlab="", las=2, notch=notch.gr, cex.lab=1.2, ylim=c(10,44) ) # xaxt="n", #axis(1,at=c(1:4),labels=c("None", "Prey 1", "Prey 2", "Prey 3")) minor.tick(ny=5, nx=0, tick.ratio=0.5) abline(v=c(4.5,8.5,12.5,16.5,20.5,24.5,28.5,32.5), lty=2, col="red") summary(time.target.sub) #Timing of the function system.time(simfxn(sims=3)) #elapsed: 355s=6min. Would take ~33hrs to run the 1000 sims. ### RUNNING THE SIMULATION - NOTE: BE SURE TO DEFINE THE scen.names APPROPRIATELY! (located before simfxn() code) #Main run #Started 3:00pm - End pm 2-11-13 SIM.MOD = simfxn(sims=100) summary(SIM.MOD) names(SIM.MOD) SIM.MOD$function.call dim(SIM.MOD$S.all); length(SIM.MOD$S.all) #Save simulation model output. #save(SIM.MOD, file="C:\\Users\\andrebuc\\VIMS\\PhD\\Research\\PredPrey_Simulations\\SimOutputs\\SimModelOutput_2-11-13.RData") #Run for decreasing F subscenarios (2-11-13) #Started 12:40am - End ?? years=80 years.2 = years #This is for a coding nuance... necessary when estimating the target SSB SIM.MOD.F.scen = simfxn(sims=100) summary(SIM.MOD.F.scen) names(SIM.MOD.F.scen) SIM.MOD.F.scen$function.call dim(SIM.MOD.F.scen$S.all); length(SIM.MOD.F.scen$S.all) #Save simulation model output. #save(SIM.MOD.F.scen, file="C:\\Users\\andrebuc\\VIMS\\PhD\\Research\\PredPrey_Simulations\\SimOutputs\\SimModelOutput_F.scen_2-12-13.RData") #Main run - now with Step scenarios added #Started 9:20pm - End 12:20 am 3-18-13 SIM.MOD = simfxn(sims=100) summary(SIM.MOD) names(SIM.MOD) SIM.MOD$function.call dim(SIM.MOD$S.all); length(SIM.MOD$S.all) #Save simulation model output. #save(SIM.MOD, file="C:\\Users\\andrebuc\\VIMS\\PhD\\Research\\PredPrey_Simulations\\SimOutputs\\SimModelOutput_3-18-13.RData") #Run for decreasing F subscenarios (3-18-13) #Started 12:28am - End ?? years=80 years.2 = years #This is for a coding nuance... necessary when estimating the target SSB SIM.MOD.F.scen = simfxn(sims=60) summary(SIM.MOD.F.scen) names(SIM.MOD.F.scen) SIM.MOD.F.scen$function.call dim(SIM.MOD.F.scen$S.all); length(SIM.MOD.F.scen$S.all) #Save simulation model output. #save(SIM.MOD.F.scen, file="C:\\Users\\andrebuc\\VIMS\\PhD\\Research\\PredPrey_Simulations\\SimOutputs\\SimModelOutput_F.scen_3-18-13.RData") #NEW Main run - now with Step scenarios added, and decrF defined as 12 years from 1.5 to 0.255. #Started 9:31am - End ?? am 9-6-13 memory.limit(2000) #Increasing a bit just in case. SIM.MOD = simfxn(sims=100) #save(SIM.MOD, file="C:\\Users\\andrebuc\\VIMS\\PhD\\Research\\PredPrey_Simulations\\SimOutputs\\SimModelOutput_9-6-13.RData") #run - Adding in a variety of different F.scenarios (e.g. mor.F, ass.F, exp.F); also, decrF defined as 12 years from 1.5 to 0.31. #Started 9:31am - End ?? am 9-19-13 memory.limit(2000) #Increasing a bit just in case. SIM.MOD = simfxn(sims=100) #save(SIM.MOD, file="C:\\Users\\andrebuc\\VIMS\\PhD\\Research\\PredPrey_Simulations\\SimOutputs\\SimModelOutput_9-19-13.RData") #RUN - Changed pulse.range to be c(2:6); P.scen = Pulse1,2,3; F.scen= hi,lo,decr,dec2,mor; I.scen=mig,mix; years=80 #Started 1230am - 9-24-13 memory.limit(2000) #Increasing a bit just in case. SIM.MOD = simfxn(sims=100) #save(SIM.MOD, file="C:\\Users\\andrebuc\\VIMS\\PhD\\Research\\PredPrey_Simulations\\SimOutputs\\SimModelOutput_9-24-13.RData") #RUN - Same as 9/24/13 above, except fixed a very minor error in simfxn() output for "all years" calculation. #Started 1230am - 9-25-13 memory.limit(2000) #Increasing a bit just in case. SIM.MOD = simfxn(sims=100) #save(SIM.MOD, file="C:\\Users\\andrebuc\\VIMS\\PhD\\Research\\PredPrey_Simulations\\SimOutputs\\SimModelOutput_9-25-13.RData") #RUN - Same as 9/25/13 above, except i am now also outputting/saving data on the recruitments by year. #Started 11:59 pm - 9-30-13 memory.limit(2000) #Increasing a bit just in case. SIM.MOD = simfxn(sims=100) #save(SIM.MOD, file="C:\\Users\\andrebuc\\VIMS\\PhD\\Research\\PredPrey_Simulations\\SimOutputs\\SimModelOutput_9-30-13.RData") summary(SIM.MOD) names(SIM.MOD) SIM.MOD$function.call dim(SIM.MOD$S.all); length(SIM.MOD$S.all) #RUN - Running 1000 simulations #Started pm - 7-25-14 memory.limit(3500) #Increasing a bit just in case. SIM.MOD = simfxn(sims=600, seed.init=100) #Rest of the simulations were run on the DELL computer #save(SIM.MOD, file="C:\\Users\\andre\\VIMS\\PhD\\Research\\PredPrey_Simulations\\SimOutputs\\SimModelOutput_600_7-25-14.RData") summary(SIM.MOD) names(SIM.MOD) SIM.MOD$function.call dim(SIM.MOD$S.all); length(SIM.MOD$S.all) ############################################################################################################ ############################################################################################################ ### SENSITIVITY FUNCTION ############################################################################################################ ############################################################################################################ #Description: # 2 different sensitivity approaches were explored: # 1) "Monte Carlo-type" senesitivity - All parameters were randomly selected from a uniform distribution # bounded by +/- 20%. New random parameters were chosen after each simfxn() run. # GOAL: look at overall output sensitivity to overall parameter uncertainty # 2) Individual parameter perturbations - the simulation function was run with all of the default # parameter values, except that a single parameter was either increased or decreased by 20%. # GOAL: To evaluate the relative sensitivity of the model results to each of the model parameters (and model # components: e.g. growth, mortality, prey, etc.), # # NOTE: Need to run the PredPreySim() and simfxn() from above for this code to work. #DEFINE GLOBAL SCENARIOS #Define the number of global scenarios to be investigated (e.g. prey 1 pulses, prey 2 pulses) ,"Step.1","Step.2","Step.3" scen.names = expand.grid(P.scen=c("Pulse1","Pulse2","Pulse3"),I.scen=c("mig","mix"),F.scen=c("highF","lowF","decrF","dec2F","mor.F")) # scen.names$scenario = c(1:nrow(scen.names)) scen.names$scen.name = paste(scen.names$P.scen, scen.names$I.scen, scen.names$F.scen, sep="-") scen.names = subset(scen.names, select=c("scenario","scen.name","P.scen","I.scen","F.scen")) #Re-order columns to match the previous version scenarios = nrow(scen.names) #For Troubleshooting (OTHERWISE: run the same scen.names as before for the simfxn(); copied above): if(FALSE){ scen.names = expand.grid(P.scen=c("Pulse1","Pulse2","Pulse3"),I.scen=c("mig"),F.scen=c("decrF", "highF")) # scen.names$scenario = c(1:nrow(scen.names)) scen.names$scen.name = paste(scen.names$P.scen, scen.names$I.scen, scen.names$F.scen, sep="-") scen.names = subset(scen.names, select=c("scenario","scen.name","P.scen","I.scen","F.scen")) #Re-order columns to match the previous version scenarios = nrow(scen.names) } #Create table with names and order of parameter value sensitivity analysis sens.names = expand.grid(parm = c("M", "alpha.SR", "K.SR", "beta.SR", "R.dev.sd", #Mortality and Recruitment "n.pulse", "pulse.range", "N.prey.2.mean", "relB.1", "relB.3", "N.prey.dev", #Prey Biomass and pulses "CA", "CB", "CTO", "CTM", "CQ","KL", "KR", "KW", #Growth submodel parms: "atck.max.1", "atck.max.2a", "atck.max.3", "atck.rho.1", "atck.rho.2a", "atck.rho.3", "atck.infl.1", "atck.infl.2a", "atck.infl.2b", "atck.infl.3"), val = c("lo","hi"), stringsAsFactors=F) sens.names = rbind(c("BASE",""),sens.names) #Add BASE run as a scenario sens.names$name = paste(sens.names[,1],sens.names[,2], sep=".") #Create sensitivity scenario names (e.g. "M.lo", "M.hi") sens.scenarios = nrow(sens.names) #Number of sensitivity scenarios (including BASE run) sens.perc = 0.2 #Fractional value used to change parameter base values (e.g. 0.2 = 20%) ## SPECIFY DEFAULT (def) VALUES FOR simfxn() AND HI/LO VALUES (commented out any non singular value parameters) BASE.def = matrix(c( "M" , 0.25, "R.dev.sd" , 0.4, "alpha.SR" , 3.4, "K.SR" , 27.4, "beta.SR" , 1, "CTO", 22 , "CTM", 35 , "CQ" , 2.5 , "CA", 0.3 , "CB", -0.2 , "KL" , 0.5 , "KR" , -0.0014 , "KW" , -600 , "n.pulse" , 5, "N.prey.2.mean" , 3e+10, "N.prey.dev" , 0.4, "relB.1" , 10, "relB.3" , 0.1, "atck.max.1" , exp(-19.5), "atck.rho.1", -1.1, "atck.infl.1", 2, "atck.max.2a", exp(-16.3), "atck.rho.2a", 0.46, "atck.rho.2b" , 0.46, "atck.infl.2a", 1, "atck.infl.2b", 4, "atck.max.3" , exp(-13.5), "atck.rho.3", 1.9, "atck.infl.3", 3.5), ncol=2,byrow=T) BASE.def = data.frame(Parm=BASE.def[,1],Value=as.numeric(BASE.def[,2])) #Calculate low and high values for each parameter BASE.def$Val.lo = BASE.def$Value*(1-sens.perc) BASE.def$Val.hi = BASE.def$Value*(1+sens.perc) ##DEFINE NUMBER OF SIMS TO BE USED FOR SENSITIVITY ANALYSIS sims.mc = 2 #For Monte-Carlo approach this needs to be >1 otherwise, the simfxn() chokes. Set to 2. ### METHOD 1 - MONTE-CARLO-TYPE SENSITIVITY - PERTURBATION OF ALL PARAMETERS (RANDOMLY PERTURBING ALL PARAMETERS SIMULATEOUSLY) #Specify number of times to run simfxn() with new random draws of parameters sens.runs = 3 #100 ## CREATE EMPTY ARRAYS TO FILL WITH SENSITIVITY RESULTS #Summary of %diff in SSB (S), Yield (Y), and weight at age (w) S.sens = array(NA, dim=c(sims.mc,scenarios,sens.runs),dimnames=list(paste("sim",1:sims.mc),scen.names[,2], paste("sens.run",1:sens.runs))) Y.sens = array(NA, dim=c(sims.mc,scenarios,sens.runs),dimnames=list(paste("sim",1:sims.mc),scen.names[,2], paste("sens.run",1:sens.runs))) R.sens = array(NA, dim=c(sims.mc,scenarios,sens.runs),dimnames=list(paste("sim",1:sims.mc),scen.names[,2], paste("sens.run",1:sens.runs))) w.sens = array(NA, dim=c(sims.mc,ages,scenarios,sens.runs),dimnames=list(paste("sim",1:sims.mc),paste("age-",0:(ages-1),sep=""),scen.names[,2], paste("sens.run",1:sens.runs))) t.sens = array(NA, dim=c(sims.mc,scenarios,2,sens.runs),dimnames=list(paste("sim",1:sims.mc),scen.names[,2], c("NoPulse","Pulse"),paste("sens.run",1:sens.runs))) dim(S.sens); length(S.sens) dim(w.sens); length(w.sens) ## RUN simfxn(), EACH TIME RANDOMLY DRAWING NEW PARAMETERS system.time( for (i in 1:sens.runs) { temp = simfxn(sims=sims.mc, #Have to run simfxn more than once for internal calculation to work... graph=F, seed.init=500*i, notch.gr=F, #Parameters to randomly select (comment out any parameters to remain constant) (Also note that the order "Val.lo" and "Val.hi" for negative parameters has to be reversed for runif() to work. M=runif(1,BASE.def[BASE.def$Parm=="M","Val.lo"], BASE.def[BASE.def$Parm=="M","Val.hi"]), #Instantaneous natural mortality rate (per year) #Fmax = 1.5, #Maximum instantaneous fishing mortality. Previously 2.0; exploring how 1.5 looks since that was more the average. (See StockRecruitFxn_8-29-13.xlsx) #Fmin = 0.31, #Minimum instantaneous fishing mortality. #decrF.time = 12, #Number of years to go from Fmax to Fmin #F.sel = c(0,0.1,0.5,rep(1.0,ages-3)), #Fishing selectivity by age R.dev.sd=runif(1,BASE.def[BASE.def$Parm=="R.dev.sd","Val.lo"], BASE.def[BASE.def$Parm=="R.dev.sd","Val.hi"]), #Standard deviation for the lognormal random recruitment error alpha.SR = runif(1,BASE.def[BASE.def$Parm=="alpha.SR","Val.lo"], BASE.def[BASE.def$Parm=="alpha.SR","Val.hi"]), #alpha parameter for Shepherd S-R function (Rothschild et al. 2012) K.SR = runif(1,BASE.def[BASE.def$Parm=="K.SR","Val.lo"], BASE.def[BASE.def$Parm=="K.SR","Val.hi"]), #K parameter for Shepherd S-R function (Rothschild et al. 2012) beta.SR = runif(1,BASE.def[BASE.def$Parm=="beta.SR","Val.lo"], BASE.def[BASE.def$Parm=="beta.SR","Val.hi"]), #beta shape parameter for Shepherd S-R function. beta.SR=1 is for Beverton Holt. CTO= runif(1,BASE.def[BASE.def$Parm=="CTO","Val.lo"], BASE.def[BASE.def$Parm=="CTO","Val.hi"]) , #Temperature for optimum consumption; assumed value (see eq in Hanson et al. 1997) CTM= runif(1,BASE.def[BASE.def$Parm=="CTM","Val.lo"], BASE.def[BASE.def$Parm=="CTM","Val.hi"]) , #Maximum temp parameter for temperature equation; assumed value (see eq in Hanson et al. 1997) CQ = runif(1,BASE.def[BASE.def$Parm=="CQ","Val.lo"], BASE.def[BASE.def$Parm=="CQ","Val.hi"]) , #Rate parameter for temperature equation; assumed value (see eq in Hanson et al. 1997) CA = runif(1,BASE.def[BASE.def$Parm=="CA","Val.lo"], BASE.def[BASE.def$Parm=="CA","Val.hi"]) , #Intercept Parameter for allometric relationship of Cmax (Hanson et al. 1997; GrowthModel_ConvEfficiency_8-25-13.xlsx) Cmax = CA*W^CB CB = runif(1,BASE.def[BASE.def$Parm=="CB","Val.hi"], BASE.def[BASE.def$Parm=="CB","Val.lo"]) , #Decay Parameter for allometric relationship of Cmax (Hanson et al. 1997; GrowthModel_ConvEfficiency_8-25-13.xlsx) Cmax = CA*W^CB KL = runif(1,BASE.def[BASE.def$Parm=="KL","Val.lo"], BASE.def[BASE.def$Parm=="KL","Val.hi"]) , #Lambda parameter for maximum conversion effiency (K) sigmoidal model(GrowthModel_ConvEfficiency_8-25-13.xlsx) KR = runif(1,BASE.def[BASE.def$Parm=="KR","Val.hi"], BASE.def[BASE.def$Parm=="KR","Val.lo"]) , #Rate of decline for sigmoidal K model (GrowthModel_ConvEfficiency_8-25-13.xlsx) KW = runif(1,BASE.def[BASE.def$Parm=="KW","Val.hi"], BASE.def[BASE.def$Parm=="KW","Val.lo"]) , #Weight at which KL/2 is attained in Sigmoidal K model (GrowthModel_ConvEfficiency_8-25-13.xlsx) n.pulse=runif(1,BASE.def[BASE.def$Parm=="n.pulse","Val.lo"], BASE.def[BASE.def$Parm=="n.pulse","Val.hi"]), #Number of prey pulses pulse.range = matrix(c(2:4,4:6),nrow=2,byrow=T)[sample(c(1,2),1,replace=T),], #This randomly makes pulse.range either c(2:4) or c(4:6) #Range of multiplier values for when a pulse in production occurs. E.g. 2-10x production in pulse years N.prey.2.mean = runif(1,BASE.def[BASE.def$Parm=="N.prey.2.mean","Val.lo"], BASE.def[BASE.def$Parm=="N.prey.2.mean","Val.hi"]) , #Number of prey 2 N.prey.dev = runif(1,BASE.def[BASE.def$Parm=="N.prey.dev","Val.lo"], BASE.def[BASE.def$Parm=="N.prey.dev","Val.hi"]), #Standard deviation for the lognormal random error of prey populations relB.1 = runif(1,BASE.def[BASE.def$Parm=="relB.1","Val.lo"], BASE.def[BASE.def$Parm=="relB.1","Val.hi"]) , #Number of prey 1 relative to prey 2; because treating each prey as 1g (for all prey), this is equivalent to relative biomass relB.3 = runif(1,BASE.def[BASE.def$Parm=="relB.3","Val.lo"], BASE.def[BASE.def$Parm=="relB.3","Val.hi"]), #Number of prey 3 relative to prey 2; because treating each prey as 1g (for all prey), this is equivalent to relative biomass atck.max.1 =runif(1,BASE.def[BASE.def$Parm=="atck.max.1","Val.lo"], BASE.def[BASE.def$Parm=="atck.max.1","Val.hi"]), atck.rho.1= runif(1,BASE.def[BASE.def$Parm=="atck.rho.1","Val.hi"], BASE.def[BASE.def$Parm=="atck.rho.1","Val.lo"]), atck.infl.1= runif(1,BASE.def[BASE.def$Parm=="atck.infl.1","Val.lo"], BASE.def[BASE.def$Parm=="atck.infl.1","Val.hi"]), #Parameters for sigmoidal Attack rate for prey 1 (max value, rate value rho (neg for decr fxn), and inflection point (age at 50%max)) atck.max.2a=runif(1,BASE.def[BASE.def$Parm=="atck.max.2a","Val.lo"], BASE.def[BASE.def$Parm=="atck.max.2a","Val.hi"]), #Parameters for double-logistic attack rate for prey 2 (max.2a = scaling parameter for max value) atck.rho.2a=runif(1,BASE.def[BASE.def$Parm=="atck.rho.2a","Val.lo"], BASE.def[BASE.def$Parm=="atck.rho.2a","Val.hi"]), atck.rho.2b = runif(1,BASE.def[BASE.def$Parm=="atck.rho.2b","Val.lo"], BASE.def[BASE.def$Parm=="atck.rho.2b","Val.hi"]), #Parameters for double-logistic attack rate for prey 2 (rho = rate parameter for ascending (2a) and descending (2b) limbs) atck.infl.2a=runif(1,BASE.def[BASE.def$Parm=="atck.infl.2a","Val.lo"], BASE.def[BASE.def$Parm=="atck.infl.2a","Val.hi"]), atck.infl.2b=runif(1,BASE.def[BASE.def$Parm=="atck.infl.2b","Val.lo"], BASE.def[BASE.def$Parm=="atck.infl.2b","Val.hi"]), #Parameters for double-logistic attack rate for prey 2 (infl = age of inflection point for ascending (2a) and descending (2b) limbs) (see Buchheister disseration Chapter 3) atck.max.3 =runif(1,BASE.def[BASE.def$Parm=="atck.max.3","Val.lo"], BASE.def[BASE.def$Parm=="atck.max.3","Val.hi"]), atck.rho.3= runif(1,BASE.def[BASE.def$Parm=="atck.rho.3","Val.lo"], BASE.def[BASE.def$Parm=="atck.rho.3","Val.hi"]), atck.infl.3= runif(1,BASE.def[BASE.def$Parm=="atck.infl.3","Val.lo"], BASE.def[BASE.def$Parm=="atck.infl.3","Val.hi"]) #Parameters for sigmoidal Attack rate for prey 3 (max value, rate value rho (neg for decr fxn), and inflection point (age at 50%max)) ) S.sens[,,i] = temp$S.diff.p.mean Y.sens[,,i] = temp$Y.diff.p.mean R.sens[,,i] = temp$R.diff.p.mean w.sens[,,,i] = temp$w.diff.p.mean t.sens[,,,i] = aperm(temp$time.target,c(2,1,3)) print(paste("Sensitivity run",i,"of",sens.runs)) } ) #Timing notes (7-27-14) #~2min 15sec (136s) for 6 scen, 3 sens.runs, = 18 rounds #NEED: 30 scen, 1000 sens.runs = 30,000 rounds x 136s = 4.08E6 s = 47 days! #Store sensitivity results data objects for future use and reference (10-6-13) 2am-? #summary(S.sens); summary(Y.sens); summary(w.sens) #Check that there are no NA's #Save summary tables SSB, Yield, Weight at age, and time to target save(S.sens, Y.sens, R.sens, w.sens, t.sens, file="C:\\Users\\andrebuc\\VIMS\\PhD\\Research\\PredPrey_Simulations\\SimOutputs\\SensitivityOutput_MC_10-5-13.RData") #backup.S = S.sens; backup.Y = Y.sens; backup.w = w.sens summary(S.sens); summary(R.sens); summary(w.sens); summary(t.sens) #JUNK TROUBLESHOOTING # S.sens[sims.mc,scenarios,sens.runs]; # w.sens[sims.mc,ages,scenarios,sens.runs] # t.sens[sims.mc,scenarios,2,sens.runs] boxplot(t(S.sens[1,,]),las=2 ) t(S.sens[1,,]) #Identical values!! S.sens[,,] S.sens[,c(1,7),1] ### METHOD 2 - SENSITIVITY TO INDIVIDUAL PARAMETER PERTURBATION (MOVING EACH PARAMETER UP OR DOWN 20%) # Re-set the number of simulations to be run sims=50 #50 #DEFINE GLOBAL SCENARIOS # - Originally, there were 36 scenarios, 100 sims per scenario, and 30 parameters to perturb (x2 for hi AND low values). # At ~8min for 100 sims of 1 scenario, that would take 288 hours (12 days)! # - So i chose to exclude "step", "mix", "highF" and "lowF" scenarios. Yielding 3 scen * 60 model runs x 4.5min (for 50 sims) ~ 14 hours # - I also ran the script in two parts (one with tinnR, the other with RStudio) to halve the time. #Define the number of global scenarios to be investigated (e.g. prey 1 pulses, prey 2 pulses) scen.names = expand.grid(P.scen=c("Pulse1","Pulse2","Pulse3"),I.scen=c("mig"),F.scen=c("decrF")) # "highF","lowF", scen.names$scenario = c(1:nrow(scen.names)) scen.names$scen.name = paste(scen.names$P.scen, scen.names$I.scen, scen.names$F.scen, sep="-") scen.names = subset(scen.names, select=c("scenario","scen.name","P.scen","I.scen","F.scen")) #Re-order columns to match the previous version scenarios = nrow(scen.names) ## CREATE EMPTY ARRAYS TO FILL WITH SENSITIVITY RESULTS #Summary of %diff in SSB (S), Yield (Y), and weight at age (w) S.sens = array(NA, dim=c(sims,scenarios,sens.scenarios),dimnames=list(paste("sim",1:sims),scen.names[,2],sens.names$name)) Y.sens = array(NA, dim=c(sims,scenarios,sens.scenarios),dimnames=list(paste("sim",1:sims),scen.names[,2],sens.names$name)) R.sens = array(NA, dim=c(sims,scenarios,sens.scenarios),dimnames=list(paste("sim",1:sims),scen.names[,2],sens.names$name)) w.sens = array(NA, dim=c(sims,ages,scenarios,sens.scenarios),dimnames=list(paste("sim",1:sims),paste("age-",0:(ages-1),sep=""),scen.names[,2],sens.names$name)) t.sens = array(NA, dim=c(sims,scenarios,2,sens.scenarios),dimnames=list(paste("sim",1:sims),scen.names[,2], c("NoPulse","Pulse"),sens.names$name)) dim(S.sens); length(S.sens) dim(w.sens); length(w.sens) ## RUN THROUGH SENSITIVITY SCENARIOS (tweaking each parameter up or down based on a given percentage, e.g. 20% 9-26-13 1:25:30 #10-8-13 - added Recruitment as an output #BASE model temp.base = simfxn(graph=F, sims=sims) S.sens[,,"BASE."] = temp.base$S.diff.p.mean Y.sens[,,"BASE."] = temp.base$Y.diff.p.mean R.sens[,,"BASE."] = temp.base$R.diff.p.mean w.sens[,,,"BASE."] = temp.base$w.diff.p.mean t.sens[,,,"BASE."] = aperm(temp.base$time.target,c(2,1,3)) #Mortality sensitivities #Natural mortality temp.lo = simfxn(M = BASE.def[BASE.def$Parm=="M","Val.lo"], graph=F, sims=sims) S.sens[,,"M.lo"] = temp.lo$S.diff.p.mean Y.sens[,,"M.lo"] = temp.lo$Y.diff.p.mean R.sens[,,"M.lo"] = temp.lo$R.diff.p.mean w.sens[,,,"M.lo"] = temp.lo$w.diff.p.mean t.sens[,,,"M.lo"] = aperm(temp.lo$time.target,c(2,1,3)) temp.hi = simfxn(M = BASE.def[BASE.def$Parm=="M","Val.hi"], graph=F, sims=sims) S.sens[,,"M.hi"] = temp.hi$S.diff.p.mean Y.sens[,,"M.hi"] = temp.hi$Y.diff.p.mean R.sens[,,"M.hi"] = temp.hi$R.diff.p.mean w.sens[,,,"M.hi"] = temp.hi$w.diff.p.mean t.sens[,,,"M.hi"] = aperm(temp.hi$time.target,c(2,1,3)) #Stock Recruitment sensitivities #Alpha, K, beta parameters temp.lo = simfxn(alpha.SR = BASE.def[BASE.def$Parm=="alpha.SR","Val.lo"], graph=F, sims=sims) S.sens[,,"alpha.SR.lo"] = temp.lo$S.diff.p.mean Y.sens[,,"alpha.SR.lo"] = temp.lo$Y.diff.p.mean R.sens[,,"alpha.SR.lo"] = temp.lo$R.diff.p.mean w.sens[,,,"alpha.SR.lo"] = temp.lo$w.diff.p.mean t.sens[,,,"alpha.SR.lo"] = aperm(temp.lo$time.target,c(2,1,3)) temp.hi = simfxn(alpha.SR = BASE.def[BASE.def$Parm=="alpha.SR","Val.hi"], graph=F, sims=sims) S.sens[,,"alpha.SR.hi"] = temp.hi$S.diff.p.mean Y.sens[,,"alpha.SR.hi"] = temp.hi$Y.diff.p.mean R.sens[,,"alpha.SR.hi"] = temp.hi$R.diff.p.mean w.sens[,,,"alpha.SR.hi"] = temp.hi$w.diff.p.mean t.sens[,,,"alpha.SR.hi"] = aperm(temp.hi$time.target,c(2,1,3)) temp.lo = simfxn(K.SR = BASE.def[BASE.def$Parm=="K.SR","Val.lo"], graph=F, sims=sims) S.sens[,,"K.SR.lo"] = temp.lo$S.diff.p.mean Y.sens[,,"K.SR.lo"] = temp.lo$Y.diff.p.mean R.sens[,,"K.SR.lo"] = temp.lo$R.diff.p.mean w.sens[,,,"K.SR.lo"] = temp.lo$w.diff.p.mean t.sens[,,,"K.SR.lo"] = aperm(temp.lo$time.target,c(2,1,3)) temp.hi = simfxn(K.SR = BASE.def[BASE.def$Parm=="K.SR","Val.hi"], graph=F, sims=sims) S.sens[,,"K.SR.hi"] = temp.hi$S.diff.p.mean Y.sens[,,"K.SR.hi"] = temp.hi$Y.diff.p.mean R.sens[,,"K.SR.hi"] = temp.hi$R.diff.p.mean w.sens[,,,"K.SR.hi"] = temp.hi$w.diff.p.mean t.sens[,,,"K.SR.hi"] = aperm(temp.hi$time.target,c(2,1,3)) temp.lo = simfxn(beta.SR = BASE.def[BASE.def$Parm=="beta.SR","Val.lo"], graph=F, sims=sims) S.sens[,,"beta.SR.lo"] = temp.lo$S.diff.p.mean Y.sens[,,"beta.SR.lo"] = temp.lo$Y.diff.p.mean R.sens[,,"beta.SR.lo"] = temp.lo$R.diff.p.mean w.sens[,,,"beta.SR.lo"] = temp.lo$w.diff.p.mean t.sens[,,,"beta.SR.lo"] = aperm(temp.lo$time.target,c(2,1,3)) temp.hi = simfxn(beta.SR = BASE.def[BASE.def$Parm=="beta.SR","Val.hi"], graph=F, sims=sims) S.sens[,, "beta.SR.hi"] = temp.hi$S.diff.p.mean Y.sens[,, "beta.SR.hi"] = temp.hi$Y.diff.p.mean R.sens[,, "beta.SR.hi"] = temp.hi$R.diff.p.mean w.sens[,,,"beta.SR.hi"] = temp.hi$w.diff.p.mean t.sens[,,,"beta.SR.hi"] = aperm(temp.hi$time.target,c(2,1,3)) #Deviate temp.lo = simfxn(R.dev.sd = BASE.def[BASE.def$Parm=="R.dev.sd","Val.lo"], graph=F, sims=sims) S.sens[,, "R.dev.sd.lo"] = temp.lo$S.diff.p.mean Y.sens[,, "R.dev.sd.lo"] = temp.lo$Y.diff.p.mean R.sens[,, "R.dev.sd.lo"] = temp.lo$R.diff.p.mean w.sens[,,,"R.dev.sd.lo"] = temp.lo$w.diff.p.mean t.sens[,,,"R.dev.sd.lo"] = aperm(temp.lo$time.target,c(2,1,3)) temp.hi = simfxn(R.dev.sd = BASE.def[BASE.def$Parm=="R.dev.sd","Val.hi"], graph=F, sims=sims) S.sens[,, "R.dev.sd.hi"] = temp.hi$S.diff.p.mean Y.sens[,, "R.dev.sd.hi"] = temp.hi$Y.diff.p.mean R.sens[,, "R.dev.sd.hi"] = temp.hi$R.diff.p.mean w.sens[,,,"R.dev.sd.hi"] = temp.hi$w.diff.p.mean t.sens[,,,"R.dev.sd.hi"] = aperm(temp.hi$time.target,c(2,1,3)) #Prey sensitivities #Pulse Parameters temp.lo = simfxn(n.pulse = BASE.def[BASE.def$Parm=="n.pulse","Val.hi"], graph=F, sims=sims) S.sens[,, "n.pulse.lo"] = temp.lo$S.diff.p.mean Y.sens[,, "n.pulse.lo"] = temp.lo$Y.diff.p.mean R.sens[,, "n.pulse.lo"] = temp.lo$R.diff.p.mean w.sens[,,,"n.pulse.lo"] = temp.lo$w.diff.p.mean t.sens[,,,"n.pulse.lo"] = aperm(temp.lo$time.target,c(2,1,3)) temp.hi = simfxn(n.pulse = BASE.def[BASE.def$Parm=="n.pulse","Val.lo"], graph=F, sims=sims) S.sens[,, "n.pulse.hi"] = temp.hi$S.diff.p.mean Y.sens[,, "n.pulse.hi"] = temp.hi$Y.diff.p.mean R.sens[,, "n.pulse.hi"] = temp.hi$R.diff.p.mean w.sens[,,,"n.pulse.hi"] = temp.hi$w.diff.p.mean t.sens[,,,"n.pulse.hi"] = aperm(temp.hi$time.target,c(2,1,3)) temp.lo = simfxn(pulse.range = c(2:4), graph=F, sims=sims) #NOTE: This was not 20% change, but instead the lower half of the range S.sens[,, "pulse.range.lo"] = temp.lo$S.diff.p.mean Y.sens[,, "pulse.range.lo"] = temp.lo$Y.diff.p.mean R.sens[,, "pulse.range.lo"] = temp.lo$R.diff.p.mean w.sens[,,,"pulse.range.lo"] = temp.lo$w.diff.p.mean t.sens[,,,"pulse.range.lo"] = aperm(temp.lo$time.target,c(2,1,3)) temp.hi = simfxn(pulse.range = c(4:6), graph=F, sims=sims) #NOTE: This was not 20% change, but instead the upper half of the range S.sens[,, "pulse.range.hi"] = temp.hi$S.diff.p.mean Y.sens[,, "pulse.range.hi"] = temp.hi$Y.diff.p.mean R.sens[,, "pulse.range.hi"] = temp.hi$R.diff.p.mean w.sens[,,,"pulse.range.hi"] = temp.hi$w.diff.p.mean t.sens[,,,"pulse.range.hi"] = aperm(temp.hi$time.target,c(2,1,3)) #Prey biomass values temp.lo = simfxn(N.prey.2.mean = BASE.def[BASE.def$Parm=="N.prey.2.mean","Val.lo"], graph=F, sims=sims) S.sens[,, "N.prey.2.mean.lo"] = temp.lo$S.diff.p.mean Y.sens[,, "N.prey.2.mean.lo"] = temp.lo$Y.diff.p.mean R.sens[,, "N.prey.2.mean.lo"] = temp.lo$R.diff.p.mean w.sens[,,,"N.prey.2.mean.lo"] = temp.lo$w.diff.p.mean t.sens[,,,"N.prey.2.mean.lo"] = aperm(temp.lo$time.target,c(2,1,3)) temp.hi = simfxn(N.prey.2.mean = BASE.def[BASE.def$Parm=="N.prey.2.mean","Val.hi"], graph=F, sims=sims) S.sens[,, "N.prey.2.mean.hi"] = temp.hi$S.diff.p.mean Y.sens[,, "N.prey.2.mean.hi"] = temp.hi$Y.diff.p.mean R.sens[,, "N.prey.2.mean.hi"] = temp.hi$R.diff.p.mean w.sens[,,,"N.prey.2.mean.hi"] = temp.hi$w.diff.p.mean t.sens[,,,"N.prey.2.mean.hi"] = aperm(temp.hi$time.target,c(2,1,3)) temp.lo = simfxn(relB.1 = BASE.def[BASE.def$Parm=="relB.1","Val.lo"], graph=F, sims=sims) S.sens[,, "relB.1.lo"] = temp.lo$S.diff.p.mean Y.sens[,, "relB.1.lo"] = temp.lo$Y.diff.p.mean R.sens[,, "relB.1.lo"] = temp.lo$R.diff.p.mean w.sens[,,,"relB.1.lo"] = temp.lo$w.diff.p.mean t.sens[,,,"relB.1.lo"] = aperm(temp.lo$time.target,c(2,1,3)) temp.hi = simfxn(relB.1 = BASE.def[BASE.def$Parm=="relB.1","Val.hi"], graph=F, sims=sims) S.sens[,, "relB.1.hi"] = temp.hi$S.diff.p.mean Y.sens[,, "relB.1.hi"] = temp.hi$Y.diff.p.mean R.sens[,, "relB.1.hi"] = temp.hi$R.diff.p.mean w.sens[,,,"relB.1.hi"] = temp.hi$w.diff.p.mean t.sens[,,,"relB.1.hi"] = aperm(temp.hi$time.target,c(2,1,3)) temp.lo = simfxn(relB.3 = BASE.def[BASE.def$Parm=="relB.3","Val.lo"], graph=F, sims=sims) S.sens[,, "relB.3.lo"] = temp.lo$S.diff.p.mean Y.sens[,, "relB.3.lo"] = temp.lo$Y.diff.p.mean R.sens[,, "relB.3.lo"] = temp.lo$R.diff.p.mean w.sens[,,,"relB.3.lo"] = temp.lo$w.diff.p.mean t.sens[,,,"relB.3.lo"] = aperm(temp.lo$time.target,c(2,1,3)) temp.hi = simfxn(relB.3 = BASE.def[BASE.def$Parm=="relB.3","Val.hi"], graph=F, sims=sims) S.sens[,, "relB.3.hi"] = temp.hi$S.diff.p.mean Y.sens[,, "relB.3.hi"] = temp.hi$Y.diff.p.mean R.sens[,, "relB.3.hi"] = temp.hi$R.diff.p.mean w.sens[,,,"relB.3.hi"] = temp.hi$w.diff.p.mean t.sens[,,,"relB.3.hi"] = aperm(temp.hi$time.target,c(2,1,3)) #Stochastic prey deviations temp.lo = simfxn(N.prey.dev = BASE.def[BASE.def$Parm=="N.prey.dev","Val.lo"], graph=F, sims=sims) S.sens[,, "N.prey.dev.lo"] = temp.lo$S.diff.p.mean Y.sens[,, "N.prey.dev.lo"] = temp.lo$Y.diff.p.mean R.sens[,, "N.prey.dev.lo"] = temp.lo$R.diff.p.mean w.sens[,,,"N.prey.dev.lo"] = temp.lo$w.diff.p.mean t.sens[,,,"N.prey.dev.lo"] = aperm(temp.lo$time.target,c(2,1,3)) temp.hi = simfxn(N.prey.dev = BASE.def[BASE.def$Parm=="N.prey.dev","Val.hi"], graph=F, sims=sims) S.sens[,, "N.prey.dev.hi"] = temp.hi$S.diff.p.mean Y.sens[,, "N.prey.dev.hi"] = temp.hi$Y.diff.p.mean R.sens[,, "N.prey.dev.hi"] = temp.hi$R.diff.p.mean w.sens[,,,"N.prey.dev.hi"] = temp.hi$w.diff.p.mean t.sens[,,,"N.prey.dev.hi"] = aperm(temp.hi$time.target,c(2,1,3)) #Growth submodel sensitivities #Cmax parameters temp.lo = simfxn(CA = BASE.def[BASE.def$Parm=="CA","Val.lo"], graph=F, sims=sims) S.sens[,, "CA.lo"] = temp.lo$S.diff.p.mean Y.sens[,, "CA.lo"] = temp.lo$Y.diff.p.mean R.sens[,, "CA.lo"] = temp.lo$R.diff.p.mean w.sens[,,,"CA.lo"] = temp.lo$w.diff.p.mean t.sens[,,,"CA.lo"] = aperm(temp.lo$time.target,c(2,1,3)) temp.hi = simfxn(CA = BASE.def[BASE.def$Parm=="CA","Val.hi"], graph=F, sims=sims) S.sens[,, "CA.hi"] = temp.hi$S.diff.p.mean Y.sens[,, "CA.hi"] = temp.hi$Y.diff.p.mean R.sens[,, "CA.hi"] = temp.hi$R.diff.p.mean w.sens[,,,"CA.hi"] = temp.hi$w.diff.p.mean t.sens[,,,"CA.hi"] = aperm(temp.hi$time.target,c(2,1,3)) temp.lo = simfxn(CB = BASE.def[BASE.def$Parm=="CB","Val.lo"], graph=F, sims=sims) S.sens[,, "CB.lo"] = temp.lo$S.diff.p.mean Y.sens[,, "CB.lo"] = temp.lo$Y.diff.p.mean R.sens[,, "CB.lo"] = temp.lo$R.diff.p.mean w.sens[,,,"CB.lo"] = temp.lo$w.diff.p.mean t.sens[,,,"CB.lo"] = aperm(temp.lo$time.target,c(2,1,3)) temp.hi = simfxn(CB = BASE.def[BASE.def$Parm=="CB","Val.hi"], graph=F, sims=sims) S.sens[,, "CB.hi"] = temp.hi$S.diff.p.mean Y.sens[,, "CB.hi"] = temp.hi$Y.diff.p.mean R.sens[,, "CB.hi"] = temp.hi$R.diff.p.mean w.sens[,,,"CB.hi"] = temp.hi$w.diff.p.mean t.sens[,,,"CB.hi"] = aperm(temp.hi$time.target,c(2,1,3)) #Temperature parameters temp.lo = simfxn(CTO = BASE.def[BASE.def$Parm=="CTO","Val.lo"], graph=F, sims=sims) S.sens[,, "CTO.lo"] = temp.lo$S.diff.p.mean Y.sens[,, "CTO.lo"] = temp.lo$Y.diff.p.mean R.sens[,, "CTO.lo"] = temp.lo$R.diff.p.mean w.sens[,,,"CTO.lo"] = temp.lo$w.diff.p.mean t.sens[,,,"CTO.lo"] = aperm(temp.lo$time.target,c(2,1,3)) temp.hi = simfxn(CTO = BASE.def[BASE.def$Parm=="CTO","Val.hi"], graph=F, sims=sims) S.sens[,, "CTO.hi"] = temp.hi$S.diff.p.mean Y.sens[,, "CTO.hi"] = temp.hi$Y.diff.p.mean R.sens[,, "CTO.hi"] = temp.hi$R.diff.p.mean w.sens[,,,"CTO.hi"] = temp.hi$w.diff.p.mean t.sens[,,,"CTO.hi"] = aperm(temp.hi$time.target,c(2,1,3)) temp.lo = simfxn(CTM = BASE.def[BASE.def$Parm=="CTM","Val.lo"], graph=F, sims=sims) S.sens[,, "CTM.lo"] = temp.lo$S.diff.p.mean Y.sens[,, "CTM.lo"] = temp.lo$Y.diff.p.mean R.sens[,, "CTM.lo"] = temp.lo$R.diff.p.mean w.sens[,,,"CTM.lo"] = temp.lo$w.diff.p.mean t.sens[,,,"CTM.lo"] = aperm(temp.lo$time.target,c(2,1,3)) temp.hi = simfxn(CTM = BASE.def[BASE.def$Parm=="CTM","Val.hi"], graph=F, sims=sims) S.sens[,, "CTM.hi"] = temp.hi$S.diff.p.mean Y.sens[,, "CTM.hi"] = temp.hi$Y.diff.p.mean R.sens[,, "CTM.hi"] = temp.hi$R.diff.p.mean w.sens[,,,"CTM.hi"] = temp.hi$w.diff.p.mean t.sens[,,,"CTM.hi"] = aperm(temp.hi$time.target,c(2,1,3)) #1/2 Way (half of all parameter perturbations 9-26-13 1:35am -?? BASE.model = temp.base #Resave the BASE run with more informative name #Save base model run, and summary tables SSB, Yield, and Weight at age sensitivities save(BASE.model, S.sens, Y.sens, R.sens, w.sens, t.sens, file="C:\\Users\\andrebuc\\VIMS\\PhD\\Research\\PredPrey_Simulations\\SimOutputs\\SensitivityOutput_Method2_part1_10-8-13.RData") temp.lo = simfxn(CQ = BASE.def[BASE.def$Parm=="CQ","Val.lo"], graph=F, sims=sims) S.sens[,, "CQ.lo"] = temp.lo$S.diff.p.mean Y.sens[,, "CQ.lo"] = temp.lo$Y.diff.p.mean R.sens[,, "CQ.lo"] = temp.lo$R.diff.p.mean w.sens[,,,"CQ.lo"] = temp.lo$w.diff.p.mean t.sens[,,,"CQ.lo"] = aperm(temp.lo$time.target,c(2,1,3)) temp.hi = simfxn(CQ = BASE.def[BASE.def$Parm=="CQ","Val.hi"], graph=F, sims=sims) S.sens[,, "CQ.hi"] = temp.hi$S.diff.p.mean Y.sens[,, "CQ.hi"] = temp.hi$Y.diff.p.mean R.sens[,, "CQ.hi"] = temp.hi$R.diff.p.mean w.sens[,,,"CQ.hi"] = temp.hi$w.diff.p.mean t.sens[,,,"CQ.hi"] = aperm(temp.hi$time.target,c(2,1,3)) #Conversion efficiency (K) parameters temp.lo = simfxn(KL = BASE.def[BASE.def$Parm=="KL","Val.lo"], graph=F, sims=sims) S.sens[,, "KL.lo"] = temp.lo$S.diff.p.mean Y.sens[,, "KL.lo"] = temp.lo$Y.diff.p.mean R.sens[,, "KL.lo"] = temp.lo$R.diff.p.mean w.sens[,,,"KL.lo"] = temp.lo$w.diff.p.mean t.sens[,,,"KL.lo"] = aperm(temp.lo$time.target,c(2,1,3)) temp.hi = simfxn(KL = BASE.def[BASE.def$Parm=="KL","Val.hi"], graph=F, sims=sims) S.sens[,, "KL.hi"] = temp.hi$S.diff.p.mean Y.sens[,, "KL.hi"] = temp.hi$Y.diff.p.mean R.sens[,, "KL.hi"] = temp.hi$R.diff.p.mean w.sens[,,,"KL.hi"] = temp.hi$w.diff.p.mean t.sens[,,,"KL.hi"] = aperm(temp.hi$time.target,c(2,1,3)) temp.lo = simfxn(KR = BASE.def[BASE.def$Parm=="KR","Val.lo"], graph=F, sims=sims) # S.sens[,, "KR.lo"] = temp.lo$S.diff.p.mean Y.sens[,, "KR.lo"] = temp.lo$Y.diff.p.mean R.sens[,, "KR.lo"] = temp.lo$R.diff.p.mean w.sens[,,,"KR.lo"] = temp.lo$w.diff.p.mean t.sens[,,,"KR.lo"] = aperm(temp.lo$time.target,c(2,1,3)) temp.hi = simfxn(KR = BASE.def[BASE.def$Parm=="KR","Val.hi"], graph=F, sims=sims) S.sens[,, "KR.hi"] = temp.hi$S.diff.p.mean Y.sens[,, "KR.hi"] = temp.hi$Y.diff.p.mean R.sens[,, "KR.hi"] = temp.hi$R.diff.p.mean w.sens[,,,"KR.hi"] = temp.hi$w.diff.p.mean t.sens[,,,"KR.hi"] = aperm(temp.hi$time.target,c(2,1,3)) temp.lo = simfxn(KW = BASE.def[BASE.def$Parm=="KW","Val.lo"], graph=F, sims=sims) S.sens[,, "KW.lo"] = temp.lo$S.diff.p.mean Y.sens[,, "KW.lo"] = temp.lo$Y.diff.p.mean R.sens[,, "KW.lo"] = temp.lo$R.diff.p.mean w.sens[,,,"KW.lo"] = temp.lo$w.diff.p.mean t.sens[,,,"KW.lo"] = aperm(temp.lo$time.target,c(2,1,3)) temp.hi = simfxn(KW = BASE.def[BASE.def$Parm=="KW","Val.hi"], graph=F, sims=sims) S.sens[,, "KW.hi"] = temp.hi$S.diff.p.mean Y.sens[,, "KW.hi"] = temp.hi$Y.diff.p.mean R.sens[,, "KW.hi"] = temp.hi$R.diff.p.mean w.sens[,,,"KW.hi"] = temp.hi$w.diff.p.mean t.sens[,,,"KW.hi"] = aperm(temp.hi$time.target,c(2,1,3)) #Functional response sensitivities #Max attack rates temp.lo = simfxn(atck.max.1 = BASE.def[BASE.def$Parm=="atck.max.1","Val.lo"], graph=F, sims=sims) S.sens[,, "atck.max.1.lo"] = temp.lo$S.diff.p.mean Y.sens[,, "atck.max.1.lo"] = temp.lo$Y.diff.p.mean R.sens[,, "atck.max.1.lo"] = temp.lo$R.diff.p.mean w.sens[,,,"atck.max.1.lo"] = temp.lo$w.diff.p.mean t.sens[,,,"atck.max.1.lo"] = aperm(temp.lo$time.target,c(2,1,3)) temp.hi = simfxn(atck.max.1 = BASE.def[BASE.def$Parm=="atck.max.1","Val.hi"], graph=F, sims=sims) S.sens[,, "atck.max.1.hi"] = temp.hi$S.diff.p.mean Y.sens[,, "atck.max.1.hi"] = temp.hi$Y.diff.p.mean R.sens[,, "atck.max.1.hi"] = temp.hi$R.diff.p.mean w.sens[,,,"atck.max.1.hi"] = temp.hi$w.diff.p.mean t.sens[,,,"atck.max.1.hi"] = aperm(temp.hi$time.target,c(2,1,3)) temp.lo = simfxn(atck.max.2a = BASE.def[BASE.def$Parm=="atck.max.2a","Val.lo"], graph=F, sims=sims) S.sens[,,"atck.max.2a.lo"] = temp.lo$S.diff.p.mean Y.sens[,,"atck.max.2a.lo"] = temp.lo$Y.diff.p.mean R.sens[,,"atck.max.2a.lo"] = temp.lo$R.diff.p.mean w.sens[,,,"atck.max.2a.lo"] = temp.lo$w.diff.p.mean t.sens[,,,"atck.max.2a.lo"] = aperm(temp.lo$time.target,c(2,1,3)) temp.hi = simfxn(atck.max.2a = BASE.def[BASE.def$Parm=="atck.max.2a","Val.hi"], graph=F, sims=sims) S.sens[,,"atck.max.2a.hi"] = temp.hi$S.diff.p.mean Y.sens[,,"atck.max.2a.hi"] = temp.hi$Y.diff.p.mean R.sens[,,"atck.max.2a.hi"] = temp.hi$R.diff.p.mean w.sens[,,,"atck.max.2a.hi"] = temp.hi$w.diff.p.mean t.sens[,,,"atck.max.2a.hi"] = aperm(temp.hi$time.target,c(2,1,3)) temp.lo = simfxn(atck.max.3 = BASE.def[BASE.def$Parm=="atck.max.3","Val.lo"], graph=F, sims=sims) S.sens[,,"atck.max.3.lo"] = temp.lo$S.diff.p.mean Y.sens[,,"atck.max.3.lo"] = temp.lo$Y.diff.p.mean R.sens[,,"atck.max.3.lo"] = temp.lo$R.diff.p.mean w.sens[,,,"atck.max.3.lo"] = temp.lo$w.diff.p.mean t.sens[,,,"atck.max.3.lo"] = aperm(temp.lo$time.target,c(2,1,3)) temp.hi = simfxn(atck.max.3 = BASE.def[BASE.def$Parm=="atck.max.3","Val.hi"], graph=F, sims=sims) S.sens[,,"atck.max.3.hi"] = temp.hi$S.diff.p.mean Y.sens[,,"atck.max.3.hi"] = temp.hi$Y.diff.p.mean R.sens[,,"atck.max.3.hi"] = temp.hi$R.diff.p.mean w.sens[,,,"atck.max.3.hi"] = temp.hi$w.diff.p.mean t.sens[,,,"atck.max.3.hi"] = aperm(temp.hi$time.target,c(2,1,3)) #Rho of attack rates temp.lo = simfxn(atck.rho.1 = BASE.def[BASE.def$Parm=="atck.rho.1","Val.lo"], graph=F, sims=sims) S.sens[,,"atck.rho.1.lo"] = temp.lo$S.diff.p.mean Y.sens[,,"atck.rho.1.lo"] = temp.lo$Y.diff.p.mean R.sens[,,"atck.rho.1.lo"] = temp.lo$R.diff.p.mean w.sens[,,,"atck.rho.1.lo"] = temp.lo$w.diff.p.mean t.sens[,,,"atck.rho.1.lo"] = aperm(temp.lo$time.target,c(2,1,3)) temp.hi = simfxn(atck.rho.1 = BASE.def[BASE.def$Parm=="atck.rho.1","Val.hi"], graph=F, sims=sims) S.sens[,,"atck.rho.1.hi"] = temp.hi$S.diff.p.mean Y.sens[,,"atck.rho.1.hi"] = temp.hi$Y.diff.p.mean R.sens[,,"atck.rho.1.hi"] = temp.hi$R.diff.p.mean w.sens[,,,"atck.rho.1.hi"] = temp.hi$w.diff.p.mean t.sens[,,,"atck.rho.1.hi"] = aperm(temp.hi$time.target,c(2,1,3)) temp.lo = simfxn(atck.rho.2a = BASE.def[BASE.def$Parm=="atck.rho.2a","Val.lo"], atck.rho.2b = BASE.def[BASE.def$Parm=="atck.rho.2a","Val.lo"], graph=F, sims=sims) S.sens[,,"atck.rho.2a.lo"] = temp.lo$S.diff.p.mean #For this case, i'm setting atck.rho.2b = atck.rho.2a, since they are assumed to be equal. Y.sens[,,"atck.rho.2a.lo"] = temp.lo$Y.diff.p.mean R.sens[,,"atck.rho.2a.lo"] = temp.lo$R.diff.p.mean w.sens[,,,"atck.rho.2a.lo"] = temp.lo$w.diff.p.mean t.sens[,,,"atck.rho.2a.lo"] = aperm(temp.lo$time.target,c(2,1,3)) temp.hi = simfxn(atck.rho.2a = BASE.def[BASE.def$Parm=="atck.rho.2a","Val.hi"], atck.rho.2b = BASE.def[BASE.def$Parm=="atck.rho.2a","Val.hi"], graph=F, sims=sims) S.sens[,,"atck.rho.2a.hi"] = temp.hi$S.diff.p.mean #For this case, i'm setting atck.rho.2b = atck.rho.2a, since they are assumed to be equal. Y.sens[,,"atck.rho.2a.hi"] = temp.hi$Y.diff.p.mean R.sens[,,"atck.rho.2a.hi"] = temp.hi$R.diff.p.mean w.sens[,,,"atck.rho.2a.hi"] = temp.hi$w.diff.p.mean t.sens[,,,"atck.rho.2a.hi"] = aperm(temp.hi$time.target,c(2,1,3)) temp.lo = simfxn(atck.rho.3 = BASE.def[BASE.def$Parm=="atck.rho.3","Val.lo"], graph=F, sims=sims) S.sens[,,"atck.rho.3.lo"] = temp.lo$S.diff.p.mean Y.sens[,,"atck.rho.3.lo"] = temp.lo$Y.diff.p.mean R.sens[,,"atck.rho.3.lo"] = temp.lo$R.diff.p.mean w.sens[,,,"atck.rho.3.lo"] = temp.lo$w.diff.p.mean t.sens[,,,"atck.rho.3.lo"] = aperm(temp.lo$time.target,c(2,1,3)) temp.hi = simfxn(atck.rho.3 = BASE.def[BASE.def$Parm=="atck.rho.3","Val.hi"], graph=F, sims=sims) S.sens[,,"atck.rho.3.hi"] = temp.hi$S.diff.p.mean Y.sens[,,"atck.rho.3.hi"] = temp.hi$Y.diff.p.mean R.sens[,,"atck.rho.3.hi"] = temp.hi$R.diff.p.mean w.sens[,,,"atck.rho.3.hi"] = temp.hi$w.diff.p.mean t.sens[,,,"atck.rho.3.hi"] = aperm(temp.hi$time.target,c(2,1,3)) #Inflection point of attack rates temp.lo = simfxn(atck.infl.1 = BASE.def[BASE.def$Parm=="atck.infl.1","Val.lo"], graph=F, sims=sims) S.sens[,,"atck.infl.1.lo"] = temp.lo$S.diff.p.mean Y.sens[,,"atck.infl.1.lo"] = temp.lo$Y.diff.p.mean R.sens[,,"atck.infl.1.lo"] = temp.lo$R.diff.p.mean w.sens[,,,"atck.infl.1.lo"] = temp.lo$w.diff.p.mean t.sens[,,,"atck.infl.1.lo"] = aperm(temp.lo$time.target,c(2,1,3)) temp.hi = simfxn(atck.infl.1 = BASE.def[BASE.def$Parm=="atck.infl.1","Val.hi"], graph=F, sims=sims) S.sens[,,"atck.infl.1.hi"] = temp.hi$S.diff.p.mean Y.sens[,,"atck.infl.1.hi"] = temp.hi$Y.diff.p.mean R.sens[,,"atck.infl.1.hi"] = temp.hi$R.diff.p.mean w.sens[,,,"atck.infl.1.hi"] = temp.hi$w.diff.p.mean t.sens[,,,"atck.infl.1.hi"] = aperm(temp.hi$time.target,c(2,1,3)) temp.lo = simfxn(atck.infl.2a = BASE.def[BASE.def$Parm=="atck.infl.2a","Val.lo"], graph=F, sims=sims) S.sens[,,"atck.infl.2a.lo"] = temp.lo$S.diff.p.mean Y.sens[,,"atck.infl.2a.lo"] = temp.lo$Y.diff.p.mean R.sens[,,"atck.infl.2a.lo"] = temp.lo$R.diff.p.mean w.sens[,,,"atck.infl.2a.lo"] = temp.lo$w.diff.p.mean t.sens[,,,"atck.infl.2a.lo"] = aperm(temp.lo$time.target,c(2,1,3)) temp.hi = simfxn(atck.infl.2a = BASE.def[BASE.def$Parm=="atck.infl.2a","Val.hi"], graph=F, sims=sims) S.sens[,,"atck.infl.2a.hi"] = temp.hi$S.diff.p.mean Y.sens[,,"atck.infl.2a.hi"] = temp.hi$Y.diff.p.mean R.sens[,,"atck.infl.2a.hi"] = temp.hi$R.diff.p.mean w.sens[,,,"atck.infl.2a.hi"] = temp.hi$w.diff.p.mean t.sens[,,,"atck.infl.2a.hi"] = aperm(temp.hi$time.target,c(2,1,3)) temp.lo = simfxn(atck.infl.2b = BASE.def[BASE.def$Parm=="atck.infl.2b","Val.lo"], graph=F, sims=sims) S.sens[,,"atck.infl.2b.lo"] = temp.lo$S.diff.p.mean Y.sens[,,"atck.infl.2b.lo"] = temp.lo$Y.diff.p.mean R.sens[,,"atck.infl.2b.lo"] = temp.lo$R.diff.p.mean w.sens[,,,"atck.infl.2b.lo"] = temp.lo$w.diff.p.mean t.sens[,,,"atck.infl.2b.lo"] = aperm(temp.lo$time.target,c(2,1,3)) temp.hi = simfxn(atck.infl.2b = BASE.def[BASE.def$Parm=="atck.infl.2b","Val.hi"], graph=F, sims=sims) S.sens[,,"atck.infl.2b.hi"] = temp.hi$S.diff.p.mean Y.sens[,,"atck.infl.2b.hi"] = temp.hi$Y.diff.p.mean R.sens[,,"atck.infl.2b.hi"] = temp.hi$R.diff.p.mean w.sens[,,,"atck.infl.2b.hi"] = temp.hi$w.diff.p.mean t.sens[,,,"atck.infl.2b.hi"] = aperm(temp.hi$time.target,c(2,1,3)) temp.lo = simfxn(atck.infl.3 = BASE.def[BASE.def$Parm=="atck.infl.3","Val.lo"], graph=F, sims=sims) S.sens[,,"atck.infl.3.lo"] = temp.lo$S.diff.p.mean Y.sens[,,"atck.infl.3.lo"] = temp.lo$Y.diff.p.mean R.sens[,,"atck.infl.3.lo"] = temp.lo$R.diff.p.mean w.sens[,,,"atck.infl.3.lo"] = temp.lo$w.diff.p.mean t.sens[,,,"atck.infl.3.lo"] = aperm(temp.lo$time.target,c(2,1,3)) temp.hi = simfxn(atck.infl.3 = BASE.def[BASE.def$Parm=="atck.infl.3","Val.hi"], graph=F, sims=sims) S.sens[,,"atck.infl.3.hi"] = temp.hi$S.diff.p.mean Y.sens[,,"atck.infl.3.hi"] = temp.hi$Y.diff.p.mean R.sens[,,"atck.infl.3.hi"] = temp.hi$R.diff.p.mean w.sens[,,,"atck.infl.3.hi"] = temp.hi$w.diff.p.mean t.sens[,,,"atck.infl.3.hi"] = aperm(temp.hi$time.target,c(2,1,3)) #Second half of all parameter perturbations #Save base model run, and summary tables SSB, Yield, and Weight at age sensitivities save(S.sens, Y.sens, R.sens, w.sens, t.sens, file="C:\\Users\\andrebuc\\VIMS\\PhD\\Research\\PredPrey_Simulations\\SimOutputs\\SensitivityOutput_Method2_part2_10-7-13.RData") #JUNK temp.lo$S.diff.mean temp.lo$S.diff.p.mean temp.base$S.diff.mean temp.base$S.diff.p.mean S.sens[,,c("alpha.SR.hi","BASE.")] #Store data objects for future use and reference #summary(S.sens); summary(Y.sens); summary(w.sens) #Check that there are no NA's #BASE.model = temp.base #Resave the BASE run with more informative name #Save base model run, and summary tables SSB, Yield, and Weight at age sensitivities #save(BASE.model, S.sens, Y.sens, w.sens, # file="C:\\Users\\andrebuc\\VIMS\\PhD\\Research\\PredPrey_Simulations\\SimOutputs\\SensitivityOutput_2-6-13.RData") #backup.S = S.sens; backup.Y = Y.sens; backup.w = w.sens ################################################ ## IMPORT SENSITIVITY RESULTS FOR PLOTTING ################################################ ## METHOD 1 - MONTE CARLO SENSITIVITY #Monte Carlo Method - S.sens, Y.sens, w.sens, t.sens load("C:\\Users\\andre\\VIMS\\PhD\\Research\\PredPrey_Simulations\\SimOutputs\\SensitivityOutput_MC_10-5-13.RData") summary(S.sens); summary(Y.sens); summary(t.sens); summary(w.sens); dim(S.sens) ## METHOD 2 - INDIVIDUAL PARAMETER PERTURBATION # - Load the two halves of the Method2 sensitivity results and merge them (I ran all the models in two sections (tinnR and Rstudio) to speed it up) #Part 1 - BASE.model, S.sens, Y.sens, w.sens, t.sens load("C:\\Users\\andre\\VIMS\\PhD\\Research\\PredPrey_Simulations\\SimOutputs\\SensitivityOutput_Method2_part1_8-4-14.RData") S.sens.1 = S.sens; Y.sens.1 = Y.sens; R.sens.1 = R.sens; w.sens.1 = w.sens; t.sens.1 = t.sens; #Part 2 - S.sens, Y.sens, w.sens, t.sens load("C:\\Users\\andre\\VIMS\\PhD\\Research\\PredPrey_Simulations\\SimOutputs\\SensitivityOutput_Method2_part2_8-4-14.RData") S.sens.2 = S.sens; Y.sens.2 = Y.sens; R.sens.2 = R.sens; w.sens.2 = w.sens; t.sens.2 = t.sens; #Merge results together library(abind) #needed for abind() S.sens = abind(S.sens.1[,,!is.na(S.sens.1[1,1,])], S.sens.2[,,!is.na(S.sens.2[1,1,])] ) #bind the two arrays after deleting the NA's in each of them. Y.sens = abind(Y.sens.1[,,!is.na(Y.sens.1[1,1,])], Y.sens.2[,,!is.na(Y.sens.2[1,1,])] ) #bind the two arrays after deleting the NA's in each of them. R.sens = abind(R.sens.1[,,!is.na(R.sens.1[1,1,])], R.sens.2[,,!is.na(R.sens.2[1,1,])] ) #bind the two arrays after deleting the NA's in each of them. w.sens = abind(w.sens.1[,,,!is.na(w.sens.1[1,1,1,])], w.sens.2[,,,!is.na(w.sens.2[1,1,1,])] ) #bind the two arrays after deleting the NA's in each of them. t.sens = abind(t.sens.1[,,,!is.na(t.sens.1[1,1,1,])], t.sens.2[,,,!is.na(t.sens.2[1,1,1,])] ) #bind the two arrays after deleting the NA's in each of them. #S.sens[sims,scen,sens]; Y.sens[sims,scen,sens]; w.sens[sims,age,scen,sens]; t.sens[sims,scen,PvNP,sens]; #################################### ## MAKE SUMMARY SENSITIVITY PLOTS #################################### #Function for extracting the right-most characters in a string substrRight <- function(x, n){ substr(as.character(x), nchar(as.character(x))-n+1, nchar(as.character(x))) } ##METHOD 1 #For MC plots, see PredPreySim_graphing_10-5-13.r ##METHOD 2 #Calculate mean and sd of %diff across simulations for SSB #SSB S.sens.mean = apply(S.sens, c(2,3), mean) S.sens.sd = apply(S.sens, c(2,3), sd) S.sens.CV = S.sens.sd / S.sens.mean *100 #present variability as CV (SD/mean) S.sens.anom = S.sens.mean - S.sens.mean[,"BASE."] #calc anomaly of mean %diff #YIELD Y.sens.mean = apply(Y.sens, c(2,3), mean) Y.sens.sd = apply(Y.sens, c(2,3), sd) Y.sens.CV = Y.sens.sd / Y.sens.mean *100 #present variability as CV (SD/mean) Y.sens.anom = Y.sens.mean - Y.sens.mean[,"BASE."] #calc anomaly of mean %diff #RECRUITMENT R.sens.mean = apply(R.sens, c(2,3), mean) R.sens.sd = apply(R.sens, c(2,3), sd) R.sens.CV = R.sens.sd / R.sens.mean *100 #present variability as CV (SD/mean) R.sens.anom = R.sens.mean - R.sens.mean[,"BASE."] #calc anomaly of mean %diff #Weight at age w.sens.mean = apply(w.sens, c(2,3,4), mean) w.sens.sd = apply(w.sens, c(2,3,4), sd) w.sens.CV = w.sens.sd / w.sens.mean *100 #present variability as CV (SD/mean) w.sens.anom = w.sens.mean for (i in 1:length(w.sens.mean[1,1,])) { w.sens.anom[,,i] = w.sens.mean[,,i] - w.sens.mean[,,"BASE."] } #calc anomaly of mean %diff #Time to target (note that the "base" values won't be identical to the simfxn() run because i only did 50 sim runs for the sensitivity analysis! 9-26-13) #With the 50 sim runs, the median delta t for prey 1 pulse is ~2% (as opposed to 0 for the 100 sims presented in the paper) #Consequently, the plot of sensitivities for delta t t.sens.diff = (t.sens[,,"NoPulse",] - t.sens[,,"Pulse",])/t.sens[,,"NoPulse",]*100 #For each run, calculate the % decrease in time to reach target SSB #t.sens.diff = t.sens[,,1,] - t.sens[,,2,] #OLD: For each run, calculate how much faster target SSB is reached (in years) by having the pulse (relative to no pulse) t.sens.diff.base = (t.sens[,,"NoPulse","BASE."] - t.sens[,,"Pulse","BASE."])/t.sens[,,"NoPulse","BASE."]*100 #For each run with default parameters, calculate the % decrease in time to target SSB t.sens.med = apply(t.sens.diff, c(2,3), function(x) median(x,na.rm=T)) t.sens.mean = apply(t.sens.diff, c(2,3), function(x) mean(x,na.rm=T)) #NOTE THAT I"M REMOVING THE NAs t.sens.sd = apply(t.sens.diff, c(2,3), function(x) sd(x,na.rm=T)) t.sens.CV = t.sens.sd / t.sens.mean *100 #present variability as CV (SD/mean) #t.sens.anom = t.sens.mean - t.sens.mean[,"BASE."] #calc anomaly using MEAN %diff (UNITS: %) t.sens.anom = t.sens.med - t.sens.med[,"BASE."] #calc anomaly of MEDIANS of %diff (UNITS: %) t.sens.anom = round(t.sens.med,2) - c(0,5.13,0) #calc anomaly of MEDIANS of %diff (UNITS: %) #Format data for easier plotting #Melt arrays into data frame library(reshape2) S.sens.plot = melt(S.sens.anom, varnames=c("Scenario","Parameter"), value.name="S.anom") #melt array to facilitate plotting by group S.sens.CV.2 = melt(S.sens.CV, varnames=c("Scenario","Parameter"), value.name="S.CV") #melt array to facilitate plotting by group Y.sens.plot = melt(Y.sens.anom, varnames=c("Scenario","Parameter"), value.name="Y.anom") #melt array to facilitate plotting by group Y.sens.CV.2 = melt(Y.sens.CV, varnames=c("Scenario","Parameter"), value.name="Y.CV") #melt array to facilitate plotting by group R.sens.plot = melt(R.sens.anom, varnames=c("Scenario","Parameter"), value.name="R.anom") #melt array to facilitate plotting by group R.sens.CV.2 = melt(R.sens.CV, varnames=c("Scenario","Parameter"), value.name="R.CV") #melt array to facilitate plotting by group w.sens.plot = melt(w.sens.anom, varnames=c("Age","Scenario","Parameter"), value.name="w.anom") #melt array to facilitate plotting by group w.sens.CV.2 = melt(w.sens.CV, varnames=c("Age","Scenario","Parameter"), value.name="w.CV") #melt array to facilitate plotting by group t.sens.plot = melt(t.sens.anom, varnames=c("Scenario","Parameter"), value.name="t.anom") #melt array to facilitate plotting by group t.sens.CV.2 = melt(t.sens.CV, varnames=c("Scenario","Parameter"), value.name="t.CV") #melt array to facilitate plotting by group #combine anom and CV data Sens.plot = merge(S.sens.plot,S.sens.CV.2) #merge anom and CV of SSB data into one data.frame. Sens.plot = merge(Sens.plot,Y.sens.plot) #merge with yield anomalies Sens.plot = merge(Sens.plot,Y.sens.CV.2) #merge with yield CV Sens.plot = merge(Sens.plot,R.sens.plot) #merge with yield anomalies Sens.plot = merge(Sens.plot,R.sens.CV.2) #merge with yield CV Sens.plot = merge(Sens.plot,t.sens.plot) #merge with yield anomalies Sens.plot = merge(Sens.plot,t.sens.CV.2) #merge with yield CV Sens.plot.w = merge(w.sens.plot,w.sens.CV.2) #merge weight at age anomalies and CV (separate dataframe because of the added age information #Organize data frame by specifying different pieces of information Sens.plot$P.scen = substr(Sens.plot$Scenario,1,6) Sens.plot$I.scen = substr(Sens.plot$Scenario,8,10) Sens.plot$F.scen = substrRight(Sens.plot$Scenario,5) Sens.plot$Parm = ifelse(Sens.plot$Parameter=="BASE.", "BASE.", substr(Sens.plot$Parameter,1,nchar(as.character(Sens.plot$Parameter))-3) ) Sens.plot$Shift2 = ifelse(Sens.plot$Parameter=="BASE.", "BASE.", substrRight(Sens.plot$Parameter,2) ) Sens.plot$Shift = ifelse(Sens.plot$Shift2=="BASE.", "BASE.", ifelse(Sens.plot$Shift2=="hi", "+20%", "-20%")) #Rename the parameter shifts for plotting Sens.plot$Shift = factor(Sens.plot$Shift, levels=c("BASE.","+20%","-20%")) Sens.plot$Plot.group = paste(Sens.plot$Parm,"_P",substrRight(Sens.plot$P.scen,1),sep="") Sens.plot$Plot.group2 = paste(Sens.plot$Parm,"_P",substrRight(Sens.plot$P.scen,1),"-",Sens.plot$I.scen,sep="") Sens.plot$Plot.group3 = paste(Sens.plot$Parameter,"_P",substrRight(Sens.plot$P.scen,1),"-",Sens.plot$I.scen,sep="") #Bring in parameter category information for arranging plots #TEMPORARY: write.csv(Sens.plot, "C:\\Users\\andrebuc\\VIMS\\PhD\\Research\\PredPrey_Simulations\\SimOutputs\\Sens.plot_9-26-13.csv") Parm.cat = read.csv("C:\\Users\\andre\\VIMS\\PhD\\Research\\PredPrey_Simulations\\SimOutputs\\ParameterCategories_11-7-13.csv") Sens.plot = merge(Sens.plot,Parm.cat); Sens.plot = Sens.plot[order(Sens.plot$Category,Sens.plot$Parm, Sens.plot$P.scen),] Sens.plot$Order = c(1:nrow(Sens.plot)) #Reorder Plot.group variables to have the factor levels grouped by category Sens.plot$Plot.group = factor(Sens.plot$Plot.group, levels=Sens.plot[!duplicated(Sens.plot$Plot.group),"Plot.group"]) Sens.plot$Plot.group2 = factor(Sens.plot$Plot.group2, levels=Sens.plot[!duplicated(Sens.plot$Plot.group2),"Plot.group2"]) Sens.plot$Plot.group3 = factor(Sens.plot$Plot.group3, levels=Sens.plot[!duplicated(Sens.plot$Plot.group3),"Plot.group3"]) #levels(Sens.plot$Plot.group3) #Sens.plot Sens.plot.w = merge(Sens.plot.w,Sens.plot[,-c(5:8)],by=c("Scenario","Parameter")) #Merge wt-at-age dataframe with the other to pull in all the ancillary variables. ## MAKE PLOTS # Notes on figures # - The red ("hi") bars are drawn behind the "lo" green bars. So, if a red bar is not visible, it means it is zero, or less than the # magnitude of the green bar (in the same positive or negative direction. library(ggplot2) #Specify parameter labels using the greek letters and subscripts as needed. This is for the parameters to match the official table and equations in Chapter 4. 10-5-13) parm.lab = c("atck.infl.1_P1" = expression(~" P1"), "atck.infl.1_P2" = expression(eta[1]~" P2"), "atck.infl.1_P3" = expression(~" P3"), "atck.infl.2a_P1" = expression(~" P1"), "atck.infl.2a_P2" = expression(eta[12]~" P2"), "atck.infl.2a_P3" = expression(~" P3"), "atck.infl.2b_P1" = expression(~" P1"), "atck.infl.2b_P2" = expression(eta[22]~" P2"), "atck.infl.2b_P3" = expression(~" P3"), "atck.infl.3_P1" = expression(~" P1"), "atck.infl.3_P2" = expression(eta[3]~" P2"), "atck.infl.3_P3" = expression(~" P3"), "atck.max.1_P1" = expression(~" P1"), "atck.max.1_P2" = expression(lambda[1]~" P2"), "atck.max.1_P3" = expression(~" P3"), "atck.max.2a_P1" = expression(~" P1"), "atck.max.2a_P2" = expression(lambda[2]~" P2"), "atck.max.2a_P3" = expression(~" P3"), "atck.max.3_P1" = expression(~" P1"), "atck.max.3_P2" = expression(lambda[3]~" P2"), "atck.max.3_P3" = expression(~" P3"), "atck.rho.1_P1" = expression(~" P1"), "atck.rho.1_P2" = expression(rho[1]~" P2"), "atck.rho.1_P3" = expression(~" P3"), "atck.rho.2a_P1" = expression(~" P1"), "atck.rho.2a_P2" = expression(rho[2]~" P2"), "atck.rho.2a_P3" = expression(~" P3"), "atck.rho.3_P1" = expression(~" P1"), "atck.rho.3_P2" = expression(rho[3]~" P2"), "atck.rho.3_P3" = expression(~" P3"), "CA_P1" = expression(~" P1"), "CA_P2" = expression(CA~" P2"), "CA_P3" = expression(~" P3"), "CB_P1" = expression(~" P1"), "CB_P2" = expression(CB~" P2"), "CB_P3" = expression(~" P3"), "CQ_P1" = expression(~" P1"), "CQ_P2" = expression(CQ~" P2"), "CQ_P3" = expression(~" P3"), "CTM_P1" = expression(~" P1"), "CTM_P2" = expression(CTM~" P2"), "CTM_P3" = expression(~" P3"), "CTO_P1" = expression(~" P1"), "CTO_P2" = expression(CTO~" P2"), "CTO_P3" = expression(~" P3"), "KL_P1" = expression(~" P1"), "KL_P2" = expression(KL~" P2"), "KL_P3" = expression(~" P3"), "KR_P1" = expression(~" P1"), "KR_P2" = expression(KR~" P2"), "KR_P3" = expression(~" P3"), "KW_P1" = expression(~" P1"), "KW_P2" = expression(KW~" P2"), "KW_P3" = expression(~" P3"), "M_P1" = expression(~" P1"), "M_P2" = expression(M~" P2"), "M_P3" = expression(~" P3"), "N.prey.2.mean_P1" = " P1", #"N.prey.2.mean_P2" = "N(bar) P2", "N.prey.2.mean_P2" = expression(bar(N)[2]~" P2"), "N.prey.2.mean_P3" = " P3", "N.prey.dev_P1" = expression(~" P1"), "N.prey.dev_P2" = expression(sigma[gamma]~" P2"), "N.prey.dev_P3" = expression(~" P3"), "n.pulse_P1" = expression(~" P1"), "n.pulse_P2" = expression(n[pulse]~" P2"), "n.pulse_P3" = expression(~" P3"), "pulse.range_P1" = expression(~" P1"), "pulse.range_P2" = expression(psi~" P2"), "pulse.range_P3" = expression(~" P3"), "relB.1_P1" = expression(~" P1"), "relB.1_P2" = expression(tau[1]~" P2"), "relB.1_P3" = expression(~" P3"), "relB.3_P1" = expression(~" P1"), "relB.3_P2" = expression(tau[3]~" P2"), "relB.3_P3" = expression(~" P3"), "alpha.SR_P1" = expression(~" P1"), "alpha.SR_P2" = expression(alpha[SR]~" P2"), "alpha.SR_P3" = expression(~" P3"), "beta.SR_P1" = expression(~" P1"), "beta.SR_P2" = expression(beta[SR]~" P2"), "beta.SR_P3" = expression(~" P3"), "K.SR_P1" = expression(~" P1"), "K.SR_P2" = expression(kappa[SR]~" P2"), "K.SR_P3" = expression(~" P3"), "R.dev.sd_P1" = expression(~" P1"), "R.dev.sd_P2" = expression(sigma[delta]~" P2"), "R.dev.sd_P3" = expression(~" P3") ) #SSB (new: 10-5-13) #anomalies S.gplot = ggplot(subset(Sens.plot, Parameter!="BASE." & I.scen!="mix"), aes(x=Plot.group, y=S.anom, fill=Shift)) + scale_x_discrete(labels= parm.lab) + #Use this line to label the X-axis parameters with the greek labels used in the manuscript; Otherwise, comment out. geom_bar(stat="identity", position="identity") + theme(axis.text.x=element_text(angle=90, vjust=0.5, hjust=1 )) + #adjusts the angle, and left/center/right alignment of the text. xlab("") + ylab(expression(paste("Mean difference in ", Delta,"S (%)",sep=""))) # + labs(title="SSB Sensitivity") #Figure title S.gplot #YIELD (new: 10-5-13) #anomalies Y.gplot = ggplot(subset(Sens.plot, Parameter!="BASE." & I.scen!="mix"), aes(x=Plot.group, y=Y.anom, fill=Shift)) + scale_x_discrete(labels= parm.lab) + #Use this line to label the X-axis parameters with the greek labels used in the manuscript; Otherwise, comment out. geom_bar(stat="identity", position="identity") + theme(axis.text.x=element_text(angle=90, vjust=0.5, hjust=1 )) + #adjusts the angle, and left/center/right alignment of the text. xlab("") + ylab(expression(paste("Mean difference in ", Delta,"Y (%)",sep=""))) # + labs(title="SSB Sensitivity") #Figure title Y.gplot #RECRUITMENT (need to re-run sensitivity analysis to get R data!) R.gplot = ggplot(subset(Sens.plot, Parameter!="BASE." & I.scen!="mix"), aes(x=Plot.group, y=R.anom, fill=Shift)) + scale_x_discrete(labels= parm.lab) + #Use this line to label the X-axis parameters with the greek labels used in the manuscript; Otherwise, comment out. geom_bar(stat="identity", position="identity") + theme(axis.text.x=element_text(angle=90, vjust=0.5, hjust=1 )) + #adjusts the angle, and left/center/right alignment of the text. xlab("") + ylab(expression(paste("Mean difference in ", Delta,"R (%)",sep=""))) # + labs(title="SSB Sensitivity") #Figure title R.gplot #TIME TO TARGET (YRS) (new: 10-5-13) #anomalies t.gplot = ggplot(subset(Sens.plot, Parameter!="BASE." & I.scen!="mix"), aes(x=Plot.group, y=t.anom, fill=Shift)) + scale_x_discrete(labels= parm.lab) + #Use this line to label the X-axis parameters with the greek labels used in the manuscript; Otherwise, comment out. geom_bar(stat="identity", position="dodge") + # , position="identity" #Use "dodge" to put bars side by side theme(axis.text.x=element_text(angle=90, vjust=0.5, hjust=1 )) + #adjusts the angle, and left/center/right alignment of the text. xlab("") + ylab(expression(paste("Median difference in ", Delta,"t (%)",sep=""))) # + labs(title="SSB Sensitivity") #Figure title #NOTES: #The large negative values for P2 indicate that these parameter perturbations lead to a median delta t value of zero... t.gplot #WEIGHT AT AGE anomalies (10-8-13) & Age=="age-2" & Age=="age-1" & Parm %in% c("BASE.","atck.max.1","atck.max.1","atck.max.1","n.pulse","pulse.range") #Anomalies w.gplot = ggplot(subset(Sens.plot.w, Parameter!="BASE." & I.scen!="mix" & Age=="age-1" ), aes(x=Plot.group, y=w.anom, fill=Shift)) + # scale_x_discrete(labels= parm.lab) + #Use this line to label the X-axis parameters with the greek labels used in the manuscript; Otherwise, comment out. facet_wrap(~Age, ncol=4) + geom_bar(stat="identity", position="identity") + theme(axis.text.x=element_text(angle=90, vjust=0.5, hjust=1 )) + xlab("") + ylab(expression(paste("Mean difference in ", Delta,"w (%)",sep=""))) #+ labs(title="WEIGHT-AT-AGE Sensitivity" ) w.gplot #NOTES: # - Age-0 - most sensitive to lambda 1 and 2, CA, CB, Nbar, phi, tau 1 (max value ~3%) # - Age-1 - most sensitive to infl 3, CA, CB, Nbar, phi (max value ~3%) # - Age-2 - most sensitive to infl 3, CTO, Nbar, phi (max value ~3%) # - OVERALL - most sensitive to Nbar, phi, inflection 3, CTO, CA, CB (but mostly <3%) ### PLOT ALL SENSITIVITY PLOTS TOGETHER AND EXPORT TO EPS library(gridExtra) #see http://www.r-bloggers.com/extra-extra-get-your-gridextra/ setEPS() #postscript("C:\\Users\\andre\\VIMS\\PhD\\Research\\PredPrey_Simulations\\SimFigures\\SensitivityFigs_8-26-14.eps", height=15,width=15*8.5/11) #pdf("C:\\Users\\andrebuc\\VIMS\\PhD\\Research\\PredPrey_Simulations\\SimFigures\\TEST.pdf", height=15,width=15*8.5/11) grid.arrange(w.gplot, S.gplot, Y.gplot, R.gplot, t.gplot, nrow=5) dev.off() ##OLD PLOTS #SSB #anomalies ggplot(subset(Sens.plot, Parameter!="BASE." & I.scen!="mix"), aes(x=Plot.group, y=S.anom, fill=Shift)) + geom_bar(stat="identity", position="identity") + theme(axis.text.x=element_text(angle=90)) + xlab("") + ylab("Change in mean SSB relative to BASE run (%)") + labs(title="SSB Sensitivity" ) #CV ggplot(subset(Sens.plot, I.scen!="mix"), aes(x=Plot.group3, y=S.CV)) + #, fill=Shift geom_bar(stat="identity", position="identity") + theme(axis.text.x=element_text(angle=90)) + xlab("") + ylab("CV of SSB differences (%)") + labs(title="SSB - CV Sensitivity" ) #Yield anomalies #anomalies ggplot(subset(Sens.plot, Parameter!="BASE."& I.scen!="mix"), aes(x=Plot.group, y=Y.anom, fill=Shift)) + geom_bar(stat="identity", position="identity") + theme(axis.text.x=element_text(angle=90)) + xlab("") + ylab("Change in mean YIELD relative to BASE run (%)") + labs(title="YIELD Sensitivity" ) #CV ggplot(subset(Sens.plot, I.scen!="mix"), aes(x=Plot.group3, y=Y.CV)) + #, fill=Shift geom_bar(stat="identity", position="identity") + theme(axis.text.x=element_text(angle=90)) + xlab("") + ylab("CV of YIELD differences (%)") + labs(title="YIELD - CV Sensitivity" ) #Weight at age anomalies & Age=="age-1" & Age=="age-1" & Parm %in% c("BASE.","atck.max.1","atck.max.1","atck.max.1","n.pulse","pulse.range") #Anomalies ggplot(subset(Sens.plot.w, Parameter!="BASE." & I.scen!="mix" ), aes(x=Plot.group, y=w.anom, fill=Shift)) + # facet_wrap(~Age, ncol=4) + geom_bar(stat="identity", position="identity") + theme(axis.text.x=element_text(angle=90)) + xlab("") + ylab("Change in mean w@age relative to BASE run (%)") + labs(title="WEIGHT-AT-AGE Sensitivity" ) #CV ggplot(subset(Sens.plot.w, I.scen!="mix" ), #Use this to "zoom in" on a few cases: & Age=="age-1" & Parm %in% c("BASE.","atck.max.1","atck.max.1","atck.max.1","n.pulse","pulse.range") aes(x=Plot.group3, y=w.CV)) + #, fill=Shift facet_wrap(~Age, ncol=4) + geom_bar(stat="identity", position="identity") + theme(axis.text.x=element_text(angle=90)) + xlab("") + ylab("CV of WEIGHT-AT-AGE differences (%)") + labs(title="WEIGHT-AT-AGE - CV Sensitivity" ) #Time to target anomalies #Anomalies ggplot(subset(Sens.plot, Parameter!="BASE." & I.scen!="mix" ), aes(x=Plot.group, y=Y.anom, fill=Shift)) + # #facet_wrap(~Age, ncol=4) + geom_bar(stat="identity", position="identity") + theme(axis.text.x=element_text(angle=90)) + xlab("") + ylab(expression(paste("Change in ",~Delta, "t (%)"))) + # "Change in mean time to target relative to BASE run (years)") labs(title=expression(paste(Delta, "t Sensitivity")) ) #CV ggplot(subset(Sens.plot, I.scen!="mix" ), #Use this to "zoom in" on a few cases: & Age=="age-1" & Parm %in% c("BASE.","atck.max.1","atck.max.1","atck.max.1","n.pulse","pulse.range") aes(x=Plot.group3, y=w.CV)) + #, fill=Shift #facet_wrap(~Age, ncol=4) + geom_bar(stat="identity", position="identity") + theme(axis.text.x=element_text(angle=90)) + xlab("") + ylab("CV of WEIGHT-AT-AGE differences (%)") + labs(title="WEIGHT-AT-AGE - CV Sensitivity" ) ##OLD JUNK = c("atck.infl.1_P1" = expression(eta[1]~"_P1"), "atck.infl.1_P2" = expression(eta[1]~"_P2"), "atck.infl.1_P3" = expression(eta[1]~"_P3"), "atck.infl.2a_P1" = expression(eta[12]~"_P1"), "atck.infl.2a_P2" = expression(eta[12]~"_P2"), "atck.infl.2a_P3" = expression(eta[12]~"_P3"), "atck.infl.2b_P1" = expression(eta[22]~"_P1"), "atck.infl.2b_P2" = expression(eta[22]~"_P2"), "atck.infl.2b_P3" = expression(eta[22]~"_P3"), "atck.infl.3_P1" = expression(eta[3]~"_P1"), "atck.infl.3_P2" = expression(eta[3]~"_P2"), "atck.infl.3_P3" = expression(eta[3]~"_P3"), "atck.max.1_P1" = expression(lambda[1]~"_P1"), "atck.max.1_P2" = expression(lambda[1]~"_P2"), "atck.max.1_P3" = expression(lambda[1]~"_P3"), "atck.max.2a_P1" = expression(lambda[2]~"_P1"), "atck.max.2a_P2" = expression(lambda[2]~"_P2"), "atck.max.2a_P3" = expression(lambda[2]~"_P3"), "atck.max.3_P1" = expression(lambda[3]~"_P1"), "atck.max.3_P2" = expression(lambda[3]~"_P2"), "atck.max.3_P3" = expression(lambda[3]~"_P3"), "atck.rho.1_P1" = expression(rho[1]~"_P1"), "atck.rho.1_P2" = expression(rho[1]~"_P2"), "atck.rho.1_P3" = expression(rho[1]~"_P3"), "atck.rho.2a_P1" = expression(rho[2]~"_P1"), "atck.rho.2a_P2" = expression(rho[2]~"_P2"), "atck.rho.2a_P3" = expression(rho[2]~"_P3"), "atck.rho.3_P1" = expression(rho[3]~"_P1"), "atck.rho.3_P2" = expression(rho[3]~"_P2"), "atck.rho.3_P3" = expression(rho[3]~"_P3"), "N.prey.2.mean_P1" = "N(bar)_P1", "N.prey.2.mean_P2" = "N(bar)_P2", "N.prey.2.mean_P3" = "N(bar)_P3", "N.prey.dev_P1" = expression(sigma[gamma]~"_P1"), "N.prey.dev_P2" = expression(sigma[gamma]~"_P2"), "N.prey.dev_P3" = expression(sigma[gamma]~"_P3"), "n.pulse_P1" = expression(n[pulse]~"_P1"), "n.pulse_P2" = expression(n[pulse]~"_P2"), "n.pulse_P3" = expression(n[pulse]~"_P3"), "pulse.range_P1" = expression(psi~"_P1"), "pulse.range_P2" = expression(psi~"_P2"), "pulse.range_P3" = expression(psi~"_P3"), "relB.1_P1" = expression(tau[1]~"_P1"), "relB.1_P2" = expression(tau[1]~"_P2"), "relB.1_P3" = expression(tau[1]~"_P3"), "relB.3_P1" = expression(tau[3]~"_P1"), "relB.3_P2" = expression(tau[3]~"_P2"), "relB.3_P3" = expression(tau[3]~"_P3"), "alpha.SR_P1" = expression(alpha[SR]~"_P1"), "alpha.SR_P2" = expression(alpha[SR]~"_P2"), "alpha.SR_P3" = expression(alpha[SR]~"_P3"), "beta.SR_P1" = expression(beta[SR]~"_P1"), "beta.SR_P2" = expression(beta[SR]~"_P2"), "beta.SR_P3" = expression(beta[SR]~"_P3"), "K.SR_P1" = expression(kappa[SR]~"_P1"), "K.SR_P2" = expression(kappa[SR]~"_P2"), "K.SR_P3" = expression(kappa[SR]~"_P3"), "R.dev.sd_P1" = expression(sigma[delta]~"_P1"), "R.dev.sd_P2" = expression(sigma[delta]~"_P2"), "R.dev.sd_P3" = expression(sigma[delta]~"_P3") ) ################################################## ### EXTRA CODE - PLOTTING OF FUNCTIONAL RESPONSE ################################################## #PLOTS FOR FUNCTIONAL RESPONSE #FxnR[year,season,region,age,prey] #N.prey[year,season,region,prey] #Calculate min and max F for graphing minFx = min(FxnR) maxFx = max(FxnR) par(mfrow=c(2,preys)) #Plot of Fx vs. prey abundance, with a different panel for each prey for (i in 1:preys) { plot(N.prey[,1,1,i],FxnR[,1,1,1,i],ylim=c(minFx,maxFx), type="l",lwd=3, col="black", ylab="Prey Consumed (kg/season)", xlab="Prey Abundance") lines(N.prey[,1,1,i],FxnR[,1,1,2,i], col="blue",lwd=3) lines(N.prey[,1,1,i],FxnR[,1,1,3,i], col="red",lwd=3) lines(N.prey[,1,1,i],FxnR[,1,1,4,i], col="green",lwd=3) lines(N.prey[,1,1,i],FxnR[,1,1,5,i], col="orange",lwd=3) legend("bottomright",c("Age0","Age1","Age2","Age3","Age4"), col=c("black","blue","red","green","orange"),lty=1,lwd=3,merge=T,bg="white") title(paste("Prey",i)) } #Plot of Fx vs. Year for all ages and prey combined plot(c(1:years),FxnR[,1,1,1,1],ylim=c(minFx,maxFx),type="l", lwd=3, lty=1, col="black", ylab="F", xlab="Year") lines(FxnR[,1,1,2,1], col="blue",lwd=3, lty=1) lines(FxnR[,1,1,3,1], col="red",lwd=3, lty=1) lines(FxnR[,1,1,1,2], col="black",lwd=3, lty=2) lines(FxnR[,1,1,2,2], col="blue",lwd=3, lty=2) lines(FxnR[,1,1,3,2], col="red",lwd=3, lty=2) legend("bottomright",c("Prey 1","Prey 2"),lty=c(1,2),lwd=3,merge=T,bg="white") title("F vs. Year") #Plot of relative prey abundances over time plot(c(1:years), N.prey[,1,1,1] / N.prey[,1,1,2] , type="l", lwd=3, xlab="Year", ylab="Relative Prey Abundance (Prey-1 / Prey-2)") title("Relative prey abundances") ################################################# # EXTRA CODE THAT IS NON CRITICAL - Trying to optimize Attack rates ################################################# #Function for optimizing attack rates attackrates = function(x) { #x is a vector of parameter values x.1 = exp(x[1]) #max attackrate for prey 1 NOTE: easier to optimize the log of the parameters because they are so small! x.2 = exp(x[2]) #max attackrate for prey 2 x.3 = x[3] #rate of reaching max attack rate (atck.rho) for prey 1 NOTE: no need to exp x.4 = x[4] #rate of reaching max attack rate (atck.rho) for prey 2 #Run the growth function to get weight at age given the attack rate parameters defined by x growthfxn(atck.max.1=x.1,atck.max.2=x.2, atck.rho.1=x.3, atck.rho.2=x.4) #UPDATE THIS AS NEEDED!!! (either have these as parameters to optimize, or fix them at certain values) #Calculate Residual sum of squares for weight at age as objective function to estimate attack rates dim(w.pred) #[years,seasons,regions,ages] dim(w) #[seasons,ages] head(LWdata) RSS = sum((w.pred[20,,1,] - w)^2) return(RSS) } #end attackrates() ### TEMPORARY JUNK - TRYING TO OPTIMIZE ALPHA'S FOR WEIGHT-AT-AGE growthfxn() attackrates(c(-21.9,-17)) attackrates(c(-21.9,-19,0,0.00000003)) #6/25 - with atck.slope attackrates(c(-152,-152,170,-)) #6/25 - with atck.slope #Optimization test = optim(c(0.000000001,0.00000001),attackrates,method="BFGS",hessian=T) test = optim(c(-12,-10),attackrates,method="BFGS",hessian=T) test = optim(c(-22,-19,-17),attackrates,method="BFGS",hessian=T) #6/25 - function now has atck rate slope, but assumed slope=0 for prey 1 test = optim(c(-22,-19,-20,-17),attackrates,method="BFGS",hessian=T) #6/25 - function now has atck rate slope, but assumed slope=0 for prey 1 test.logistic = optim(c(-20,-20,0.5,0.5),attackrates,method="BFGS",hessian=T) #-19.3, -19, -3, -0.44 test.logistic = optim(c(-19,-19),attackrates,method="BFGS",hessian=T) #With rho1=-1,rho2=0.5, optim gave max1=-28.27, max2=-15.7 test.logistic = optim(c(-19,-19,-1,0.5),attackrates,method="L-BFGS-B",lower=c(-25,-25,-3,0),upper=c(-10,-10,0,3)) #With rho1=-1,rho2=0.5, optim gave max1=-28.27, max2=-15.7 growthfxn(atck1=exp(-21.9),atck2=exp(-17), #Pre "attack rate slope" addition to the model Cmax.perc=0.05*91, ce.val=0.1, N.prey.1.cv = .8, N.prey.2.cv = .8 ) growthfxn(atck1=exp(-32.4),atck2=exp(-27.8),atck2.slope=exp(-16.3), #Optimized parameters for atck1,atck2, and atck2.slope (assume atck1.slope=0) Cmax.perc=0.05*91, ce.val=0.1, N.prey.dev=0.4) growthfxn(atck1=exp(-30.7),atck2=exp(-26.1), #6/26/12 - Optimized parameters for atck1,atck2, and both slopes (but these could potentially be a local minimum. atck1.slope=exp(-21.4), atck2.slope=exp(-17.5), Cmax.perc=0.05*91, ce.val=0.1, N.prey.dev=0.4) atck #Notes: # - 6/25/12 - After changing the growth function to the one that Wilberg uses (relying on Gmax), # i was able to obtain an optimized value for the attack rates on each of the prey. This was # done by assuming prey abundances for each prey population (from Ches Bay), and temporarily # assuming that attack rate remains constant across ages. It was also critical to optimize the log # of the actual values because they were so small (e.g. 1E-14). But, the consumption at age relative # to Cmax was very low, due to the constact attack rate assumption. It is more realistic to have # the attack rates change over time, and increase for fish prey as piscivory and capture efficiency # increase. # - 6/27/12 - Initial attempts at making the attack rate of function of age involved making it a linear # function of age. This caused problems when attack rates went negative. Also, the "optimized" values for the # two attack rates and slopes represented a local minimum, and it caused consumption at age-0 to be essentially zero # (ie there was no growth for age-0 fish. This is a problem that needs to be fixed. One potential solution is to recode # attack rates as a logistic function with age. # - 7/10/12 - The logistic function was applied. After spending a lot of time trying to "optimize" the attack rates, i decided # to pick reasonable values and move on. The effect of the attack rate assumptions can be evaluated later through sensitivity # analyses or multiple scenarios. ###CODE FROM PAT AS LOOPING EXAMPLE### for (y in 1:nrow(N2.Ma)) { if (y==1) #initiate population in first year { #initiate age-1s as a function of virgin recruitment N2.Fa[y,1]=R0/2; N2.Ma[y,1]=R0/2 #then other ages in year 1 following mortality equation for (j in 2:(length(ages)-1)) { N2.Ma[y,j]=N2.Ma[y,j-1]*exp(-(M.Ma2[j-1]+testF.M[j])) N2.Fa[y,j]=N2.Fa[y,j-1]*exp(-(M.Fa2[j-1]+testF.F[j])) } #then the plus group....CHECK THIS! N2.Ma[y,max(ages)]=N2.Ma[y,max(ages)-1]*(exp(-(M.Ma2[max(ages)-1]+ testF.M[j]))/(1-exp(-(M.Ma2[max(ages)]+testF.M[j])))) N2.Fa[y,max(ages)]=N2.Fa[y,max(ages)-1]*(exp(-(M.Fa2[max(ages)-1]+ testF.F[j]))/(1-exp(-(M.Fa2[max(ages)]+testF.F[j])))) #Biomass at each age is then related to weight at age B2.Ma[y,]=N2.Ma[y,]*Wa.M B2.Fa[y,]=N2.Fa[y,]*Wa.F #Spawning biomass SB2[y]=sum(B2.Ma[y,]*Mat.Ma+B2.Fa[y,]*Mat.Fa) } #Now for remaining years if (y>1) { for (j in 1:max(ages)) { if (j==1) { #for avg relationship, exclude any recruitment deviations rec=((0.8*R0*Steep*SB2[y-1])/(0.2*phi0*R0*(1-Steep)+(Steep-0.2) *SB2[y-1])) N2.Ma[y,j]=0.5*rec N2.Fa[y,j]=0.5*rec } if (j>1 & j