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

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

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

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

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

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

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

44. Let’s try it out…..

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

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

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

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

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

50. 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)