set.seed(123321)
nrep <- 10000
dat <- matrix(c(0:10, 0:10 + 1, 0:10 + 2, 0:10 + 3, 0:10 + 4), nrow=5, byrow=TRUE)
n_all <- length(dat)
n_group <- ncol(dat)
n_each <- nrow(dat)
usedgroup <- 2
usedeach <- 2
usedall <- usedgroup * usedeach

sdp <- function(x) sqrt(mean((x-mean(x))^2))

resultmean <- matrix(NA, nrep, 2)
resultse <- matrix(NA, nrep, 2)
for(i in 1:nrep) {
  # Simple Random Sampling
  index <- sample(1:n_all, usedall)
  tempdat1 <- dat[index]
  resultmean[i, 1] <- mean(tempdat1)
  resultse[i, 1] <- sd(tempdat1)/sqrt(usedall)
  
  # Multistage Random Sampling
  index_group <- sample(1:n_group, usedgroup)
  tempdat2g1 <- dat[sample(1:n_each, usedeach), index_group[1]]
  tempdat2g2 <- dat[sample(1:n_each, usedeach), index_group[2]]
  tempdat2 <- c(tempdat2g1, tempdat2g2)
  resultmean[i, 2] <- mean(tempdat2)
  resultse[i, 2] <- sd(tempdat2)/sqrt(usedall)
}

apply(resultmean, 2, mean)
apply(resultmean, 2, sdp)
apply(resultse, 2, mean)

par(mfrow=c(2,2))
hist(resultmean[,1], xlim=c(0, 14), main="Samp Dist of M (Simple)", xlab="M")
abline(v=7, col="red")
hist(resultmean[,2], xlim=c(0, 14), main="Samp Dist of M (Multistage)", xlab="M")
abline(v=7, col="red")

apply(resultmean, 2, mean) - 7 # Bias
(apply(resultmean, 2, mean) - 7)/7 # Relative Bias
apply(resultse, 2, mean) - apply(resultmean, 2, sdp) # Bias in SE
(apply(resultse, 2, mean) - apply(resultmean, 2, sdp))/apply(resultmean, 2, sdp) # Relative Bias in SE



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


set.seed(123321)
nrep <- 10000
dat <- matrix(c(0:10, 0:10 + 3, 0:10 + 4, 0:10 + 5, 0:10 + 8), nrow=5, byrow=TRUE)
n_all <- length(dat)
n_group <- ncol(dat)
n_each <- nrow(dat)
usedgroup <- 2
usedeach <- 2

sdp <- function(x) sqrt(mean((x-mean(x))^2))

resultmean <- matrix(NA, nrep, 2)
resultse <- matrix(NA, nrep, 2)
for(i in 1:nrep) {
  # Single step
  index <- sample(1:n_group, usedgroup)
  tempdat1 <- dat[3, index]
  resultmean[i, 1] <- mean(tempdat1)
  resultse[i, 1] <- sd(tempdat1)/sqrt(usedgroup)
  
  # Double Step
  tempdat2m1 <- mean(dat[sample(1:n_each, usedeach), index[1]])
  tempdat2m2 <- mean(dat[sample(1:n_each, usedeach), index[2]])
  tempdat2 <- c(tempdat2m1, tempdat2m2)
  resultmean[i, 2] <- mean(tempdat2)
  resultse[i, 2] <- sd(tempdat2)/sqrt(usedgroup)
}

apply(resultmean, 2, mean)
apply(resultmean, 2, sdp)
apply(resultse, 2, mean)

par(mfrow=c(2,1))
hist(resultmean[,1], xlim=c(2, 16), main="Samp Dist of M (Single)", xlab="M", breaks=11)
abline(v=9, col="red")
hist(resultmean[,2], xlim=c(2, 16), main="Samp Dist of M (Double)", xlab="M", breaks=11)
abline(v=9, col="red")

apply(resultmean, 2, mean) - 9 # Bias
(apply(resultmean, 2, mean) - 9)/9 # Relative Bias
apply(resultse, 2, mean) - apply(resultmean, 2, sdp) # Bias in SE
(apply(resultse, 2, mean) - apply(resultmean, 2, sdp))/apply(resultmean, 2, sdp) # Relative Bias in SE