`panel10pt2.fn` <- function(ni=3000,nb=2000,nthin=1,nchain=3,nz=160){ # # Jolly-Seber Model w/Individual Heterogeneity -- This R script fits the # occupancy formulation of the Jolly-Seber model induced by data # augmentation to the dipper data, and allows for individual # heterogeneity in survival probability. This script produced the # results in Table 10.3. The Markov chains mix slowly, so a # long run is necessary to reduce Monte Carlo error. Note also that # results in Table 10.3 were run with nz=100 which probably was not high enough 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) } sigma~dunif(0,10) tau<-1/(sigma*sigma) for(j in 1:M){ eta[j]~dnorm(0,tau) for(i in 1:(T-1)){ logit(phi[j,i])<-lphi[i] + eta[j] } } # one parameterization is to give each gamma a uniform(0,1) prior # the other way is to put the dirichlet prior on the multinomial equivalent 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:T){ b[i]<-(1.4-a[i]) gamma[i]~dbeta(a[i],b[i]) lgamma[i]<- log(gamma[i]/(1-gamma[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]~dbern(gamma[1]) mu[i]<-z[i,1]*ptmp[1] y[i,1]~dbern(mu[i]) 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]~dbern(muz[i,j]) muy[i,j]<-z[i,j]*ptmp[j] y[i,j]~dbern(muy[i,j]) } } 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,sigma=1) } 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) }