# install.packages("../PUMP_1.0.2.tar.gz", repos = NULL, type = "source")

library(PUMP)
library(dplyr)
library(ggplot2)
library(knitr)
library(kableExtra)
library(tidyr)

set.seed(0440044)

#------------------------------------------------------------#
# Calculating MDES
#------------------------------------------------------------#

m <- pump_mdes(d_m = "d3.2_m3fc2rc", MTP = "HO", target.power = 0.80,
    power.definition = "D1indiv", M = 5, J = 3, K = 21, nbar = 258, Tbar = 0.50,
    alpha = 0.05, numCovar.1 = 5, numCovar.2 = 3, R2.1 = 0.1, R2.2 = 0.7,
    ICC.2 = 0.05, ICC.3 = 0.4, rho = 0.4
)

# for table in manuscript
knitr::kable(m, digits = 3, booktabs = TRUE, position = "h!",
  caption = "MDES Estimate") %>%
kableExtra::kable_styling(position = "center")

m2 <- update(m, power.definition = "min1")
print(m2)

m3 <- update(m2, numZero = 2)
print(m3)


#------------------------------------------------------------#
# Determining necessary sample size
#------------------------------------------------------------#

smp <- pump_sample(d_m = "d3.2_m3fc2rc", MTP = "HO", typesample = "K",
    target.power = 0.80, power.definition = "min1", tol = 0.01, MDES = 0.10,
    M = 5, nbar = 258, J = 3, Tbar = 0.50, alpha = 0.05, numCovar.1 = 5,
    numCovar.2 = 3, R2.1 = 0.1, R2.2 = 0.7, ICC.2 = 0.05, ICC.3 = 0.40,
    rho = 0.4)
print(smp)

p_check <- update(smp, type = "power", tnum = 50000, long.table = TRUE)

knitr::kable(p_check, digits = 2, booktabs = TRUE, row.names = FALSE,
  position = "h!", caption = "Power table") %>%
kableExtra::kable_styling(position = "center")


plot( smp )

#------------------------------------------------------------#
# Comparing adjustment procedures
#------------------------------------------------------------#

p2 <- update(p_check, MTP = c("BF", "HO", "WY-SD"), tnum = 5000,
  parallel.WY.cores = 1)
plot(p2)

#------------------------------------------------------------#
# Exploring sensitivity to design parameters
#------------------------------------------------------------#

pow <- update(p_check, tnum = 10000)

#-------------------------------------#
# Exploring power with update()
#-------------------------------------#

p_b <- update(pow, ICC.2 = 0.20, ICC.3 = 0.25)
print(p_b)

p_d <- update(pow, R2.1 = c(0.1, 0.3, 0.1, 0.2, 0.2),
  R2.2 = c(0.4, 0.8, 0.3, 0.2, 0.2))
print(p_d)
summary(p_d)

#-------------------------------------#
# Exploring the impact of the ICC
#-------------------------------------#

gridICC <- update_grid(pow, ICC.2 = seq(0, 0.3, 0.05),
  ICC.3 = seq(0, 0.60, 0.20))

plot(gridICC, power.definition = "min1")


#-------------------------------------#
# Exploring the impact of rho
#-------------------------------------#

gridRho <- update_grid(pow, MTP = c("BF", "WY-SD"),
  rho = seq( 0, 0.9, by = 0.15 ), tnum = 500, B = 3000)

plot( gridRho )

#-------------------------------------#
# Exploring the impact of null outcomes
#-------------------------------------#

gridZero <- update_grid(pow, numZero = 0:4, M = 5)
plot(gridZero, nrow = 1)


#-------------------------------------#
# Checking the correlation between test statistics
#-------------------------------------#

covariate.corr.matrix <- gen_corr_matrix(M = 5, rho.scalar = 1)
cor.tstat <- check_cor(pow, rho.V = covariate.corr.matrix,
    rho.X = covariate.corr.matrix, rho.C = covariate.corr.matrix, n.sims = 500)
est.cor <- mean(cor.tstat[lower.tri(cor.tstat)])
print( est.cor )