#SIMULATION STUDY:

#1) Fix \beta_0 and \beta_1
set.seed(1) #Seed for randomness.
beta_0 <- round(runif(n = 1, min = -2, max = 2), 1)    #-0.9
beta_1 <- round(runif(n = 1, min = 0, max = 0.05), 2)  #0.02


#2) Create p_1 = \frac{e^{\beta_0 + \beta_1 x}}{1 + e^{\beta_0 + \beta_1 x}}
n <- 10000 #Number of observations.
X <- runif(n = n, min = 0, max = 10)
X <- floor(X) 
p_1 <- exp(beta_0 + beta_1 * X) / (1 + exp(beta_0 + beta_1 * X))
range(p_1)
#[1] 0.2890505 0.3273930


#3) Simulate Y_1 \sim Bernouilli(p_1)
Y_1 <- rbinom(n = n, size = 1, prob = p_1)


#4) Fit the glm (check you get \beta_0, \beta_1 and p_1 correct)
glm_lev1 <- glm(Y_1 ~ X, family = binomial) 
coef(glm_lev1)
#(Intercept)           X 
#-0.89642433  0.01796293
#They are very similar values to beta_0 and beta_1, as expected:
c(beta_0, beta_1)
#-0.9         0.02

glm_probs = predict(glm_lev1, type = "response")
#glm_probs are very similar to p_1, as expected, see for example:
glm_probs[1:10]
p_1[1:10]


#Rows of each top category:
aux0 <- X[which(Y_1 == 0)]
aux1 <- X[which(Y_1 == 1)]


#5) Fix \sigma, \beta_{00}, \beta_2 and the random factors
beta_00 <- round(runif(n = 1, min = -3, max = 3), 1)  #-2.1
beta_2 <- round(runif(n = 1, min = 0, max = 2), 1)    #0.6
sigma <- round(runif(n = 1, min = 0, max = 1), 1)     #0.2

r0  <- rnorm(length(aux0), 0, sigma) #random factor
r1  <- rnorm(length(aux1), 0, sigma) #random factor


#6) Create p_{2ij} = \frac{e^{\beta_{00} + \beta_2 x}}{1 + e^{\beta_{00} + \beta_2 x}}  
#and draw Y_{2ij} \sim Bernouilli(p_{2ij}) for j=0,1
p <- function(x, a, b) exp(a + b * x) / (1 + exp(a + b * x))
Y_2i0 <- rbinom(length(aux0), 1, prob = p(aux0, beta_00 + r0, beta_2))
Y_2i1 <- rbinom(length(aux1), 1, prob = p(aux1, beta_00 + r1, beta_2))


Y_2 <- c(Y_2i0, Y_2i1)
Y_1_aux <- c(Y_1[which(Y_1 == 0)], Y_1[which(Y_1 == 1)])
data <- data.frame(X, Y_1_aux, Y_2)
names(data) <- c("X", "Y1", "Y2")
data[c(1:6,9995:10000),]
#X Y1 Y2
#5  0  1
#9  0  0
#2  0  1
#8  0  0
#9  0  0
#6  0  0
#4  1  0
#5  1  1
#5  1  1
#6  1  1
#0  1  1
#6  1  0


#7) Fit the glmer (check you get \beta_00, \beta_2, \sigma correct)
library(lme4)
glm_lev2 <- glmer(Y2 ~ X + (1|Y1), data = data, family=binomial)
coef(glm_lev2)
#$Y1
#    (Intercept)           X
#0   0.3872395  -0.003517575
#1   0.3872395  -0.003517575
#These values aren't similar to the generated beta_00 and beta_2.
VarCorr(glm_lev2)
#Groups Name        Std.Dev.  
#Y1     (Intercept) 1.9464e-07
#This value isn't similar to the generated sigma.