salamanderRemovals.Bayes = function() {

# read data from file
  df = read.csv('salamanderRemovals2001.csv')
  y = as.vector(df[,'Desmognathus'])
  #  y = as.vector(df[,'Eurycea'])
  site = as.vector(df[, 'site'])
  pass = as.vector(df[, 'pass'])


# model specification in WinBUGS
modelFilename = 'salamanderBayes.txt'
cat('
model {

# priors
pMean ~ dunif(0,1)
etaMean <- log(pMean) - log(1-pMean)
etaSigma ~ dunif(0,10)
etaTau <- pow(etaSigma, -2)

beta ~ dnorm(0., 0.01)
lambdaMean <- exp(beta)

sigma ~ dunif(0,10)
tau <- pow(sigma, -2)


# models of heterogeneity
for (i in 1:nsites) {
    eta[i] ~ dnorm(etaMean, etaTau)
    logit(p[i]) <- eta[i]
    b[i] ~ dnorm(beta, tau)
}

# marginal model of counts
for (i in 1:n) {
    log(lambda[i]) <- b[site[i]]
    piProb[i] <- p[site[i]] * pow(1-p[site[i]], pass[i]-1)
    mu[i] <- lambda[i] * piProb[i]
    y[i] ~ dpois(mu[i])
}

}
', fill=TRUE, file=modelFilename)



# arguments for bugs()
  nsites = max(site)
  data = list(n=length(y), nsites=nsites, y=y, site=site, pass=pass)
  
  lambdaGuess = mean(y)


  params = list('pMean', 'lambdaMean')
  inits = function() {
    list(pMean=runif(1,0,1), beta=rnorm(1,mean=log(lambdaGuess),sd=1), sigma=runif(1,0,2), etaSigma=runif(1,0,2) )
  }


# call to bugs()
  library(R2WinBUGS)
  fit = bugs(data, inits, params, model.file=modelFilename,
    n.chains=1, n.iter=300000, n.burnin=100000, n.thin=20,  bugs.seed=sample(1:9999,size=1), debug=FALSE, DIC=FALSE)

  fit
}