model { # PRIORS and fixed values ------------------------------------------------------------ #...residual errors for pop observations tau.mlp ~ dgamma(0.001,0.001) # May LP residual error tau.m ~ dgamma(0.001,0.001) # May count residual error tau.n ~ dgamma(0.001,0.001) # Nov count residual error #...November count bias #......random-year effect # mu.bias ~ dnorm(lbiasmu,lbiastau) # sigma.bias ~ dunif(0,0.2) # tau.bias = 1/sigma.bias^2 # for (t in 1:nyears) { # epsilon[t] ~ dnorm(0,tau.bias) # log(bias[t]) = mu.bias + epsilon[t] # } #......fixed-year effect # for (t in 1:nyears) { bias[t] ~ dnorm(0.87,1/0.05^2) } # from random-year model: conv. issues for (t in 1:nyears) { bias[t] ~ dunif(0.2,1.8) } #...overdispersion for reproductivity betabinomial #......year-specfic values for (t in 1:nyears) { size[t] ~ dgamma(muS[t]^2/varS[t],muS[t]/varS[t]) } #...differential vulnerability mud <- 2.22 # weighted mean & SS from 1994-96,2008-09,2012-17 (Kevin) sigmad <- 0.033 # SS dv ~ dgamma(mud^2/sigmad,mud/sigmad) #...demographic rates for (t in 1:nyears) { theta[t] ~ dbeta(49,7) hn[t] ~ dunif(0,0.3) # harvest rate of adults in Norway hd[t] ~ dunif(0,0.3) # harvest rate of adults in Denmark ho[t] ~ dunif(0,0.3) # harvest rate of adults in Denmark, Sept-Oct s[t] ~ dbeta(B.par.a[t],B.par.b[t]) # annual survival CMR priors } #...intital May population size Nm[1] ~ dlnorm(lmean.m,ltau.m) log.nm[1] <- log(Nm[1]) #...regression coefficients for reproduction = f(thaw days) beta0 ~ dnorm(0,0.001) beta1 ~ dnorm(0,0.001) # System process ----------------------------------------------------------------------------------- for (t in 2:(nyears+1)) { Nm[t] <- max(0.01,Nm[t-1] * (s[t-1] + r[t-1]*theta[t-1]*(1-dv*hn[t-1]-dv*hd[t-1]))) } for (t in 1:nyears) { Nn[t] <- max(0.01,Nm[t] * pow(theta[t],6/12) * (1-(hn[t]+ho[t]) + r[t]*(1-dv*hn[t]-dv*ho[t]))) } # Likelihood for pops ------------------------------------------------------------------------------ #... May LP population estimates [missing 1999] ################### for (t in 2:7) { mlp[t] ~ dlnorm(log(Nm[t])-0.5/tau.mlp,tau.mlp) # Post predictive check exp.mlp[t] = Nm[t] mlp_pred[t] ~ dlnorm(log(Nm[t])-0.5/tau.mlp,tau.mlp) mlp_fit_pred[t] = pow(pow(mlp_pred[t],0.5) - pow(exp.mlp[t],0.5),2) mlp_fit[t] = pow(pow(mlp[t],0.5) - pow(exp.mlp[t],0.5),2) } for (t in 9:(nyears)) { mlp[t] ~ dlnorm(log(Nm[t])-0.5/tau.mlp,tau.mlp) # Post predictive check exp.mlp[t] = Nm[t] mlp_pred[t] ~ dlnorm(log(Nm[t])-0.5/tau.mlp,tau.mlp) mlp_fit[t] = pow(pow(mlp[t],0.5) - pow(exp.mlp[t],0.5),2) mlp_fit_pred[t] = pow(pow(mlp_pred[t],0.5) - pow(exp.mlp[t],0.5),2) } mlp_fit_stat = sum(mlp_fit[2:7],mlp_fit[9:(nyears)]) mlp_fit_pred_stat = sum(mlp_fit_pred[2:7],mlp_fit_pred[9:(nyears)]) #... May counts 2010 onward (assumes lognormal residual error) ############################# for (t in 19:(nyears+1)) { mc[t] ~ dlnorm(log(Nm[t])-0.5/tau.m,tau.m) # Post predictive check exp.mc[t] = Nm[t] mc_pred[t] ~ dlnorm(log(Nm[t])-0.5/tau.m,tau.m) mc_fit[t] = pow(pow(mc[t],0.5) - pow(exp.mc[t],0.5),2) mc_fit_pred[t] = pow(pow(mc_pred[t],0.5) - pow(exp.mc[t],0.5),2) } mc_fit_stat = sum(mc_fit[19:(nyears+1)]) mc_fit_pred_stat = sum(mc_fit_pred[19:(nyears+1)]) #... November counts 1992 onward (assumes lognormal residual error) ####################### for (t in 1:nyears) { n[t] ~ dlnorm(log(Nn[t]*bias[t])-0.5/tau.n,tau.n) # Post predictive check exp.n[t] = Nn[t]*bias[t] n_pred[t] ~ dlnorm(log(Nn[t]*bias[t])-0.5/tau.n,tau.n) n_fit[t] = pow(pow(n[t],0.5) - pow(exp.n[t],0.5),2) n_fit_pred[t] = pow(pow(n_pred[t],0.5) - pow(exp.n[t],0.5),2) } n_fit_stat = sum(n_fit[]) n_fit_pred_stat = sum(n_fit_pred[]) # Likelihood for pre-November harvest in Denmark ########################################## for (t in 1:nyears) { alpha[t] <- min(0.99999,max(0.00001,ho[t]/hd[t])) w[t] ~ dbin(alpha[t],W[t]) # Post predictive check exp.w[t] = alpha[t]*W[t] w_pred[t] ~ dbin(alpha[t],W[t]) w_fit[t] = pow(pow(w[t],0.5) - pow(exp.w[t],0.5),2) w_fit_pred[t] = pow(pow(w_pred[t],0.5) - pow(exp.w[t],0.5),2) } w_fit_stat = sum(w_fit[]) w_fit_pred_stat = sum(w_fit_pred[]) # Likelihood for proportion of young in November ########################################## for (t in 1:nyears) { logit.gamma[t] <- beta0 + beta1*days[t] gamma[t] <- exp(logit.gamma[t])/(1+exp(logit.gamma[t])) r[t] <- gamma[t]/(1-gamma[t]) psi[t] <- r[t]*(1-dv*hn[t]-dv*ho[t]) / (1-(hn[t]+ho[t]) + r[t]*(1-dv*hn[t]-dv*ho[t])) # } prob[t] ~ dbeta(size[t]*psi[t],size[t]*(1-psi[t])) p[t] ~ dbin(prob[t],P[t]) # Post predictive check exp.p[t] = prob[t]*P[t] # expectation for p p_pred[t] ~ dbin(prob[t],P[t]) p_fit[t] = pow(pow(p[t],0.5) - pow(exp.p[t],0.5),2) p_fit_pred[t] = pow(pow(p_pred[t],0.5) - pow(exp.p[t],0.5),2) } p_fit_stat = sum(p_fit[]) p_fit_pred_stat = sum(p_fit_pred[]) # Likelihood of harvests ################################################################# for (t in 1:nyears) { #... Norway harvn[t] <- 10*(max(0.00001,hn[t]*Nm[t]*pow(theta[t],4/12)*(1+dv*r[t]))) hnor_lam[t] = harvn[t] #*(1-cr*rate[t]) hnor[t] ~ dpois(hnor_lam[t]) # Post predictive check hnor_pred[t] ~ dpois(hnor_lam[t]) hnor_fit[t] = pow(pow(hnor[t],0.5) - pow(hnor_lam[t],0.5),2) hnor_fit_pred[t] = pow(pow(hnor_pred[t],0.5) - pow(hnor_lam[t],0.5),2) #... Denmark harvd[t] <- 10*(max(0.00001,hd[t]*Nm[t]*pow(theta[t],4/12)*(1+dv*r[t]))) hden_lam[t] = harvd[t] hden[t] ~ dpois(hden_lam[t]) # Post predictive check hden_pred[t] ~ dpois(hden_lam[t]) hden_fit[t] = pow(pow(hden[t],0.5) - pow(hden_lam[t],0.5),2) hden_fit_pred[t] = pow(pow(hden_pred[t],0.5) - pow(hden_lam[t],0.5),2) } hnor_fit_stat = sum(hnor_fit[]) hnor_fit_pred_stat = sum(hnor_fit_pred[]) hden_fit_stat = sum(hden_fit[]) hden_fit_pred_stat = sum(hden_fit_pred[]) # Derived parms -------------------------------------------------------------------------------------- for (t in 1:nyears) { H[t] <- Nm[t]*pow(theta[t],4/12)*((hn[t]+hd[t]) + r[t]*(dv*hn[t]+dv*hd[t])) h[t] <- hn[t]+hd[t] } } # End model