###     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 <- 20
n2high <- 300
n2step <- 20
npred <- 2
randsize <- 1
beta <- c(-0.226000, 0.257000, 0.146000)
betasize <- length(beta)
effectbeta <- abs(beta)
randcolumn <- 0
meanpred <- c(0, 0.600000, 0.462000)
varpred <- matrix(c(0.120000, 0.000000, 0.000000, 0.000000), npred, npred)
varpred2 <- matrix(c(0.125000, 0.045000, 0.045000, 0.249000), npred, npred)
sigma2u <- matrix(c(0.156000), 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)
rowcount <- 1
##-----------------        Inputs for model fitting       -----------------##

fixname <- c("x0", "x1", "x2")
fixform <- "1+x1+x2"
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) <- c("b0", "b1", "b2")
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       +++++++++++++++++++  ##

      micpred <- mvrnorm(length, meanpred[-1], varpred)
      macpred <- mvrnorm(n2, rep(0, npred), varpred2)
      x[, (2:dim(x)[2])] <- micpred+macpred[l2id, ]
      ##--------------------------------------------------------------## 
      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]
      ###~~~~~~~~~~      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] <- 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", "zLb1", "zpb1", "zUb1", "sLb1", "spb1", "sUb1", "zLb2", "zpb2", "zUb2", "sLb2", "spb2", "sUb2")
write.table(output, "powerout.txt", sep="\t ", quote=FALSE, eol="\n", dec=".", col.names=TRUE, row.names=FALSE, qmethod="double")