`panel10pt1b.fn` <- function(ni=3000,nb=2000,nthin=1,nchain=3,nz=160){ #Jolly-Seber Model -- occupancy formulation induced by data augmentation -- # using the alternative prior specification (see Section 10.3.7). # Produces results consistent with those in Table 10.1 ("prior 2"). library("R2WinBUGS") y<-dipper.data\$y # set pre-capture NA to 0 -- this is important y[is.na(y)]<-0 y<-rbind(y,matrix(rep(0,ncol(y)*nz),ncol=ncol(y),byrow=TRUE)) M<-nrow(y) T<-ncol(y) z<-y sink("jsmodel.txt") cat(" model { for(j in 1:(T-1)){ phi0[j]~dunif(0,1) lphi[j]<-log(phi0[j]/(1-phi0[j])) } ptmp[1]<- 1 ptmp[T]<- 1 for(i in 2:(T-1)){ ptmp[i]~dunif(0,1) } for(j in 1:M){ for(i in 1:(T-1)){ logit(phi[j,i])<-lphi[i] } } a[1]<- 1.4/14 a[2]<- 1.4/13 a[3]<- 1.4/12 a[4]<- 1.4/11 a[5]<- 1.4/10 a[6]<- 1.4/9 a[7]<- 1.4/8 for(i in 1:7){ b[i]<-(1.4-a[i]) } for(i in 1:T){ gamma[i]~dbeta(a[i],b[i]) gammaSA[i]<- cprob[i]/psi.imp } cprob[1]<- gamma[1] cprob[2]<- gamma[2]*(1-gamma[1]) cprob[3]<- gamma[3]*(1-gamma[2])*(1-gamma[1]) cprob[4]<- gamma[4]*(1-gamma[3])*(1-gamma[2])*(1-gamma[1]) cprob[5]<- gamma[5]*(1-gamma[4])*(1-gamma[3])*(1-gamma[2])*(1-gamma[1]) cprob[6]<- gamma[6]*(1-gamma[5])*(1-gamma[4])*(1-gamma[3])*(1-gamma[2])*(1-gamma[1]) cprob[7]<- gamma[7]*(1- gamma[6])*(1-gamma[5])*(1-gamma[4])*(1-gamma[3])*(1-gamma[2])*(1-gamma[1]) psi.imp<- sum(cprob[1:7]) for(i in 1:M){ z[i,1]~dbin(gamma[1],1) mu[i]<-z[i,1]*ptmp[1] y[i,1]~dbin(mu[i],1) recruitable[i,1]<-1 for(j in 2:T){ survived[i,j]<- phi[i,j-1]*z[i,j-1] recruitable[i,j]<- recruitable[i,(j-1)]*(1-z[i,j-1]) muz[i,j]<- survived[i,j] + gamma[j]*recruitable[i,j] z[i,j]~dbin(muz[i,j],1) muy[i,j]<-z[i,j]*ptmp[j] y[i,j]~dbin(muy[i,j],1) } } for(j in 1:M){ recruit[j,1]<-z[j,1] for(i in 2:T){ recruit[j,i]<-z[j,i-1]*(1-z[j,i]) } } for(i in 1:T){ N[i]<-sum(z[1:M,i]) B[i]<-sum(recruit[1:M,i]) } for(i in 1:M){ Nind[i]<-sum(z[i,1:T]) Nalive[i]<-1-equals(Nind[i],0) } Nsuper<-sum(Nalive[1:M]) } ",fill=TRUE) sink() data <- list ("y","M","T") inits <- function(){ list ( phi0=runif(T-1,.2,.9),z=z) } parameters <- c("ptmp","phi0","N","Nsuper","gamma","gammaSA","cprob") out <- bugs (data, inits, parameters, "jsmodel.txt", n.thin=nthin,n.chains=nchain, n.burnin=nb,n.iter=ni,debug=TRUE) return(out) }