Distance sampling in BUGS

[camundson]

I was hoping to sneak in a winBUGS question as it relates to unmarked
- sorry if it's inappropriate.

I ran analyses in unmarked and am looking to compare estimates in
winBUGS. I'm wondering if it's possible to include covariates on the
detection function when attempting to estimate abundance using
distance sampling (circular point counts) and the N mixture model in
BUGS. I'm using code provided by Andy at the recent Patuxent
workshop, but can't quite wrap my head around where I would apply
covariates on detection (f[i]?, pcap?). I'm using simulated data
tweaked from what was provided. Does anyone have experience with this
type of problem?

Also, can anyone help me define where the "dclass" and "nind" portions
of code work into estimates of detection or abundance? I tried
running it without and it appears to affect beta0 and the deviance -
but I'm also not running enough simulations to get convergence so it
could be my imagination.

Thanks!

Courtney

Modified from code provided by A. Royle
#######################################################
set.seed(1234)
nsites<- 840
x<-rnorm(nsites) # a covariate value
lambda<- exp(2 + .5*x) # density per "square" poisson -
#providing relationship with abundance
N<-rpois(nsites,lambda) # site-specific abundances
sigma<- 1.0 # distance function parameter - shape
parameter
B<- 5.6 # 2 x max count distance, simulating on a
square, scaled = d/100 m

data<-NULL
for(i in 1:nsites){

if(N[i]==0){
data<-rbind(data,c(i,NA,NA,NA,NA))
next
}
### Simulation data on a square
u<-runif(N[i],0,B) # generating individuals uniformly
v<-runif(N[i],0,B)
d<- sqrt( (u - B/2)^2 + (v-B/2)^2)
### But we can only count individuals on a circle so we truncate p
here
p<- ifelse(d< (B/2),1,0)*exp(-d*d/(2*(sigma^2)))
y<-rbinom(N[i],1,p)

u<-u[y==1] # this subsets to "captured" individuals only
v<-v[y==1]
d<-d[y==1]
y<-y[y==1]

if(sum(y)>0)
data<-rbind(data,cbind(rep(i,sum(y)),y,u,v,d))
else
data<-rbind(data,c(i,NA,NA,NA,NA))
}

## some examination of the dataset we created
summary(data)
data[1:10,]
plot(data[,c("u","v")])


ncap<- table(data[,1]) # equals "1" if no individuals were captured
sites0<-data[is.na(data[,2]),][,1]
ncap[as.character(sites0)]<-0
y<-as.vector(ncap)

maxd<-B/2
delta<- .1 # bin size, very small to facilitate integration
xg<-seq(delta/2,maxd,delta) #generates midpoints up to max distance
dclass<-data[,5]%/%delta + 1 # decides what distance class the
observation goes in,
# starts at 1 so it's not zero
midpt<-xg # renames the midpoints
prob<- 2*midpt*delta/(maxd*maxd)
nG<-length(xg)
dclass<-dclass[!is.na(data[,2])]
nind<-length(dclass)

sink("model1.txt")
cat("
model{
sigma~dunif(0,10)
beta0 ~ dunif(-10,10)
beta1 ~ dunif(-10,10)
for(i in 1:nG){

# for the following Pr(detection|x)Pr(x) =
# exp(-(x^2/2*sigma^2)f(x))) - here, p[i]*pi[i]
log(p[i])<- -xg[i]*xg[i]/(2*sigma*sigma) # first part of
trapezoidal #approximation to integrate over delta
pi[i]<- ( 2*xg[i]*delta )/(maxd*maxd) # f(x) - second part of
#trap. approximation - area,
f[i]<- p[i]*pi[i] #account for the spatial coverage bias?
(Chelgren et al. 2011)
fc[i]<- f[i]/pcap # fi/sum(fi)
}
pcap<-sum(f[]) # overall detection probability

#How are the following 4 lines linked to the model?

for(i in 1:nind){
dclass[i] ~ dcat(fc[]) #single binomial trial with categorical
distribution - here
## same as multinomial distribution since it's a single trial
}

for(i in 1:nsites){

# binomial model for # of captured individuals
y[i]~ dbin(pcap,N[i]) #number of ind. detected within effective
distance at site i

## abundance model
N[i]~dpois(lambda[i])
log(lambda[i])<- beta0 + beta1*x[i]
}




}
",fill=TRUE)
sink()

Nst<- y + 1
ni<-2000; nb<-1000; nthin<-2; nc<-3
bugs.data <- list
("y","xg","delta","maxd","nind","nsites","nG","dclass","x")
inits <- function(){ list
(sigma=runif(1,1,10),beta0=0,beta1=0,N=Nst) }
parameters <- c("sigma","beta0","beta1")
out <- bugs (bugs.data, inits, parameters, "model1.txt",
n.thin=nthin,n.chains=nc, n.burnin=nb,n.iter=ni,debug=TRUE)

[Jeffrey Royle]

hi Courtney,
 You can put covariates in the sigma parameter -- but you have to then define it for each point or individual depending on the nature of the covariate. Is the covatiate a point level or individual level covariate?  If you send me your data and R script I can try to hack that for you.
In the R/BUGS script from the workshop, the encounter probability loop would need to be nested within a larger loop (over points or individuals depending on the nature of the covariate).

nind in the model specification is the total number of _captured_ individuals.

dclass is the "distance class interval" that each observation is in -- this is an interval from 1 to however many distance classes there are.

 Sorry that the script is not very clear -- this was very late going into the workshop materials.

 Anyhow, email me specific questions and I can help get  some covariates included in the model.
regards
andy