`panel10pt1.fn` <- function(ni=3000,nb=2000,nthin=1,nchain=3,nz=160){ # Jolly-Seber model -- occupancy formulation induced by data augmentation. # R script for fitting model to the dipper data. Requires the file # dipper.data -- load it into your R workspace using # source("dipper.data.R"). Results are consistent with "prior 1" results # in Table 10.1. There is a slight difference because M was set too low # for the results reported in Table 10.1 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] } } # cprob = multinomial entrance probs, gammaSA = Schwarz and Arnason entrance # probs. gamma = colonization probabilities for(i in 1:T){ gamma[i]~dunif(0,1) 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]<-(1-z[j,i-1])*(z[j,i]) } } # time specific estimates of population size and births 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) }