# Simulate outcome, random numbers and groups
set.seed(283)
Out_1 <- rbinom(500, 1, 0.4)
Random <- rnorm(1:length(Out_1))
ID <- seq(1,length(Out_1), 1)
trial_db <- as.data.frame(cbind(ID, Out_1, Random))

# Make outcome logical so I can set 'paired' outcome to same or opposite
trial_db$Out_1log <- trial_db$Out_1 == 1

result <-  list()
i <- 0.05
rindex <- 0
results <- repeat {

  trial_db$Out_2 <- with(trial_db, ifelse(Random <= unname(quantile(trial_db$Random, i)), as.numeric(Out_1log), as.numeric(!Out_1log)))

  # Make in to long format Random effects

  one_db <- trial_db %>% dplyr::select(ID, Out_1) %>% rename(OUT = Out_1)
  two_db <- trial_db %>% dplyr::select(ID, Out_2) %>% rename(OUT = Out_2)
  long_db <- rbind(one_db, two_db) %>% arrange(ID)
  long_db$ID <- as.factor(as.character(long_db$ID))

  mod.adp <- GLMMadaptive::mixed_model(fixed = OUT ~ 1, random = ~ 1 | ID, data = long_db, family = binomial(link = "logit"))
  mod.tmb <- glmmTMB::glmmTMB(OUT ~ (1|ID), data = long_db, family = binomial(link = "logit"))
  mod.lme <- lme4::glmer(OUT ~ (1|ID), data = long_db, family = binomial(link = "logit"))
  ADAP <- round(icc(mod.adp)[[1]], 3)
  TMB <- round(icc(mod.tmb)[[1]], 3)
  LME4 <- round(icc(mod.lme)[[1]], 3)

  estimates <- cbind(i, ADAP, TMB, LME4)
  rindex <- rindex + 1
  result[[rindex]] <- estimates

  i <- i + 0.05

  if(i == 0.95) {
    break
  }
}

table <- do.call(rbind.data.frame, result)

i  ADAP   TMB  LME4
1  0.05 0.000    NA    NA
2  0.10 0.000    NA    NA
3  0.15 0.000    NA    NA
4  0.20 0.001    NA    NA
5  0.25 0.001    NA    NA
6  0.30 0.001    NA    NA
7  0.35 0.003    NA    NA
8  0.40 0.002    NA    NA
9  0.45 0.003    NA    NA
10 0.50 0.017    NA    NA
11 0.55 0.126 0.066 0.066
12 0.60 0.267 0.154 0.154
13 0.65 0.412 0.274 0.274
14 0.70 0.551 0.421 0.420
15 0.75 0.678 0.570 0.570
16 0.80 0.789 0.700 0.700
17 0.85 0.877 0.806 0.805
18 0.90 0.936 0.897 0.897
19 0.95 0.959 0.995 0.995