###     A programme to obtain the power of parameters in a
#             XC 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 <- 20
n1high <- 100
n1step <- 20
n2low <- 10
n2high <- 30
n2step <- 10
npred <- 0
randsize <- 2
rand1size <- 1
rand2size <- 1
repliclow <- 1
replichigh <- 1
replicstep <- 1
beta <- c(0.500)
betasize <- length(beta)
effectbeta <- abs(beta)
sigma2u <- matrix(c(1.200), rand1size, rand1size)
sigma2v <- matrix(c(0.400), rand2size, rand2size)
sigmae <- sqrt(8.000)
n1range <- seq(n1low, n1high, n1step)
n2range <- seq(n2low, n2high, n2step)
n1size <- length(n1range)
n2size <- length(n2range)
reprange <- seq(repliclow, replichigh, replicstep)
repsize <- length(reprange)
totalsize <- n1size*n2size*repsize
finaloutput <- matrix(0, totalsize, 6*betasize)
rowcount <- 1
##~~~~~~~~~~~~~~~~~~          Creating grouped data       ~~~~~~~~~~~~~~~~~~##

fixname <- "x0"
fixform <- "1"
rand1form <- "(1|l1id)"
rand2form <- "(1|l2id)"
expression <- paste(c(fixform, rand1form, rand2form), collapse="+")
modelformula <- formula(paste("y ~", expression))
data <- vector("list", 3+length(fixname))
names(data) <- c("l1id", "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)) {
    for (replication in seq(repliclow, replichigh, replicstep)) {

      l1id <- rep(c(1:n1), each=n2*replication)
      l2id <- rep(c(1:n2), times=n1, each=replication)
      length <- n1*n2*replication
      x <- matrix(1, length, betasize)
   	  z1 <- matrix(1, length, rand1size)
      z2 <- matrix(1, length, rand2size)
      powaprox[1:betasize] <- rep(0, betasize)
      powsde <- powaprox
      sdepower <- matrix(0, betasize, simus)
      ###-----------------------     Simulation step    -----------------------###

      cat("Start of simulation for ", n2, " XC2, ", n1, " XC1 and ", replication, " 1-level units\n")
      for (iter in 1:simus) {

        if (iter/10 == floor(iter/10)) {
          cat(" Iteration remain=", simus-iter, "\n")
        }
        ##---------------------------------------------------------------------## 
        e <- rnorm(length, 0, sigmae)
        u <- mvrnorm(n1, rep(0, rand1size), sigma2u)
        v <- mvrnorm(n2, rep(0, rand2size), sigma2v)
        fixpart <- x*beta
        rand1part <- rowSums(z1*u[l1id, ])
        rand2part <- rowSums(z2*v[l2id, ])
        y <- fixpart+rand1part+rand2part+e

        ##-------------------        Inputs for model fitting       ---------------##

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

        (fitmodel <- lmer(modelformula, data, REML=FALSE))

        #~~~~~~~~~~         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] <- as.numeric(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)

      }

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

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

finaloutput <- as.data.frame(round(finaloutput, 3))
output <- data.frame(cbind(rep(n2range, each=n1size*repsize), rep(n1range, n2size, each=repsize), rep(reprange,n1size*n2size), finaloutput))
names(output) <- c("#XC2", "#XC1", "#repeat", "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")