n.groups <- 20
n.sample <-  3
pH<- runif(20,3,7)
n <- n.groups* n.sample
x <- rep(1:n.groups, rep(n.sample, n.groups) )
ph.data <- rep(0,n)
for (i in 1:n) {
  ph.data[i] <- pH[x[i]] }
occ <- c(rep(0,30),rep(1,30))

error.sd <- rnorm(n,0,0.61)
pop <- factor(x, labels = c(1:20))
Xmat <- model.matrix(~ ph.data + occ)
#print(Xmat,dig=2)

X <- list(c=1:20)
beta1.vec <- rapply(X,function(x) {rnorm(20,-8.6,1.05)},how ="replace")
beta2.vec <- rapply(X,function(x) {rnorm(20,1.76,0.18)},how ="replace")
beta3.vec <- rapply(X,function(x) {rnorm(20,0,0.20)},how ="replace")
beta.vec <- rbind(beta1.vec$c,beta2.vec$c,beta3.vec$c)
alfa <- rapply(X,function(x) {rnorm(20,0,0.37)},how ="replace")
r.eff <- NULL
for (i in 1:20){
  r.eff <- c(r.eff,rep(alfa$c[i],3))
}
lin.pred <- Xmat[,]%*%beta.vec[,sample(1:20,1)]
lin.pred <- as.data.frame(lin.pred) + r.eff

exp.p <- exp(lin.pred)/(1+exp(lin.pred))
alive <- rep(0,length(exp.p$V1))
for (i in 1:length(exp.p$V1)){
  alive[i] <- rbinom(1,3,exp.p$V1[i])
}

sel <- seq(1,60,3)
snail = rep(0,60)
for (i in 1:20) {
  if (alive[sel[i]]== 3) { snail[sel[i]:(sel[i]+2)] <- c(1,1,1)}
  if (alive[sel[i]]== 2) { snail[sel[i]:(sel[i]+2)] <- c(0,1,1)}
  if (alive[sel[i]]== 1) { snail[sel[i]:(sel[i]+2)] <- c(0,0,1)}
   }
snail


occupied <- as.data.frame(cbind(snail,ph.data,occ,pop))
occupied
library(nlme)
library(bbmle)
library(lme4)


occupied$occ <- as.factor(occupied$occ)


####################################

#######  Data for students  ########

occupied