##### RUN LINES 4-27 TO GENERATE ONE SIMULATED DATA SET FOR EXERCISE 9 #####
##### (Based on Kery 2010) #####

n.groups <- 2
n.sample <- 100 
n <- n.groups* n.sample
x <- rep(1:n.groups, rep(n.sample, n.groups))
pop <- factor(x, labels = c("NOrth","South"))
wet1 <- sort(runif(n.sample,0,1))
wet2 <- sort(runif(n.sample,0,1))
habitat <- c(wet1,wet2)
N <- round(runif(n,10,50))
Xmat <- model.matrix(~pop+habitat)
Xmat.p <- model.matrix(~habitat)
beta.vec <- c(-2,4)
lin.pred <- Xmat.p[,]%*%beta.vec
exp.p <- exp(lin.pred)/(1+exp(lin.pred))
C <- rbinom(n=n, size=N, prob=exp.p)
P_bi <- rbinom(n=n, size=1, prob=exp.p)
beta.vec_p <- c(-2,2,3)
lin.pred_p <- Xmat[,]%*%beta.vec_p
lambda <- exp(lin.pred_p)
C_p <- rpois(n=n, lambda=lambda)
flowers <- P_bi*C_p
roses <- as.data.frame(cbind(habitat,pop,flowers))
roses$pop <- factor(roses$pop)
levels(roses$pop) <- c("North","South")
roses