source("MAMMA.R") #################################################################################### # Computing Expected Standard Errors for Sample Size Determination for # (Unlinked) Repeated Measurements Data # #################################################################################### # set model parameters alpha <- 11.35 mu.b <- 311.81 lsigmasq.b <- 7.61 lsigmasq.bI <- 4.57 lsigmasq.e1 <- 4.45 lsigmasq.e2 <- 4.38 param <- c(alpha, mu.b, lsigmasq.b, lsigmasq.bI, lsigmasq.e1, lsigmasq.e2) names(param) <- c("alpha", "mu.b", "lsigmasq.b", "lsigmasq.bI", "lsigmasq.e1", "lsigmasq.e2") true.meas <- true.measures(param, prob) # > true.meas # tdi tdi.1 tdi.2 ccc ccc.1 ccc.2 # 36.2734171 21.5252140 20.7848632 0.8922672 0.9610866 0.9636217 # > # > c(log(true.meas[1:3]), atanh(true.meas[4:6])) # tdi tdi.1 tdi.2 ccc ccc.1 ccc.2 # 3.591085 3.069225 3.034225 1.432939 1.959958 1.994287 # > # set probability for TDI and confidence level prob <- 0.90 conf.level <- 0.95 # need to do simulations for various combinations of number of subjects n # and number of replications m # repeat the following for several (n, m) combinations (see below) and save results n <- 100 # number of subjects m <- 2 # number of replications nsim <- 500 # number Monte Carlo repetitions bds <- replicate(nsim, simulate.bounds()) write.csv(t(bds), paste0("nsub=", n, "-nrep=", m, ".csv"), row.names=F) # read these results n30m2 <- read.table("nsub=30-nrep=2.csv", header = T, sep = ",") n30m3 <- read.table("nsub=30-nrep=3.csv", header = T, sep = ",") n40m2 <- read.table("nsub=40-nrep=2.csv", header = T, sep = ",") n40m3 <- read.table("nsub=40-nrep=3.csv", header = T, sep = ",") n50m2 <- read.table("nsub=50-nrep=2.csv", header = T, sep = ",") n50m3 <- read.table("nsub=50-nrep=3.csv", header = T, sep = ",") n60m2 <- read.table("nsub=60-nrep=2.csv", header = T, sep = ",") n60m3 <- read.table("nsub=60-nrep=3.csv", header = T, sep = ",") n70m2 <- read.table("nsub=70-nrep=2.csv", header = T, sep = ",") n70m3 <- read.table("nsub=70-nrep=3.csv", header = T, sep = ",") n80m2 <- read.table("nsub=80-nrep=2.csv", header = T, sep = ",") n80m3 <- read.table("nsub=80-nrep=3.csv", header = T, sep = ",") n90m2 <- read.table("nsub=90-nrep=2.csv", header = T, sep = ",") n90m3 <- read.table("nsub=90-nrep=3.csv", header = T, sep = ",") n100m2 <- read.table("nsub=100-nrep=2.csv", header = T, sep = ",") n100m3 <- read.table("nsub=100-nrep=3.csv", header = T, sep = ",") n125m2 <- read.table("nsub=125-nrep=2.csv", header = T, sep = ",") n125m3 <- read.table("nsub=125-nrep=3.csv", header = T, sep = ",") n150m2 <- read.table("nsub=150-nrep=2.csv", header = T, sep = ",") n150m3 <- read.table("nsub=150-nrep=3.csv", header = T, sep = ",") n175m2 <- read.table("nsub=175-nrep=2.csv", header = T, sep = ",") n175m3 <- read.table("nsub=175-nrep=3.csv", header = T, sep = ",") n200m2 <- read.table("nsub=200-nrep=2.csv", header = T, sep = ",") n200m3 <- read.table("nsub=200-nrep=3.csv", header = T, sep = ",") # > colnames(n50m2) # [1] "ltdi" "se.ltdi" "ucb.ltdi" "tdi" "ucb.tdi" "ltdi.1" # [7] "se.ltdi.1" "ucb.ltdi.1" "tdi.1" "ucb.tdi.1" "ltdi.2" "se.ltdi.2" # [13] "ucb.ltdi.2" "tdi.2" "ucb.tdi.2" "zccc" "se.zccc" "lcb.zccc" # [19] "ccc" "lcb.ccc" "zccc.1" "se.zccc.1" "lcb.zccc.1" "ccc.1" # [25] "lcb.ccc.1" "zccc.2" "se.zccc.2" "lcb.zccc.2" "ccc.2" "lcb.ccc.2" # > # compute the means n30m2.means <- colMeans(n30m2[, c( "ltdi", "se.ltdi", "ltdi.1", "se.ltdi.1", "ltdi.2", "se.ltdi.2", "zccc", "se.zccc", "zccc.1", "se.zccc.1", "zccc.2", "se.zccc.2")]) n30m3.means <- colMeans(n30m3[, c( "ltdi", "se.ltdi", "ltdi.1", "se.ltdi.1", "ltdi.2", "se.ltdi.2", "zccc", "se.zccc", "zccc.1", "se.zccc.1", "zccc.2", "se.zccc.2")]) n40m2.means <- colMeans(n40m2[, c( "ltdi", "se.ltdi", "ltdi.1", "se.ltdi.1", "ltdi.2", "se.ltdi.2", "zccc", "se.zccc", "zccc.1", "se.zccc.1", "zccc.2", "se.zccc.2")]) n40m3.means <- colMeans(n40m3[, c( "ltdi", "se.ltdi", "ltdi.1", "se.ltdi.1", "ltdi.2", "se.ltdi.2", "zccc", "se.zccc", "zccc.1", "se.zccc.1", "zccc.2", "se.zccc.2")]) n50m2.means <- colMeans(n50m2[, c( "ltdi", "se.ltdi", "ltdi.1", "se.ltdi.1", "ltdi.2", "se.ltdi.2", "zccc", "se.zccc", "zccc.1", "se.zccc.1", "zccc.2", "se.zccc.2")]) n50m3.means <- colMeans(n50m3[, c( "ltdi", "se.ltdi", "ltdi.1", "se.ltdi.1", "ltdi.2", "se.ltdi.2", "zccc", "se.zccc", "zccc.1", "se.zccc.1", "zccc.2", "se.zccc.2")]) n60m2.means <- colMeans(n60m2[, c( "ltdi", "se.ltdi", "ltdi.1", "se.ltdi.1", "ltdi.2", "se.ltdi.2", "zccc", "se.zccc", "zccc.1", "se.zccc.1", "zccc.2", "se.zccc.2")]) n60m3.means <- colMeans(n60m3[, c( "ltdi", "se.ltdi", "ltdi.1", "se.ltdi.1", "ltdi.2", "se.ltdi.2", "zccc", "se.zccc", "zccc.1", "se.zccc.1", "zccc.2", "se.zccc.2")]) n70m2.means <- colMeans(n70m2[, c( "ltdi", "se.ltdi", "ltdi.1", "se.ltdi.1", "ltdi.2", "se.ltdi.2", "zccc", "se.zccc", "zccc.1", "se.zccc.1", "zccc.2", "se.zccc.2")]) n70m3.means <- colMeans(n70m3[, c( "ltdi", "se.ltdi", "ltdi.1", "se.ltdi.1", "ltdi.2", "se.ltdi.2", "zccc", "se.zccc", "zccc.1", "se.zccc.1", "zccc.2", "se.zccc.2")]) n80m2.means <- colMeans(n80m2[, c( "ltdi", "se.ltdi", "ltdi.1", "se.ltdi.1", "ltdi.2", "se.ltdi.2", "zccc", "se.zccc", "zccc.1", "se.zccc.1", "zccc.2", "se.zccc.2")]) n80m3.means <- colMeans(n80m3[, c( "ltdi", "se.ltdi", "ltdi.1", "se.ltdi.1", "ltdi.2", "se.ltdi.2", "zccc", "se.zccc", "zccc.1", "se.zccc.1", "zccc.2", "se.zccc.2")]) n90m2.means <- colMeans(n90m2[, c( "ltdi", "se.ltdi", "ltdi.1", "se.ltdi.1", "ltdi.2", "se.ltdi.2", "zccc", "se.zccc", "zccc.1", "se.zccc.1", "zccc.2", "se.zccc.2")]) n90m3.means <- colMeans(n90m3[, c( "ltdi", "se.ltdi", "ltdi.1", "se.ltdi.1", "ltdi.2", "se.ltdi.2", "zccc", "se.zccc", "zccc.1", "se.zccc.1", "zccc.2", "se.zccc.2")]) n100m2.means <- colMeans(n100m2[, c( "ltdi", "se.ltdi", "ltdi.1", "se.ltdi.1", "ltdi.2", "se.ltdi.2", "zccc", "se.zccc", "zccc.1", "se.zccc.1", "zccc.2", "se.zccc.2")]) n100m3.means <- colMeans(n100m3[, c( "ltdi", "se.ltdi", "ltdi.1", "se.ltdi.1", "ltdi.2", "se.ltdi.2", "zccc", "se.zccc", "zccc.1", "se.zccc.1", "zccc.2", "se.zccc.2")]) n125m2.means <- colMeans(n125m2[, c( "ltdi", "se.ltdi", "ltdi.1", "se.ltdi.1", "ltdi.2", "se.ltdi.2", "zccc", "se.zccc", "zccc.1", "se.zccc.1", "zccc.2", "se.zccc.2")]) n125m3.means <- colMeans(n125m3[, c( "ltdi", "se.ltdi", "ltdi.1", "se.ltdi.1", "ltdi.2", "se.ltdi.2", "zccc", "se.zccc", "zccc.1", "se.zccc.1", "zccc.2", "se.zccc.2")]) n150m2.means <- colMeans(n150m2[, c( "ltdi", "se.ltdi", "ltdi.1", "se.ltdi.1", "ltdi.2", "se.ltdi.2", "zccc", "se.zccc", "zccc.1", "se.zccc.1", "zccc.2", "se.zccc.2")]) n150m3.means <- colMeans(n150m3[, c( "ltdi", "se.ltdi", "ltdi.1", "se.ltdi.1", "ltdi.2", "se.ltdi.2", "zccc", "se.zccc", "zccc.1", "se.zccc.1", "zccc.2", "se.zccc.2")]) n175m2.means <- colMeans(n175m2[, c( "ltdi", "se.ltdi", "ltdi.1", "se.ltdi.1", "ltdi.2", "se.ltdi.2", "zccc", "se.zccc", "zccc.1", "se.zccc.1", "zccc.2", "se.zccc.2")]) n175m3.means <- colMeans(n175m3[, c( "ltdi", "se.ltdi", "ltdi.1", "se.ltdi.1", "ltdi.2", "se.ltdi.2", "zccc", "se.zccc", "zccc.1", "se.zccc.1", "zccc.2", "se.zccc.2")]) n200m2.means <- colMeans(n200m2[, c( "ltdi", "se.ltdi", "ltdi.1", "se.ltdi.1", "ltdi.2", "se.ltdi.2", "zccc", "se.zccc", "zccc.1", "se.zccc.1", "zccc.2", "se.zccc.2")]) n200m3.means <- colMeans(n200m3[, c( "ltdi", "se.ltdi", "ltdi.1", "se.ltdi.1", "ltdi.2", "se.ltdi.2", "zccc", "se.zccc", "zccc.1", "se.zccc.1", "zccc.2", "se.zccc.2")]) # set up a data frame to make plotting easier all.means <- rbind(n30m2.means, n30m3.means, n40m2.means, n40m3.means, n50m2.means, n50m3.means, n60m2.means, n60m3.means, n70m2.means, n70m3.means, n80m2.means, n80m3.means, n90m2.means, n90m3.means, n100m2.means, n100m3.means, n125m2.means, n125m3.means, n150m2.means, n150m3.means, n175m2.means, n175m3.means, n200m2.means, n200m3.means) all.means <- data.frame(n = rep(c(seq(f = 30, t = 100, b = 10), 125, 150, 175, 200), each = 2), m = rep(c(2, 3), 12), all.means) # > str(all.means) # 'data.frame': 24 obs. of 14 variables: # $ n : num 30 30 40 40 50 50 60 60 70 70 ... # $ m : num 2 3 2 3 2 3 2 3 2 3 ... # $ ltdi : num 3.58 3.58 3.58 3.59 3.58 ... # $ se.ltdi : num 0.0992 0.091 0.0864 0.0799 0.0773 ... # $ ltdi.1 : num 3.05 3.06 3.05 3.06 3.06 ... # $ se.ltdi.1: num 0.129 0.0913 0.1117 0.079 0.1 ... # $ ltdi.2 : num 3.02 3.03 3.03 3.03 3.03 ... # $ se.ltdi.2: num 0.129 0.0913 0.1118 0.079 0.0999 ... # $ zccc : num 1.41 1.41 1.42 1.41 1.42 ... # $ se.zccc : num 0.159 0.154 0.138 0.134 0.123 ... # $ zccc.1 : num 1.95 1.93 1.95 1.95 1.94 ... # $ se.zccc.1: num 0.179 0.153 0.155 0.133 0.138 ... # $ zccc.2 : num 1.97 1.96 1.97 1.98 1.97 ... # $ se.zccc.2: num 0.179 0.154 0.155 0.133 0.138 ... # > # plot the results par(mfrow=c(1,2)) plot(se.ltdi ~ n, subset = (m == 2), data = all.means, ylim = c(0.035, 0.13), type = "l", col = "black", ylab = "SE of estimated $\\log$(TDI)", xlab = "$n$") lines(se.ltdi ~ n, subset = (m == 3), data = all.means, lty = 1, col = gray(0.5)) lines(se.ltdi.1 ~ n, subset = (m == 2), data = all.means, lty = 2, col = "black") lines(se.ltdi.1 ~ n, subset = (m == 3), data = all.means, lty = 2, col = gray(0.5)) grid(lty = "solid", col = gray(0.9)) plot(se.zccc ~ n, subset = (m == 2), data = all.means, ylim = c(0.06, 0.18), type = "l", col = "black", ylab = "SE of estimated $z$(CCC)", xlab = "$n$") lines(se.zccc ~ n, subset = (m == 3), data = all.means, lty = 1, col = gray(0.5)) lines(se.zccc.1 ~ n, subset = (m == 2), data = all.means, lty = 2, col = "black") lines(se.zccc.1 ~ n, subset = (m == 3), data = all.means, lty = 2, col = gray(0.5)) grid(lty = "solid", col = gray(0.9)) #########################################################################################