### OpenBUGS model code ### # Created by Jeff Moore, last edited 2 March 2014 # Code for estimating parameters in JOINT MODEL (detection-groupSize-counts) # # Detection model detail: # Data are detection distances and associated covariates for sperm whales # ESW for small groups (1,2) set to 1/W. For larger groups (>2), detection function is estimated using an intercept-only model # # Group size model detail: # count data 'cs' are cs = 1, 2,... (with some non-integers since each datum is an avg from multiple observers) # Mean group size of small groups fixed at 1.14. Mean group size of larger group sizes in each j,t is modeled # model 'gs' as generalized Poisson variable with overdispersion parameter fixed at alpha=0.12 so that variance ests are similar to variance of data # # Number of groups detected are assumed to be generalized Poisson distributed with constant variance inflation (estimated from separate bootstrap exercise) # # Process model is exponential growth, of regression or Markov form ##################### model{ ###### # PROCESS MODELS FOR ABUNDANCE IN CALIFORNIA CURRENT ###### #### MODEL SET 1: SIMPLE REGRESSION MODEL FOR EXPONENTIAL GROWTH #for (t in 1:allYrs){ #N[t] <- N0*exp(r*t)*exp(err.p.n[t]) #err.p.n[t] ~ dnorm(0,tau.pe.n) #} #### MODEL SET 2: MARKOV MODEL FOR EXPONENTIAL GROWTH N[1] <- N0*exp(r)*exp(err.p.n[1]) err.p.n[1] ~ dnorm(0,tau.pe.n) for(t in 2:allYrs){ N[t] <- N[t-1]*exp(r)*exp(err.p.n[t]) err.p.n[t] ~ dnorm(0,tau.pe.n) } #### DERIVED PARAMETERS FOR WHOLE STUDY AREA for(t in 1:allYrs){D[t] <- N[t]/sum(A[])} # Mean density each year for CA Current study area for(t in 1:nYrs){ N.surveyYrs[t] <- N[T[t]] } # CA Current abundance for study years only #### DISTRIBUTION OF N ACROSS SMALL VS. LARGE GROUPS (SMALL GROUP SIZE = 1 OR 2 INDIVIDUALS) AND ACROSS STRATA, IN SURVEY YEARS for (t in 1:nYrs){ N.small[t] <- N.surveyYrs[t] * f.small[t] N.large[t] <- N.surveyYrs[t] * f.large[t] logit.fsmall[t] ~ dnorm(mu.fsmall, tau.fsmall) f.small[t] <- exp(logit.fsmall[t])/(1+exp(logit.fsmall[t])) # fraction of total abundance occurring in small groups f.large[t] <- 1 - f.small[t] # fraction occurring in large groups for(j in 1:4){ N.strata[t,j] <- sum(N.strata.gs[t,j,]) D.strata[t,j] <- N.strata[t,j] / A[j] # Calculate stratum densities # stratum abundance and density by group size class N.strata.gs[t,j,1] <- N.small[t] * f.strata.size[1,t,j] N.strata.gs[t,j,2] <- N.large[t] * f.strata.size[2,t,j] D.strata.gs[t,j,1] <- N.strata.gs[t,j,1] / A[j] D.strata.gs[t,j,2] <- N.strata.gs[t,j,2] / A[j] } f.strata.size[1,t,1:4] ~ ddirich(alpha.dirich.sm[1:4]) f.strata.size[2,t,1:4] ~ ddirich(alpha.dirich.lg[1:4]) } #### PRIORS FOR PROCESS MODEL N0 <- exp(logN0) # Back-projected population in 1990 (year before 1st survey) logN0 ~ dunif(4,10) # range on real scale is 50 to 20,000 r ~ dunif(-1, 1) # growth rate tau.pe.n <- 1/var.pe.n var.pe.n <- sd.pe.n * sd.pe.n # process variance for N[t], log scale sd.pe.n ~ dunif(0,3) for(j in 1:4){ # dirichlet parms for fractional distribution of abundance across strata alpha.dirich.sm[j] ~ dunif(0,100) # for small-group abundance alpha.dirich.lg[j] ~ dunif(0,100) # for large-group abundance } mu.fsmall ~ dnorm(0,0.00001) # mean logit(proportion of N that occurs in small groups) tau.fsmall <- 1/(sd.fsmall * sd.fsmall) sd.fsmall ~ dunif(0,10) # variance across years for proportion of N occurring in small groups #################### ###### # OBSERVATION MODEL ###### for(t in 1:nYrs){ # for the six survey years for(j in 1:4){ # for the four strata for(g in 1:2){ # for large and small groups #### MODEL FOR EXPECTED NUMBER OF GROUPS # Assume, within a group size class, that group size is independent of the detection process) mu[j,t,g] <- 2*L[j,t]*W * p.avg[t,j,g] * D.strata.gs[t,j,g]/es[g,j,t] #### LIKELIHOOD EXPRESSIONS # Generalized Poisson [Famoye 1993, 2004] with variable alpha[j,t] to maintain constant overdispersion (i.e., assume all j,t have same variance:mean ratio) # Likelihood logL.gp.n[j,t,g] <- n[g,j,t]*log(mu[j,t,g]/(1+alpha.n[j,t,g]*mu[j,t,g])) + (n[g,j,t]-1)*log(1+alpha.n[j,t,g]*n[g,j,t]) - logfact(n[g,j,t]) - mu[j,t,g]*(1+alpha.n[j,t,g]*n[g,j,t])/(1+alpha.n[j,t,g]*mu[j,t,g]) # overdispersion parameter alpha.n[j,t,g] <- (pow(vif[g],0.5)-1) / mu[j,t,g] # Derived from eqn for variance of the generalized poisson: var/mu = vif = (1+alpha*mu)^2. Find alpha[j,t] so that var[j,t]/mu[j,t] = vif # zeros trick nlogL.gp.n[j,t,g] <- -1 * logL.gp.n[j,t,g] zeros.n[j,t,g] <- 0 zeros.n[j,t,g] ~ dpois(nlogL.gp.n[j,t,g]) }}} # Variance inflaction factor, used in the Generalized Poisson approach (estimated using separate bootstrap approach as described in paper supplement) vif[1] <- 1.3 # small groups vif[2] <- 2.3 # large groups ###### DETECTION PROBABILITY COMPONENT OF OBSERVATION MODEL for(t in 1:nYrs){ for(j in 1:4){ for(g in 1:2){ # g1 = small groups, g2 = large groups #avg detectability in j,t is weighted by amount of effort in each Beaufort class p.avg[t,j,g] <- 1/(W*f0.size[g]*L[j,t]) * (Lb[1,j,t]*g0[1]+Lb[2,j,t]*g0[2]+Lb[3,j,t]*g0[3]+Lb[4,j,t]*g0[4]+Lb[5,j,t]*g0[5]) }}} f0.size[1] <- mean.f0.sm # f0 for small groups (null model fit to detection data, see below) f0.size[2] <- 1/W # f0 for large groups, assumed equal to 1/W (i.e., esw = W for large groups) ##### g0 COMPONENT ##### log.g0slope ~ dnorm(-0.125,g0slope.tau)C(,0) # slope (log scale) of declining g0 with beaufort (from Jay Barlow's g0 analysis, unpublished) g0slope.tau <- 1/(.1*.1) # SD = 0.1 on slope gives CVs on beaufort-specific g0 consistent with Barlow's jacknife analysis g0[1] ~ dbeta(15.179, 2.268) # Beaufort 0 & 1 => mean 0.87, cv 0.09 (sd 0.08) g0[2] <- g0[1] * exp(log.g0slope * 2) # Beaufort 2 g0[3] <- g0[1] * exp(log.g0slope * 3) # Beaufort 3 g0[4] <- g0[1] * exp(log.g0slope * 4) # Beaufort 4 g0[5] <- g0[1] * exp(log.g0slope * 5) # Beaufort 5 # Notes for g0 # Barlow's analysis provides the following relative g0 values for sperm whale (relative to g0 in Bft 0): # bft 1 = 0.883, bft 2 = 0.78, bft 3 = 0.688, bft 4 = 0.608, bft 5 = 0.536 # Another analysis of sperm whale g0 (by Barlow and Rankin) provides absolute estimate of 0.63 for visual observers # avg Beaufort conditions for sightings in that study was around 3.5 to 4 (pers comm, Jay Barlow) # If visual g0 in beaufort 3.5 to 4 is 0.63, and if g0 in beaufort 3/4 is 0.69/0.61 of g0 in Beaufort 0, then g0 in beaufort 0 is around 0.9 - 1. # Barlow & Taylor (2005) and Barlow & Forney (2007) use value of 0.87 (CV 0.09) in Beaufort 0/1 # Approach here: Use 0.87 (CV 0.09) for Beaufort 0 and 1. Use estimated slope-decline to estimate g0 in Beaufort 2 - 5 ##### f0 COMPONENT ##### ### Constants pi <- 3.141593 a <- 0.147 # approximation constant used in the error function (erf) for estimating esw (Winitski 2008) ### PRIORS for fixed effects on half-normal detection parameter sigma b.df.0 ~ dunif(-10,4) ### Model detection function separately for large (group size > 2) and small groups ## SMALL GROUPS ONLY: Estimate detection as function of covariates for (i in 1:nSmGrps){ # Covariate models # null (use this for sperm whales -- data don't support fitting of more complex models) mu.df[i] <- b.df.0 # estimate sd and var, given model coeffs above sig.df[i] <- exp(mu.df[i]) sig2.df[i] <- sig.df[i]*sig.df[i] # estimate effective strip width (esw, which is also the integral from 0 to W of g(y)) x[i] <- W/sqrt(2*sig2.df[i]) # argument of the 'error func' in the integrand for detect function g(y) erf[i] <- sqrt(1-exp(-x[i]*x[i]*(4/pi + a*x[i]*x[i])/(1 + a*x[i]*x[i]))) # approximation 1 of the 'error function' esw.sm[i] <- sqrt(pi * sig2.df[i] / 2) * erf[i] f0.sm[i] <- 1/esw.sm[i] # LIKELIHOOD, truncated half-normal likelihood y[smGrp.idx[i]] ~ dnorm(0, tau.df[i]) T(0,W) tau.df[i] <- 1/sig2.df[i] } # end i loop mean.esw.sm <- mean(esw.sm[]) mean.f0.sm <- mean(f0.sm[]) ##### GROUP SIZE COMPONENT OF OBSERVATION MODEL ### COVARIATE MODELS (for group size when distance = 0) for(j in 1:nStrata){ for(t in 1:nYrs){ ## Small groups es[1,j,t] <- 1.141 # mean group size for groups of size <= 2 (SE for small groups is 0.05 -- ignore) ## Large groups es[2,j,t] <- mu.cs[j,t] + 2 # mean group size for group with size > 2 # Covariate models for true group size (minus 2), large groups mu.cs[j,t] <- exp(B.cs.0 + b.cs.jt[j,t]) # random j x t # Model for observed group size mu.cs.m[j,t] <- mu.cs[j,t] * exp(b.cs.m * m[t]) }} ### LIKELIHOOD, Negative Binomial (on transformed group sizes, large groups) for (i in 1:nLgGrps){ ### Generalized Poisson approach # Likelihood logL.s[i] <- gs[lgGrp.idx[i]]*log(mu.cs.m[covars[lgGrp.idx[i],1], covars[lgGrp.idx[i],8]]/(1+alpha.s*mu.cs.m[covars[lgGrp.idx[i],1], covars[lgGrp.idx[i],8]])) + (gs[lgGrp.idx[i]]-1)*log(1+alpha.s*gs[lgGrp.idx[i]]) - logfact(gs[lgGrp.idx[i]]) - mu.cs.m[covars[lgGrp.idx[i],1], covars[lgGrp.idx[i],8]]*(1+alpha.s*gs[lgGrp.idx[i]])/(1+alpha.s*mu.cs.m[covars[lgGrp.idx[i],1], covars[lgGrp.idx[i],8]]) # zeros trick nlogL.s[i] <- -1 * logL.s[i] zeros.s[i] <- 0 zeros.s[i] ~ dpois(nlogL.s[i]) } ### PRIORS for group size component of model # Intercept B.cs.0 ~ dunif(-5,5) # random effect for each j,t sig.cs.jt ~ dunif(0.01,4) tau.cs.jt <- 1/(sig.cs.jt * sig.cs.jt) for(j in 1:nStrata){for (t in 1:nYrs){b.cs.jt[j,t] ~ dnorm(0,tau.cs.jt)}} # year (era) effect (2 eras) b.cs.m ~ dunif(-10,0) # overdispersion alpha.s <- 0.12 # Empirical estimate. VIF (var/mean) of raw group size data corresponds to alpha = 0.13 first 3 years and 0.11 for later 3 years. } # end model