`panel6pt4.fn` <-
function(ni=3000,nb=1000,nthin=1,nc=2,nz=125){
# R script to fit closed population model with behavioral response and
# individual covariates using data augmentation in WinBUGS. Data are the
# microtus data with individual covariate body mass.  This model is 
# similar to Panel 6.4
library("R2WinBUGS")

y<-microtus.data[,1:5]
mass<-microtus.data[,6]
nind<-nrow(y)
J<-ncol(y)

first<-rep(NA,nind)
prevcap<-matrix(0,nrow=nind,ncol=J)
for(i in 1:nind){
first[i]<- min((1:J)[y[i,]==1])
if(first[i]<J)
  prevcap[i,(first[i]+1):J]<-1
}

# Augmenting the data, including covariates

y<-rbind(y,matrix(0,nrow=nz,ncol=J))
prevcap<-rbind(prevcap,matrix(0,nrow=nz,ncol=J))
mass<- (mass-mean(mass))/sqrt(var(mass))
mass<-c(mass,rep(NA,nz))

lagY<-matrix(NA,nrow=nrow(y),ncol=ncol(y))
lagY[,1]<-0
lagY[,2:J]<-y[,1:4]


sink("model.txt")
cat("
model {

mu~dunif(-10,10)
tau~dgamma(.01,.01)
sigma<-sqrt(1/tau)
psi~dunif(0,1)
a0~dunif(-10,10)
a1~dunif(-10,10)
a2~dunif(-10,10)
a3~dunif(-10,10)

for(i in 1:(nind+nz)){
 mass[i]~dnorm(mu,tau)
 z[i]~dbin(psi,1)
for(t in 1:J){
logit(p[i,t])<-a0*(1-prevcap[i,t]) + a1*prevcap[i,t]  + a2*lagY[i,t] + a3*mass[i]
muY[i,t]<-p[i,t]*z[i]
y[i,t] ~ dbern(muY[i,t])
}
}
 N<-sum(z[1:(nind+nz)])
}
",fill=TRUE)
sink()

data <- list ("y","nind","nz","prevcap","J","lagY","mass")
inits <- function(){
  list (a0=rnorm(1),a1=rnorm(1),a2=rnorm(1),a3=rnorm(1),z=rbinom(nind+nz,1,.5),mu=rnorm(1,0,.3),tau=runif(1,1,2))
}

parameters <- c("N","a0","a1","a2","a3","mu","sigma","psi")

out <- bugs (data, inits, parameters, "model.txt", n.thin=nthin,n.chains=nc, n.burnin=nb,n.iter=ni,debug=TRUE)


out


}