source("MAMMA.R")

####################################################################################
# Computing Expected Standard Errors for Sample Size Determination for Paired Data #
####################################################################################

# model parameters

mu.1 <- 311.8
mu.2 <- 323.2
sigmasq.1 <- 2204.4
sigmasq.2 <- 2198.7
rho <- 0.92

param <- c(mu.1, mu.2, log(sigmasq.1), log(sigmasq.2), atanh(rho))
names(param) <- c("mu.1", "mu.2", "lsigmasq.1", "lsigmasq.2", "zrho")

# probabilty for TDI and confidence level

prob <- 0.90
conf.level <- 0.95 

# true values for agreement measures

true.meas.bvn <- true.measures.bvn(param, prob)

# > true.meas.bvn
# tdi        ccc 
# 36.0764417  0.8936234 
# > 

# > c(log(true.meas.bvn[1]), atanh(true.meas.bvn[2]))
# tdi      ccc 
# 3.585640 1.439631 
# > 


# need to do simulations for various values of number of subjects n

# repeat the following for several n (see below) and save results

n <- 100 # number of subjects
nsim <- 500 # number Monte Carlo repetitions
bds.bvn <- replicate(nsim, simulate.bounds.bvn())
write.csv(t(bds.bvn), paste0("nsub=", n, ".csv"), row.names=F)

# read these results

n30 <- read.table("nsub=30.csv", header = T, sep = ",")
n40 <- read.table("nsub=40.csv", header = T, sep = ",")
n50 <- read.table("nsub=50.csv", header = T, sep = ",")
n60 <- read.table("nsub=60.csv", header = T, sep = ",")
n70 <- read.table("nsub=70.csv", header = T, sep = ",")
n80 <- read.table("nsub=80.csv", header = T, sep = ",")
n90 <- read.table("nsub=90.csv", header = T, sep = ",")
n100 <- read.table("nsub=100.csv", header = T, sep = ",")
n125 <- read.table("nsub=125.csv", header = T, sep = ",")
n150 <- read.table("nsub=150.csv", header = T, sep = ",")
n175 <- read.table("nsub=175.csv", header = T, sep = ",")
n200 <- read.table("nsub=200.csv", header = T, sep = ",")

# compute the means

n30.means <- colMeans(n30[, c("ltdi", "se.ltdi", "zccc", "se.zccc")])
n40.means <- colMeans(n40[, c("ltdi", "se.ltdi", "zccc", "se.zccc")])
n50.means <- colMeans(n50[, c("ltdi", "se.ltdi", "zccc", "se.zccc")])
n60.means <- colMeans(n60[, c("ltdi", "se.ltdi", "zccc", "se.zccc")])
n70.means <- colMeans(n70[, c("ltdi", "se.ltdi", "zccc", "se.zccc")])
n80.means <- colMeans(n80[, c("ltdi", "se.ltdi", "zccc", "se.zccc")])
n90.means <- colMeans(n90[, c("ltdi", "se.ltdi", "zccc", "se.zccc")])
n100.means <- colMeans(n100[, c("ltdi", "se.ltdi", "zccc", "se.zccc")])
n125.means <- colMeans(n125[, c("ltdi", "se.ltdi", "zccc", "se.zccc")])
n150.means <- colMeans(n150[, c("ltdi", "se.ltdi", "zccc", "se.zccc")])
n175.means <- colMeans(n175[, c("ltdi", "se.ltdi", "zccc", "se.zccc")])
n200.means <- colMeans(n200[, c("ltdi", "se.ltdi", "zccc", "se.zccc")])

# set up a data frame to make plotting easier

all.means <- rbind(n30.means, n40.means, n50.means, n60.means,  
	n70.means, n80.means, n90.means, n100.means, n125.means,
	n150.means, n175.means, n200.means)

all.means <- data.frame(n = c(seq(f = 30, t = 100, b = 10), 125, 150, 175, 200), all.means)

# > str(all.means)
# 'data.frame':	12 obs. of  5 variables:
 # $ n      : num  30 40 50 60 70 80 90 100 125 150 ...
 # $ ltdi   : num  3.58 3.59 3.6 3.58 3.59 ...
 # $ se.ltdi: num  0.1198 0.1045 0.0939 0.0857 0.0794 ...
 # $ zccc   : num  1.41 1.4 1.41 1.42 1.42 ...
 # $ se.zccc: num  0.171 0.149 0.134 0.122 0.113 ...
# > 

# > all.means
             # n     ltdi    se.ltdi     zccc    se.zccc
# n30.means   30 3.583338 0.11975577 1.410719 0.17074261
# n40.means   40 3.593061 0.10448647 1.403952 0.14900380
# n50.means   50 3.599496 0.09389555 1.407199 0.13375491
# n60.means   60 3.583962 0.08573655 1.421095 0.12232680
# n70.means   70 3.588496 0.07939987 1.419575 0.11341482
# n80.means   80 3.593901 0.07431427 1.419237 0.10617649
# n90.means   90 3.595481 0.07017993 1.421348 0.10026461
# n100.means 100 3.591873 0.06651463 1.424399 0.09510179
# n125.means 125 3.590974 0.05968996 1.422572 0.08524847
# n150.means 150 3.593019 0.05457199 1.426532 0.07793687
# n175.means 175 3.598210 0.05048631 1.421745 0.07213936
# n200.means 200 3.596789 0.04731258 1.426630 0.06757045
# >


# plot the results

par(mfrow=c(1,2))
plot(all.means$n, all.means$se.ltdi, type = "l", 
	panel.first = grid(lty = "solid", col = gray(0.9)), xlab = "$n$", 
	ylab = "SE of estimated $\log(TDI)")

plot(all.means$n, all.means$se.zccc, type = "l", 
	panel.first = grid(lty = "solid", col = gray(0.9)), xlab = "$n$", 
	ylab = "SE of estimated z(CCC)")

######################################################################################