MultiSpeciesOcc = function() {

  Ymat = as.matrix(read.csv('detectionFreq.NH17.csv'))
  # Ymat = as.matrix(read.csv('detectionFreq.AL17.csv'))
  nrepls = 11

  start.time = Sys.time()
  
  # augment data matrix with an arbitrarily large number of zero row vectors
  nzeros = 151
  n = dim(Ymat)[1]
  nsites = dim(Ymat)[2]
  Yaug = rbind(Ymat, matrix(0, nrow=nzeros, ncol=nsites))

  # create arguments for bugs()
  data = list(n=n, nzeros=nzeros, R=nsites, J=nrepls, Y=Yaug)

  params = list('alpha', 'beta', 'rho', 'sigma.u', 'sigma.v', 'omega', 'N', 'phi','eta')

  inits = function() {
    omegaGuess = runif(1, n/(n+nzeros), 1)
    psi.meanGuess = runif(1, .25,1)
    p.meanGuess = runif(1, .25,1)
    rhoGuess = runif(1, 0,1)
    sigma.uGuess =  runif(1,0,1.5)
    sigma.vGuess =  runif(1,0,1.5)
    list(omega=omegaGuess, psi.mean=psi.meanGuess, p.mean=p.meanGuess,
         sigma.u=sigma.uGuess, sigma.v=sigma.vGuess, rho=rhoGuess,
         w=c(rep(1, n), rbinom(nzeros, size=1, prob=omegaGuess)),
         phi=rnorm(n+nzeros, log(psi.meanGuess/(1-psi.meanGuess)), sigma.uGuess),
         eta=rnorm(n+nzeros, log(p.meanGuess/(1-p.meanGuess)), sigma.vGuess),
         Z = matrix(rbinom((n+nzeros)*nsites, size=1, prob=psi.meanGuess), nrow=(n+nzeros))
         )
  }
  
  modelFilename = "MultiSpeciesOccModel.txt"

  cat("
model {

omega ~ dunif(0,1)

psi.mean ~ dunif(0,1)
beta <- log(psi.mean) - log(1-psi.mean)

p.mean ~ dunif(0,1)
alpha <- log(p.mean) - log(1-p.mean)

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

rho ~ dunif(-1,1)
var.eta <- tau.v/(1.-pow(rho,2))


for (i in 1:(n+nzeros)) {
    w[i] ~ dbern(omega)
    phi[i] ~ dnorm(beta, tau.u)I(-20,20)

    mu.eta[i] <- alpha + (rho*sigma.v/sigma.u)*(phi[i] - beta)
    eta[i] ~ dnorm(mu.eta[i], var.eta)I(-20,20)


    logit(psi[i]) <- phi[i]
    logit(p[i]) <- eta[i]

    mu.psi[i] <- psi[i]*w[i]
    for (k in 1:R) {
	Z[i,k] ~ dbern(mu.psi[i])
	mu.p[i,k] <- p[i]*Z[i,k]
	Y[i,k] ~ dbin(mu.p[i,k], J)
    }
}

n0 <- sum(w[(n+1):(n+nzeros)])
N <- n + n0
}
", fill=TRUE, file=modelFilename)
  

  # fit model to data using WinBUGS code
  library(R2WinBUGS)
  fit = bugs(data, inits, params, modelFilename,
    n.chains=5, n.iter=55000, n.burnin=5000, n.thin=50,  debug=FALSE, DIC=FALSE)
  
  end.time = Sys.time()
  elapsed.time = difftime(end.time, start.time, units='mins')
  cat(paste(paste('Posterior computed in ', elapsed.time, sep=''), ' minutes\n', sep=''))
  
  fit
}