1.24k likes | 1.42k Vues
Welcome to our second (advanced) course on using AD Model Builder. Instructors Brian Linton and Travis Brenden Special thanks to Jim Bence, Jenny Nelson, and Weihai Liu And to MSU, GLFC, MDNR and CLC Agency Partners. Statistical Catch-At-Age Analysis. New ADMB concepts and techniques
E N D
Welcome to our second (advanced) course on using AD Model Builder • Instructors Brian Linton and Travis Brenden • Special thanks to Jim Bence, Jenny Nelson, and Weihai Liu • And to MSU, GLFC, MDNR and CLC Agency Partners
Statistical Catch-At-Age Analysis • New ADMB concepts and techniques • Array and matrix functions • Phases • Advanced conditional statements • Some parameterization issues • Overview of catch-at-age model • Population submodel • Observation submodel • Negative log likelihood
Array and Matrix Functions • The operator * provides matrix multiplication • For vector objects x and y, and number z z=x*y; // returns z=x1y1+ . . . +xnyn • For matrix objects x, y and z z=x*y; //returns zi,j=xi,1y1,j+ . . . +xi,nyn,j
Array and Matrix Functions • Functions elem_prod and elem_div provide elementwise multiplication and division • For vector objects x, y and z z=elem_prod(x,y); //returns zi=xiyi z=elem_div(x,y); //returns zi=xi/yi • For matrix objects x, y and z z=elem_prod(x,y); //returns zi,j=xi,jyi,j z=elem_div(x,y); //returns zi,j=xi,j/yi,j
Phases • Minimization of objective function can be carried out in phases • Parameter remains fixed at starting value until its phase is reached, then it become active • Allows difficult parameters to be estimated when other parameters are “almost” estimated
Phases • Specified in PARAMETER_SECTION init_number x //estimate in phase 1 init_number x(1) //estimate in phase 1 init_number x(-1) //remains fixed init_vector x(1,n,2) //estimate in phase 2 init_matrix x(1,n,1,m,3) //estimate in phase 3
if condition is true then run this code Advanced Conditional Statements • Runs code if conditions met if (condition) { . . . . . . ; }
Advanced Conditional Statements • active(parameter) • Returns true if parameter is being estimated • last_phase() • Returns true if in last phase of estimation • mceval_phase() • Returns true if –mceval switch is used • sd_phase() • Returns true if in SD report phase
Parameterization Issues • How do you estimate highly correlated parameters? • Catchabilities for multiple fisheries • Annual recruitments • Deviation method • Difference method • Random walk
Deviation Method • Estimate one free parameter X • Estimate m parameters as bounded_dev_vector w1, . . . , wm • Then logX1 = logX + w1 . . . . logXn = logX + wn
bounded_dev_vector • Specified in PARAMETER_SECTION init_bounded_dev_vector x(1,m,-10,10) • Each element must take value between lower and upper bounds -10 < xi < 10 • All elements must sum to zero
Difference Method • Estimate n free parameters: X1, y1, . . . , yn-1 • Then logX1 logX2 = logX1 + y1 . . . . . logXn = logXn-1 + yn-1
Random Walk • Estimate n parameters X1, w1, . . . , wn-1 • Then
CAA estimates of population dynamics CAA predictions of observed data Overview of Catch-At-Age Observed data Negative log likelihood
Observed Data • Total annual fishery catch • Proportion of catch-at-age • Auxiliary data • Fishing effort • Survey index of relative abundance
Initial numbers at age Recruitment Population Submodel . . . . logR+w1 logN1,1+y1 logN1,n-1+yn-1 logR+w2 . . . . logR+wm
Numbers of fish Survival Total mortality Fishing mortality Natural mortality Population Submodel
Selectivity Effort Catchability Effort error Population Submodel
Observation Submodel • Baranov’s catch equation
Total catch Observed total catch Observation error Observation Submodel
Proportion of catch-at-age Numbers sampled at age Proportions Effective sample size Obs. proportion of catch-at-age Observation Submodel
Model Parameters R w1, . . . , wm y1, . . . , yn-1 q s1, . . . , sn-1 z1, . . . , zm se Ratio of relative variances (assumed known)
Sensitivity to Parameter Starting Values • Why do we care about sensitivity to parameter starting values? • Methods for specifying starting values • Default values • In tpl file • In dat file • In pin file • Precedence between the methods
Why do we care about sensitivity to starting values? • Avoiding local minimums in the likelihood surface • If different starting values lead to solution with lower obj. function value, then you were at a local minimum • Identifying sensitive parameters • If small change to parameter starting value causes large change in results, then you may want to reparameterize model
Default Starting Values • Parameter with unspecified starting value has default starting value of zero • Bounded parameter has default starting value which is midway between lower and upper bounds
Specify starting values in tpl file INITIALIZATION_SECTION log_q -1.0 • Must recompile tpl file everytime starting values are changed
Specify starting values in dat file DATA_SECTION init_number start_log_q PRELIMINARY_CALCS_SECTION log_q = start_log_q; • Can change starting values without recompiling tpl file
Specify starting values in pin file #Example pin file for model with 15 parameters 0 0 0 -1 0 0 0 0 0 0 0 0 0 0 0 • Can change starting values without recompiling tpl file • Must specify a starting value for each parameter
Precedence Between Methods • Specifying starting values in dat file takes precedence over pin file and INITIALIZATION_SECTION • Specifying starting values in pin file takes precedence over INITIALIZATION_SECTION
Asymptotic results • Done automatically for parameters ADMB or for any calculating quantity of type sdreport_* or likeprof_number • Results are in *.std and *.cor • These are based on:
Asymptotic standard errors can produce misleading inferences • When sample sizes are small • and the curvature of the likelihood surface changes substantially within the range of plausible estimates – i.e., “near to the maximum likelihood estimates • We will now explore alternative approaches to assessing uncertainty using the catch-at-age model • Likelihood profile method • Bayesian inference
How to use the profile method • Declare a variable you would like to profile as type likeprof_number in the parameter section, and assign it the correct value in the procedure section. • When you run your program use the lprof switch: myprog -lprof • Results are saved in myvar.plt where myvar is the name of your likeprof_number variable • Your variable is varied over a “profile” of values and the best fit constrained to match each value of your variable is found
PLT File contains list of point (x,y) x is value (say biomass) y is associated prob density Plot of Y vs X gives picture of prob distribution ADMB manual says estimate probability x in (xr,xs) by
Profile likelihood options Switch -prsave This saves the parameter values associated with each step of profile in myvar.pv Options set in tpl (preliminary calcs section): e.g., for lprof var myvar: PRELIMINARY_CALCS_SECTION myvar.set_stepnumber(10); // default is 8 myvar.set_stepsize(0.2); //default is 0.5 Note manuals says stepsize is in estimated standard deviations but this appears to be altered adaptively during the profile WARNING -- LOTS OF STEPS CAN TAKE LOTS OF TIME!
Profile Likelihood Method • This is NOT inverting a likelihood ratio test in ADMB land! • This is Bayesian in philosophy (in the same way that MCMC is). Can also be motivated by likelihood theory (support intervals) • Idea is to use the profile for g() to approximate the probability density function for g.
When g() = 1-- i.e, we are interested in the distribution of a parameter -- ADMB approximates the marginal distribution of1 : with
More generally g() (say biomass) is a complex function of many parameters and we want a pdf for g() • The approximation needs to be modified: • biomass (or other derived quantities) will change at different rates with respect to changes in parameters in different parts of the parameter space. • i.e., some “biomasses” might represent a small part of the parameter space and others might represent a larger part. • The modified approximation is:
PLT File contains list of point (x,y) x is value (say biomass) y is associated prob density Plot of Y vs X gives picture of prob distribution ADMB manual says estimate probability x in (xr,xs) by
MCMC • MCMC is a way to generate samples from a complex multivariate pdf. • In practice the pdf is usually the posterior in a Bayesian analysis. • This is useful in looking at marginal distributions of derived quantities. • These marginal distributions are the same thing the profile likelihood method was approximating.
We need priors to do this! Prior Likelihood Posterior ADMB presumes we are going to start by finding the parameters that maximize the posterior density (called highest posterior or modal estimates), so just miminize the log posterior. Just like a negative log-likelihood but with new terms for priors
Two examples • If prior on M were log-normal with median of 0.2 and with sd for ln(M)=0.1, then just add to your likelihood: • For special case of diffuse prior ln(p()) is constant inside the bounds, so a bounded diffuse prior can be specified just by setting bounds on parameters.
The Basic Idea Behind MCMC • Transition fromitoi+1is stochastic • Transition only depends uponi • For suitable transition rule approaches target pdf • Approximate marginal pdf by output of chain after burn-in 100000 1 0 2 ....
Reversibility =s t(r,s) t(s,r) =r (r)t(r,s)= (s)t(s,r)
Approach to Achieve Reversibility (Metropolis algorithm) • Generate a trial value of * using “candidate” probability distribution • Candidate distribution is easy to simulate • Here we assume candidate distributions propose symmetric transitions • Calculate (i) and (*) • If (*) > (i) make the transition i+1= * • If (*) < (i) generate u~uniform(0,1) • If u< make the transition i+1= * • otherwise stay puti+1= i • = (*)/ (i)