###     A programme to obtain the power of parameters in an
#             XC unbalanced model with  Normal response
#    NB: unbalanced design based on input file
#                    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 <- 148
n1high <- 148
n1step <- 1
n2low <- 19
n2high <- 19
n2step <- 1
npred <- 0
randsize <- 2
rand1size <- 1
rand2size <- 1
repliclow <- 1
replichigh <- 1
replicstep <- 1
totalsamplelow <- 200
totalsamplehigh <- 4000
totalsamplestep <- 200
xctable <- matrix(scan("fife.txt"), n1low, n2low, byrow=TRUE)
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)
reprange <- seq(repliclow, replichigh, replicstep)
n1size <- length(n1range)
n2size <- length(n2range)
repsize <- length(reprange)
Tsamplerange <- seq(totalsamplelow, totalsamplehigh, totalsamplestep)
Tsamplesize <- length(Tsamplerange)
totalsize <- n1size*n2size*repsize*Tsamplesize
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)) {
  n2run <- n2
  for (n1 in seq(n1low, n1high, n1step)) {
    n1run <- n1
    for (replication in seq(repliclow, replichigh, replicstep)) {

      cat("Start of simulation for ", n2run, " XC2, ",n1run, " XC1 and ", replication, " 1-level units\n")
      for(totalsample in Tsamplerange) {

        cat("Start of simulation for ", totalsample, " observations\n")
        powaprox[1:betasize] <- rep(0, betasize)
        powsde <- powaprox

        sdepower <- matrix(0, betasize, simus)

        ###-----------------------     Simulation step    -----------------------###

        for (iter in 1:simus) {
          if (iter/10 == floor(iter/10)) {
            cat(" Iteration remain=", simus-iter, "\n")
          }
          ###---------   Joint samling from XC1 crossed XC2   ----###

          xcunbal <- rmultinom(1, totalsample, c(t(xctable)))
          l1id <- rep(rep(1:n1, each=n2), xcunbal)
          l2id <- rep(rep(1:n2, n1), xcunbal)
          cumxc <- c(0, cumsum(xcunbal))
          length <- replication*cumxc[length(cumxc)]
          x<-matrix(1, length, betasize)
          z1<-matrix(1, length, rand1size)
          z2<-matrix(1, length, rand2size)

          #######----------------------------------------------------######## 

          u <- mvrnorm(n1, rep(0, rand1size), sigma2u)
          v <- mvrnorm(n2run, rep(0, rand2size), sigma2v)
          fixpart <- x*beta
          rand1part <- rowSums(z1*u[l1id, ])
          rand2part <- rowSums(z2*v[l2id, ])
          e <- rnorm(length, 0, sigmae)
          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 total sample units
    } ## 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(Tsamplerange,repsize*n1size*n2size), finaloutput))
names(output) <- c("#XC2", "#XC1", "#Tsample", "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")