###     A programme to obtain the power of parameters in 2 level
#       balanced model  with Binary response
#                    generated on 28/07/23
###~~~~~~~~~~~~~~~~~    Required packages  ~~~~~~~~~~~~~~~~~~~~~###
library(MASS)
library(lme4)
###~~~~~~~~~~~~~~~~~~~     Initial inputs    ~~~~~~~~~~~~~~~~~~~~###

set.seed(1)
siglevel <- 0.025
z1score <- abs(qnorm(siglevel))
simus <- 1000
n1low <- 10
n1high <- 10
n1step <- 1
n2low <- 10
n2high <- 50
n2step <- 5
npred <- 0
randsize <- 1
beta <- c(-0.405500)
betasize <- length(beta)
effectbeta <- abs(beta)
randcolumn <- 0
sigma2u <- matrix(c(0.250000), randsize, randsize)
n1range <- seq(n1low, n1high, n1step)
n2range <-seq(n2low, n2high, n2step)
n1size <- length(n1range)
n2size <- length(n2range)
totalsize <- n1size*n2size
finaloutput <- matrix(0, totalsize, 6*betasize)
rowcount <- 1
##-----------------        Inputs for model fitting       -----------------##

fixname <- "x0"
fixform <- "1"
randform <- "(1|l2id)"
expression <- paste(c(fixform, randform), collapse="+")
modelformula <- formula(paste("y ~", expression))
data <- vector("list", 2+length(fixname))
names(data) <- c("l2id", "y", fixname)

#####--------- Initial input for power in two approaches ----------------#####

powaprox <- vector("list", betasize)
names(powaprox) <- "b0"
powsde <- powaprox

cat("     The programme was executed at", date(),"\n")
cat("--------------------------------------------------------------------\n")

for (n2 in seq(n2low, n2high, n2step)) {
  for (n1 in seq(n1low, n1high, n1step)) {

    length <- n1*n2
    x <- matrix(1, length, betasize)
    z <- matrix(1, length, randsize)
    l2id <- rep(c(1:n2), each=n1)
    sdepower <- matrix(0, betasize, simus)
    powaprox[1:betasize] <- rep(0,betasize)
    powsde <- powaprox

    cat(" Start of simulation for sample sizes of ", n1, " micro and ", n2, "macro units\n")
    for (iter in 1:simus) {

      if (iter/10 == floor(iter/10)) {
        cat(" Iteration remain=", simus-iter,"\n")
      }
      ##--------------------------------------------------------------## 
      u <- mvrnorm(n2, rep(0, randsize), sigma2u)
      fixpart <- x * beta
      randpart <- rowSums(z*u[l2id, ])
      binomprob <- exp(fixpart+randpart)/(1+exp(fixpart+randpart))
      y <- rbinom(length, 1, binomprob)
      ##-------------------        Inputs for model fitting       ---------------##

      data$l2id <- l2id
      data$y <- y
      data$x0 <- x
      ###~~~~~~~~~~      Fitting the model using lmer funtion    ~~~~~~~~~~###

      (fitmodel <- glmer(modelformula, data, family=binomial(link="logit"), nAGQ="0"))

      ######~~~~~~~~~~   To obtain the power of parameter(s) ~~~~~~~~~~######

      estbeta <- fixef(fitmodel)
      sgnbeta <- sign(estbeta)
      sdebeta <- sqrt(diag(vcov(fitmodel)))
      for (l in 1:betasize)
      {
        cibeta <- estbeta[l]-sgnbeta[l]*z1score*sdebeta[l]
        if (estbeta[l]*cibeta > 0) powaprox[[l]] <- powaprox[[l]]+1
        sdepower[l, iter] <- sdebeta[l]
      }

      ##-------------------------------------------------------------------------##
    } ##  iteration end here

    ###---------                  Powers and their CIs             ---------###

    for (l in 1:betasize) {

      meanaprox <- powaprox[[l]]/simus
      Laprox <- meanaprox-z1score*sqrt(meanaprox*(1-meanaprox)/simus)
      Uaprox <- meanaprox+z1score*sqrt(meanaprox*(1-meanaprox)/simus)
      meansde <- mean(sdepower[l,])
      varsde <- var(sdepower[l,])
      USDE <- meansde-z1score*sqrt(varsde/simus)
      LSDE <- meansde+z1score*sqrt(varsde/simus)
      powLSDE <- pnorm(effectbeta[l]/LSDE-z1score)
      powUSDE <- pnorm(effectbeta[l]/USDE-z1score)
      powsde[[l]] <- pnorm(effectbeta[l]/meansde-z1score)


      ###---------   Restrict the CIs within 0 and 1   ---------##
      if (Laprox < 0) Laprox <- 0
      if (Uaprox > 1) Uaprox <- 1
      if (powLSDE < 0) powLSDE <- 0
      if (powUSDE > 1) powUSDE <- 1

      finaloutput[rowcount, (6*l-5):(6*l-3)] <- c(Laprox, meanaprox, Uaprox)
      finaloutput[rowcount, (6*l-2):(6*l)] <- c(powLSDE, powsde[[l]], powUSDE)

    } 

    ###~~~~~~~~~~    Set out the results in a data frame    ~~~~~~~~~~###

    rowcount <- rowcount+1
    cat("--------------------------------------------------------------------\n")
  } ## end of the loop  over the first level
} ## end of the loop  over the second level

###---------         Export output in a file                      ---------###

finaloutput <- as.data.frame(round(finaloutput, 3))
output <- data.frame(cbind(rep(n2range, each=n1size), rep(n1range, n2size), finaloutput))
names(output) <- c("N", "n", "zLb0", "zpb0", "zUb0", "sLb0", "spb0", "sUb0")
write.table(output, "powerout.txt", sep="\t ", quote=FALSE, eol="\n", dec=".", col.names=TRUE, row.names=FALSE, qmethod="double")