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

set.seed(1)
siglevel <- 0.025
z1score <- abs(qnorm(siglevel))
simus <- 1000
n1low <- 50
n1high <- 1500
n1step <- 50
npred <- 1
beta <- c(-0.140000, 0.234000)
betasize <- length(beta)
effectbeta <- abs(beta)
xprob <- c(0, 0.600000)
meanpred <- c(0, 0.000000)
sigmae <- sqrt(0.985000)
n1range <- seq(n1low, n1high, n1step)
n1size <- length(n1range)
totalsize <- n1size
finaloutput <- matrix(0, totalsize, 6*betasize)
rowcount <- 1
##-----------------        Inputs for model fitting       -----------------##

fixname <- c("x0", "x1")
fixform <- "1+x1"
data <- vector("list", 1+length(fixname))
names(data) <- c("y", fixname)
modelformula <- as.formula(paste("y ~", fixform, sep=""))


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

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

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

for (n1 in seq(n1low, n1high, n1step)) {
  cat(" Start of simulation for sample sizes of ", n1, " units\n")
  length <- n1
  x <- matrix(1, length, betasize)
  sdepower <- matrix(0, betasize, simus)
  powaprox[1:betasize] <- rep(0, betasize)
  powsde <- powaprox

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

  for (iter in 1:simus) {

    if (iter/10 == floor(iter/10)){
      cat(" Iteration remain=", simus-iter, "\n")
    }
    ## +++++++++++++++++++       Data generating      +++++++++++++++++++  ##

    x[, 2] <- rbinom(length, 1, xprob[2])
    fixpart <- x %*% beta
    e <- rnorm(length, 0, sigmae)
    data$y <- fixpart+e
    ###-------------    Fitting the model using glm funtion   -------------###

    for (i in 2:length(data)) data[[i]] <- x[, i-1]
    fitmodel <- glm(modelformula, data, family=gaussian)

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

    estbeta <- coef(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    ~~~~~~~~~~###

  cat("                  power of parameter(s) for the sample size",n1,"\n") 
  rowcount <- rowcount+1
  cat("--------------------------------------------------------------------\n")
} ## end of the loop  over the sample size (n1)

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

finaloutput <- as.data.frame(round(finaloutput, 3))
output <- data.frame(cbind(n1range, finaloutput))

names(output) <- c("n", "zLb0", "zpb0", "zUb0", "sLb0", "spb0", "sUb0", "zLb1", "zpb1", "zUb1", "sLb1", "spb1", "sUb1")
write.table(output, "powerout.txt", sep="\t ", quote=FALSE, eol="\n", dec=".", col.names=TRUE, row.names=FALSE, qmethod="double")