data { int nNeighbors; int ncells; int nzones; int ntime; int y[ncells,ntime]; vector[ntime] time; int nrep[ncells,ntime]; int zone[ncells]; int m[ncells,ntime]; real search[ncells,ntime]; int known[ncells,ntime]; //int photos[ntime]; int hpd[ncells]; int air[ncells,ntime]; real distSq[ncells,nNeighbors]; real muFPC[ncells]; real muIMP[ncells]; //int adj[ncells,nNeighbors]; } parameters { real ptdet; real pknown; real baiteff; real pdelin; real alpha1; real delta1; real beta0; real beta1; real beta2; real beta3; real sigma; real rho0; real muN[ncells,ntime]; real N[ncells,ntime]; real alpha0[nzones]; //real alpha01[nzones]; //real alpha02[nzones]; real alpha2[nzones]; real delta0[nzones]; real muhpd[nzones]; real lambdaN[nzones,ntime]; real theta[nzones,ntime]; real psi[nzones,ntime]; real scN[nzones]; } transformed parameters{ real Np[ncells,ntime]; matrix[nzones,ntime] murep; matrix[nzones,ntime] lambda; real efflambda[ncells,ntime]; real Ez[ncells,ntime]; vector[ncells] Hgamma; matrix[ncells,ntime] gamma; matrix[ncells,nNeighbors] gammaDistPairs[ntime];//b is an array; maybe need to move this to the model? matrix[ncells,ntime] lkill; matrix[ncells,ntime] phi; matrix[ncells,ntime] lfind; for (t in 1:ntime){ for (j in 1:ncells){ lfind[j,t]=poisson_lcdf(m[j,t]|N[j,t]); lkill[j,t]=1-((1-(lfind[j,t] * search[j,t])) * ((1-(baiteff * pdelin * search[j,t])))); phi[j,t]=1-((1-lkill[j,t])*(1-(baiteff*air[j,t]))); } } for (t in 1:ntime){ for (l in 1:nzones){ //murep[l,t]=exp(alpha0[l] + alpha01[l] * (1-photos[t]) + alpha02[l] * photos[t] + alpha1*muhpd[l] + alpha2[l] * time[t]); murep[l,t]=exp(alpha0[l] + alpha1*muhpd[l] + alpha2[l] * time[t]); lambda[l,t]=exp(delta0[l] + delta1 * time[t]); } for (j in 1:ncells){ Np[j,t]=N[j,t]*ptdet; efflambda[j,t] = murep[zone[j],t]*lambda[zone[j],t]; } } for (t in 1:ntime){ for (c in 1:nNeighbors){ for (j in 1:ncells){ if (t<2){ gammaDistPairs[t,j,c]=0; } else { //gammaDistPairs[t,j,c]=(1 - 1/(1+(distSq[j,c] * adj[j,c]/sigma^2))) * Ez[c,t-1]; gammaDistPairs[t,j,c]=(1 - 1/(1+(distSq[j,c]/sigma^2))) * Ez[c,t-1]; } } } } for (t in 1:ntime){ for (j in 1:ncells){ Hgamma[j]=inv_logit(beta0 + beta1*muFPC[j] + beta2*muIMP[j] + beta3*muFPC[j]*muIMP[j]); if (t < 2){ Ez[j,t] = Hgamma[j]; } else { gamma[j,t] = (1 - ((prod(gammaDistPairs[t,j,1:nNeighbors])) * (1-rho0)))*Hgamma[j]; Ez[j,t] = Ez[j,t-1] * (1-(phi[j,t-1] * (1-gamma[j,t]))) + (1 - Ez[j,t-1]) * gamma[j,t];//gamma is p(col); phi is p(extinction) } } } } model { // priors alpha1~normal(0,100); sigma~gamma(0.6,0.75); rho0~beta(1,99); ptdet~beta(6.5,3.5); pknown~uniform(0,1); baiteff~beta(16,4); delta1~normal(0.5,10); beta0~normal(0,100); beta1~normal(0,100); beta2~normal(0,100); beta3~normal(0,100); pdelin~beta(6.5,3.5); alpha2~normal(0,100); for (t in 1:ntime){ for (l in 1:nzones){ lambdaN[l,t]~gamma(0.001,0.001); psi[l,t]~uniform(0,1); theta[l,t]~uniform(0,1);////prior for 16cellzone-level prevalence alpha0[l]~normal(0,100); //alpha01[l]~normal(0,100); //alpha02[l]~normal(0,100); alpha2[l]~normal(0,100); delta0[l]~normal(0,100); muhpd[l]~gamma(0.001,0.001); muN[l,t]~normal(5,1.5); scN[l]~uniform(0,0.1); } } for (t in 1:ntime){ for (j in 1:ncells){ N[j,t]~lognormal(muN[zone[j],t],scN[zone[j]]); target += poisson_lpmf(hpd[j]|muhpd[zone[j]]); if (m[j,t] == 0){ target += log(theta[zone[j],t]); } else { target += log1m(theta[zone[j],t]) + poisson_lpmf(m[j,t] | Np[j,t]) - poisson_lccdf(0 | Np[j,t]); } } } for (t in 1:ntime){ for (j in 1:ncells){ if (nrep[j,t]==0){ target += log_sum_exp(bernoulli_lpmf(1 | psi[zone[j],t]),bernoulli_lpmf(0 | psi[zone[j],t]) + poisson_lpmf(nrep[j,t] | murep[zone[j],t])); } else { target += bernoulli_lpmf(0 | psi[zone[j],t]) + poisson_lpmf(nrep[j,t] | murep[zone[j],t]); } if (y[j,t]==0){ target += log_sum_exp(bernoulli_lpmf(1 | Ez[j,t]),bernoulli_lpmf(0 | Ez[j,t]) + poisson_lpmf(y[j,t] | efflambda[j,t])); } else { target += bernoulli_lpmf(0 | Ez[j,t]) + poisson_lpmf(y[j,t] | efflambda[j,t]); } if (known[j,t]==0){ target += (log_sum_exp(log(Ez[j,t]) + bernoulli_lpmf(known[j,t]|pknown),log(1-Ez[j,t]))); } else { target += (log(Ez[j,t]) + bernoulli_lpmf(known[j,t]|pknown)); } } } }