# Script which generates a simulated dataset (note the use of "seed")

# 	R: number of surveyed sites
# 	T: number of surveys at each site
# 	lambda: expected abundance
# 	p0: intercept of detection probability (logit scale)
# 	p1: slope of detection probability on covariate x (logit scale)
# 	pSD: standard deviation of normal distribution from 
#     which observer effects u are drawn
# 	nG: number of different observers (called K in section 2)
# 	which covariate values x are generated

R <- 100
T <- 10

lambda <- 5 #2
p0 <- 1 #-1
p1 <- 1

pSD <- 0.5
nG <- 40

gID <- rep(1:nG, each=(R*T)/nG)
gID <- matrix(gID, R, T, byrow=TRUE) # matrix of observer IDs

set.seed(11)

x <- rnorm(R) # covariate of p
u0 <- rnorm(nG, 0, pSD) # random effects vector
u <- rep(u0, each=(R*T)/nG)
u <- matrix(u, R, T, byrow=TRUE)
p <- plogis(p0 + p1*x + u) # detection probs
dim(p)

N <- rpois(R, lambda)
y <- matrix(NA, R, T) # matrix of counts

y[] <- rbinom(R*T, N, p)
table(y)