#############################
### Randomization schemes ###
#############################

# (ii)-(a) We expect to randomize 100 patients
# What's the probability of getting equal sample size?
?dbinom
N = 100
dbinom(x = N/2, size = N,p = 0.5)

# (ii)-(b) Simulate an example of simple randomization
p = 0.5
N = 100
data = data.frame(id = paste("id",1:N,sep = ""), trt = NA)
head(data)

# Set the seed if you would like to reproduce your results
set.seed(125)

# Simulate 100 uniform random variable
u = runif(N)
hist(u)

# Make assignments
data$trt[u<0.5] = "A"
data$trt[u>0.5] = "B"
data[1:10,]

table(data$trt)
table(data$trt)/N

# (ii)-c: Repeating simple randomization 1000 times
K = 1000
set.seed(2342)

results = data.frame(replication = 1:K, NumberA = NA, NumberB = NA)

for (i in 1:K)
{
  data = data.frame(id = paste("id",1:N,sep = ""), trt = NA)
  u = runif(N)
  
  # Make assignments
  data$trt[u<0.5] = "A"
  data$trt[u>0.5] = "B"
  
  results[i,"NumberA"] = table(data$trt)["A"]
  results[i,"NumberB"] = table(data$trt)["B"]
}
head(results)

# Number of replications yielding imbalance greater than 60/40
table(results[,2] > 60 | results[,2] < 40)/K

# See means and variances
# Note the equal number on each treatment in the long run 
apply(results[,2:3],2,mean)
apply(results[,2:3],2,var)

setwd("C:/Clinical trials/R labs/lab3/graphs")
# Notice that Na converges to 50
pdf("simple_rand.pdf",height = 5.5,width = 5.5)
plot(cumsum(results[,2])/1:K,type = "l",xlab = "Replication index",
     ylab = "Mean number of patients on treatment A")
dev.off()

# (iii) Randomization in blocks
# (a) block size 4 and two treatments A and B
# Number of permutations
choose(4,2)

# See the permutations
library(Deducer)
t(perm(c("A","A","B","B")))

# (c) Randomize 100 subjects using random block sizes
library(blockrand)
?blockrand 

set.seed(1322)
block_rand = blockrand(n = 100, num.levels = 2, levels = c("A","B"), 
                       block.sizes = 1:4)
block_rand[block_rand$block.id %in% 1:2,]

# Actual number of randomizations made
nrow(block_rand)

# Number of patients in each group
table(block_rand$treatment[1:100])

# (iv): Blocking within strata
# (b) Suppose we stratify by sex
set.seed(2435)
male=blockrand(n=100,id.prefix='M', block.prefix='M',stratum='Male')
female=blockrand(n=100,id.prefix='F', block.prefix='F',stratum='Female')

# 55 men and 45 women at the end of the study
male = male[1:55,]
female = female[1:45,]

# Combine into one dataset
mystudy = rbind(male,female)

# Number of subject in the two groups
table(mystudy$treatment)

# Two-way contingency table 
x = table(mystudy$treatment,mystudy$stratum)
x
chisq.test(x,simulate.p.value = T,B = 10000)

# Urn schemes
#library(installr)
#updateR()
set.seed(12)
library(randomizeR)
x = udPar(N=2500, ini = 1, add = 1, groups = LETTERS[1:2])
y = genSeq(x,r = 1,seed = 12)
y
  
c(getRandList(y))[1:10]
table(c(getRandList(y)))

urnModel = function(n)
{
  na = 1
  nb = 1
  
  prob = rep(NA,n)
  trt = rep(NA,n)
  u = runif(n)
  
  for (i in 1:n)
  {
    # Probability of getting A
    prob[i] = na/(na+nb)
    
    if (u[i]<prob[i])
    {
      trt[i] = "A"
      nb = nb + 1
    } else{trt[i] = "B";na = na + 1}
  }
  
  return(data.frame(treatment = trt,probA = prob,na = na,nb = nb))
}

draws = urnModel(n = 2500)
head(draws)
table(draws$treatment)

setwd("C:/Clinical trials/R labs/lab3/graphs")

pdf("urn.pdf",height = 5.5,width = 5.5)
plot(1:2500,draws$probA,ylab = "Probability",xlab = "Number of subjects",
     main = "Probablity of getting A based on the urn model",type = "l")
abline(h = 0.5,col = "red")
dev.off()

################
### Optional ###
################

# Number of Subjects
set.seed(5)
N = 14
it = 200

# Results
results = data.frame(SR = rep(NA,it),BR = NA,urn = NA)

for (i in 1:it)
{
  ########################
  # Simple randomization #
  ########################
  data = data.frame(id = paste("id",1:N,sep = ""), trt = NA)
  
  # Simulate u to decide
  u = runif(N)
  
  # Assign patients to treatment arms
  data$trt[u<0.5] = "A"
  data$trt[u>0.5] = "B"
  
  results$SR[i] = table(data$trt)[1]/N
  
  ###########################
  # Randomization in blocks #
  ###########################
  block_rand = blockrand(n = N, num.levels = 2, levels = c("A","B"), 
                         block.sizes = 1:4)[1:N,]
  
  results$BR[i] = table(block_rand$treatment)[1]/N
  
  ##############
  # Urn scheme #
  ##############
  data = urnModel(N)
  
  results$urn[i] = table(data$treatment)[1]/N
}

# 50% chance of getting each treament in the long run
apply(results,2,mean)

# The variance is particularly high in simple randomization
apply(results,2,sd)

pdf("sim_results.pdf",height = 5.5,width = 5.5)
plot(cumsum(results$SR)/1:it,type = "l",lty = 1,col = "black",
     ylim = c(0.40,0.60),ylab = "Proportion of assigning A",xlab = "Draws")
lines(cumsum(results$BR)/1:it,lty = 2,col = "red")
lines(cumsum(results$urn)/1:it,lty = 3,col = "blue")

legend("topright",lty = 1:3,col = c("black","red","blue"),
       legend = c("Simple randomization","Randomization in blocks",
                  "Urn scheme"),bty = "n")
dev.off()