###     A programme to obtain the power of parameters in 2 level
#       balanced model  with  Normal 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 <- 40
n1high <- 40
n1step <- 10
n2low <- 60
n2high <- 300
n2step <- 20
npred <- 3
randsize <- 1
beta <- c(-0.228000, 0.262000, 0.191000, 0.123000)
betasize <- length(beta)
effectbeta <- abs(beta)
randcolumn <- 0
xprob <- c(0, 0.600000, 0.000000, 0.000000)
meanpred <- c(0, 0.000000, 0.150000, 0.300000)
varpred <- c(0, 0.000000, 0.000000, 0.000000)
varpred2 <- c(0, 0.000000, 0.130000, 0.210000)
sigma2u <- matrix(c(0.155000), randsize, randsize)
sigmae <- sqrt(0.839000)
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+1)
rowcount <- 1
##-----------------        Inputs for model fitting       -----------------##

fixname <- c("x0", "x1", "x2", "x3")
fixform <- "1+x1+x2+x3"
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)
modelformula1 <- formula(y ~ 1+x1+(1|l2id))
devtestsim <- rep(0, simus) 

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

powaprox <- vector("list", betasize)
names(powaprox) <- c("b0", "b1", "b2", "b3")
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")
      }
      ## +++++++++++++++++++       Set up X matrix       +++++++++++++++++++  ##

      x[, 2] <- rbinom(length, 1, 0.5)
      macpred <- rmultinom(n2, 1, c(0.15, 0.30, 0.55))
      x[, 3] <- macpred[1, ][l2id]
      x[, 4] <- macpred[2, ][l2id] 
      x[, 2] <- x[,4]+x[, 2]*(x[, 3]==0&x[, 4]==0) 
      ##--------------------------------------------------------------## 
      e <- rnorm(length, 0, sigmae)
      u <- mvrnorm(n2, rep(0, randsize), sigma2u)
      fixpart <- x %*% beta
      randpart <- rowSums(z*u[l2id, ])
      y <- fixpart+randpart+e
      ##-------------------        Inputs for model fitting       ---------------##

      data$l2id <- l2id
      data$y <- y
      data$x0 <- x[, 1]
      data$x1 <- x[, 2]
      data$x2 <- x[, 3]
      data$x3 <- x[, 4]
      ###~~~~~~~~~~      Fitting the model using lmer funtion    ~~~~~~~~~~###

      (fitmodel <- lmer(modelformula, data, REML=FALSE))
      (fitmodel1 <- lmer(modelformula1, data, REML=FALSE))
      devtestsim[iter] <- deviance(fitmodel1) - deviance(fitmodel)
      ######~~~~~~~~~~   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)
      finaloutput[rowcount, 6*l+1] <- mean(devtestsim > 7.38)

    } 

    ###~~~~~~~~~~    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", "zLb1", "zpb1", "zUb1", "sLb1", "spb1", "sUb1", "zLb2", "zpb2", "zUb2", "sLb2", "spb2", "sUb2", "zLb3", "zpb3", "zUb3", "sLb3", "spb3", "sUb3", "devtest")
write.table(output, "powerout.txt", sep="\t ", quote=FALSE, eol="\n", dec=".", col.names=TRUE, row.names=FALSE, qmethod="double")