1 / 123

Welcome to our second (advanced) course on using AD Model Builder

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

jag
Download Presentation

Welcome to our second (advanced) course on using AD Model Builder

An Image/Link below is provided (as is) to download presentation Download Policy: Content on the Website is provided to you AS IS for your information and personal use and may not be sold / licensed / shared on other websites without getting consent from its author. Content is provided to you AS IS for your information and personal use only. Download presentation by click this link. While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server. During download, if you can't get a presentation, the file might be deleted by the publisher.

E N D

Presentation Transcript


  1. 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

  2. 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

  3. 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

  4. 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

  5. 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

  6. 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

  7. if condition is true then run this code Advanced Conditional Statements • Runs code if conditions met if (condition) { . . . . . . ; }

  8. 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

  9. Parameterization Issues • How do you estimate highly correlated parameters? • Catchabilities for multiple fisheries • Annual recruitments • Deviation method • Difference method • Random walk

  10. Deviation Method • Estimate one free parameter X • Estimate m parameters as bounded_dev_vector w1, . . . , wm • Then logX1 = logX + w1 . . . . logXn = logX + wn

  11. 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

  12. Difference Method • Estimate n free parameters: X1, y1, . . . , yn-1 • Then logX1 logX2 = logX1 + y1 . . . . . logXn = logXn-1 + yn-1

  13. Random Walk • Estimate n parameters X1, w1, . . . , wn-1 • Then

  14. CAA estimates of population dynamics CAA predictions of observed data Overview of Catch-At-Age Observed data Negative log likelihood

  15. Observed Data • Total annual fishery catch • Proportion of catch-at-age • Auxiliary data • Fishing effort • Survey index of relative abundance

  16. Initial numbers at age Recruitment Population Submodel . . . . logR+w1 logN1,1+y1 logN1,n-1+yn-1 logR+w2 . . . . logR+wm

  17. Numbers of fish Survival Total mortality Fishing mortality Natural mortality Population Submodel

  18. Selectivity Effort Catchability Effort error Population Submodel

  19. Observation Submodel • Baranov’s catch equation

  20. Total catch Observed total catch Observation error Observation Submodel

  21. Proportion of catch-at-age Numbers sampled at age Proportions Effective sample size Obs. proportion of catch-at-age Observation Submodel

  22. Negative Log Likelihoodfor Multinomial

  23. Model Parameters R w1, . . . , wm y1, . . . , yn-1 q s1, . . . , sn-1 z1, . . . , zm se Ratio of relative variances (assumed known)

  24. Negative Log Likelihood (ignoring constants)

  25. Now to look at the code…..

  26. 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

  27. 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

  28. 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

  29. Specify starting values in tpl file INITIALIZATION_SECTION log_q -1.0 • Must recompile tpl file everytime starting values are changed

  30. 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

  31. 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

  32. 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

  33. 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:

  34. 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

  35. 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

  36. 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

  37. 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!

  38. 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.

  39. When g() = 1-- i.e, we are interested in the distribution of a parameter -- ADMB approximates the marginal distribution of1 : with

  40. 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:

  41. 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

  42. Let’s try it out…..

  43. 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.

  44. 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

  45. 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.

  46. The Basic Idea Behind MCMC • Transition fromitoi+1is stochastic • Transition only depends uponi • For suitable transition rule approaches target pdf • Approximate marginal pdf by output of chain after burn-in 100000 1 0 2 ....

  47. Reversibility =s t(r,s) t(s,r) =r (r)t(r,s)= (s)t(s,r)

  48. 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 puti+1= i •  = (*)/ (i)

More Related