library(MASS)
library(glmm)
library(GLMMadaptive)
library(lme4)

n <- 100 # Sample size
T <- 5 # Number of time points

beta <- c(-1,1,-0.2)  #True beta vector
Sig <- diag(1,T+1,T+1)  #Variance covariance matrix
Z <- matrix(c(1,1,0,0,0,0,1,0,1,0,0,0,1,0,0,1,0,0,1,0,0,0,1,0,1,0,0,0,0,1),ncol=6,byrow=TRUE)  #Z Design matrix

  Xit1 <- cbind(rbinom(n,1,0.5),rbinom(n,1,0.5),rbinom(n,1,0.5), rbinom(n,1,0.5),rbinom(n,1,0.5))  #Generate Bernoulli covariates
  Xit2 <- cbind(rnorm(n,0,1),rnorm(n,0,1),rnorm(n,0,1),rnorm(n,0,1), rnorm(n,0,1))  #Generate normal covariates
  Y <- matrix(rep(NA,n*T),ncol=T)  #Matrix to store simulated measurements. Each individual has 5 measurements in 1 row.
  for(row in 1:n){
    Ui <- mvrnorm(1,mu=rep(0,T+1),Sigma=Sig)  #Generate random intercepts U_i
    for(col in 1:T){
      z1 <- beta[1] + beta[2]*Xit1[row,col] + beta[3]*Xit2[row,col] + sum(Z[col,]*Ui)  #The RHS of model equation
      pi <- exp(z1)/(1+exp(z1))  #Predicted probability given covariates
      Y[row,col] <- runif(1) < pi  #Simulate binary response with predicted probability
    }
  }
  y <- as.vector(t(Y))  #Transforms Y matrix into vector
  xit1 <- as.vector(t(Xit1)) #Transforms Xit1 matrix into vector
  xit2 <- as.vector(t(Xit2)) #Transforms Xit2 matrix into vector
  t1 <- rep(c(1,0,0,0,0),n)  #Binary variable, = 1 if measurement was at time 1
  t2 <- rep(c(0,1,0,0,0),n)  #Binary variable, = 1 if measurement was at time 2
  t3 <- rep(c(0,0,1,0,0),n)  #Binary variable, = 1 if measurement was at time 3
  t4 <- rep(c(0,0,0,1,0),n)  #Binary variable, = 1 if measurement was at time 4
  t5 <- rep(c(0,0,0,0,1),n)  #Binary variable, = 1 if measurement was at time 5
  t <- rep(c(1,2,3,4,5),n)   #Stores time in non-binary form
  i <- rep(1:n,each=5)       #Subject label
  dat <- data.frame(y=y,xt1=xit1,xt2=xit2,t1=t1,t2=t2,t3=t3,t4=t4,t5=t5,i=factor(i),t=factor(t))
  model <- glmm(y ~ xt1 + xt2, random=list(y~0+i,~0+t1,~0+t2,~0+t3,~0+t4,~0+t5),varcomps.names=c("u","t1","t2","t3","t4","t5"),
                data=dat,family.glm=bernoulli.glmm)