#########################################################################
## Permuted-block Randomization Scheme   N=80 patients grand total
## Generate N=160 to account for loss and to allow a later decision
## to recruit fewer in the high BMI group and more in low BMI.
##
## Stratified randomization
## 2 strata (D2R, DAT imaging), 2 treatment arms (Insulin; No Insulin)
## 40 patients for each stratum: block sizes randomly chosen among 2,3 *2
##
## Update: Imaging modalities (D2R, DAT) need to be treated as strata so
## that these can be filled according to lab logistics
## Add Doubly-labeled H2O as if additional treatment factor
##
## Usage:
# R --no-save --no-restore < randomization.r && gnumeric rTable.csv && cp -p rTable.csv /tmp
# This is attached to http://biostat.mc.vanderbilt.edu/RS
#########################################################################
 
require(blockrand)
type <- c('dryrun', 'real')[2]
seed <- 1317
approach <- c('image bmi strat', 'bmi strat image tx',
              'image bmi strat h2o tx')[3]
options(width=90)
sink('randomization.lst')
cat('Approach:', approach, '      Type:', type, '      seed:', seed, '\n\n')

## Old approach using 2 stratification factors
if(approach == 'image bmi strat')
  {
    nPerStrat <- 40
    strat  <- c('D2R BMI<35', 'D2R BMI>=35', 'DAT BMI<35', 'DAT BMI>=35')
    stratn <- paste('IDM', c('A','B','C','D'), sep='')
    
    ## make a randomization table for each strata, then stack
    set.seed(seed)
    d <- NULL
    for(i in 1:length(strat))
      {
        r <- blockrand(n=nPerStrat, levels=c('Insulin', 'No Insulin'),
                       id.prefix=stratn[i], block.prefix='B', stratum=strat[i],
                       block.sizes = 2:3)
        d <- rbind(d, r)
      }
  }

# Note: number of patients in the strata may be different.
# This is because of the use of random block sizes.

#---------------------------------------------------------------------

if(approach == 'bmi strat image tx')
  {
    nPerStrat <- 80
    strat  <- c('BMI<35', 'BMI>=35')
    stratn <- c('L','H')
    
    ## make a randomization table for each strata, then stack
    set.seed(11)  # used 11 for test randomization
    d <- NULL
    for(i in 1:length(strat))
      {
        r <- blockrand(n=nPerStrat,
                       levels=c('Insulin D2R', 'No Insulin D2R',
                         'Insulin DAT', 'No Insulin DAT'),
                       id.prefix=stratn[i], block.prefix='B', stratum=strat[i],
                       block.sizes = 2:3)
        d <- rbind(d, r)
      }
  }

if(approach == 'image bmi strat h2o tx')
  {
    strat  <- c('D2R BMI<35', 'D2R BMI>=35', 'DAT BMI<35', 'DAT BMI>=35')
    stratn <- paste('IDM', c('A','B','C','D'), sep='')
    nPerStrat <- 40
    
    ## make a randomization table for each strata, then stack
    set.seed(seed)
    d <- NULL
    for(i in 1:length(strat))
      {
        r <- blockrand(n=nPerStrat,
                       levels=c('Insulin, DL H2O',   'Insulin, No DL H2O',
                             'No Insulin, DL H2O','No Insulin, No DL H2O'),
                       id.prefix=stratn[i], block.prefix='B', stratum=strat[i],
                       block.sizes = c(2,2))
        d <- rbind(d, r)
      }
  }

cat('\nStratum x Treatment Frequencies     DL = double labeled\n\n')
with(d, table(stratum, treatment))
## with(d, table(paste(stratum, block.id), treatment))

cat('\n\nSubject Randomizations\n\n')
print(d)
cat('\n\nNote: Block IDs and counts will not be shown in real randomizations\n')
sink()

d <- d[c('id', 'stratum', 'treatment')]
write.table(d, 'rTable.csv', sep=",", row.names=FALSE, col.names=TRUE)
## system('gnumeric rTable.csv &')

# Note: number of patients in the strata may be different.
# This is because of the use of random block sizes.