###############################################################################
#                                                                             #
#   Reproducible files to the paper 'Pitfalls of hypothesis tests             #
#   and model selection on bootstrap samples: causes and consequences in      #
#   biometrical applications'  (2014)                                         #
#   by Janitza, Binder and Boulesteix                                         #
#                                                                             #
#   Contact: S. Janitza <janitza@ibe.med.uni-muenchen.de>                     #
#                                                                             #
###############################################################################

# NOTE: running this code requires installation of R packages snowfall and mboost

# uncomment the following line and set the working directory where the objects shall be stored
# setwd("...") 

require(snowfall)

sfInit(parallel = TRUE, cpus = 10, type = "MPI") # if you won't do parallel computing set parallel = FALSE



# function to investigate tuning parameter selection using cross-validation on bootstrap samples

# input parameters: 
#   - seed: a seed number for reproducibility purposes (permits that the same results are obtained when running in sequential mode)
#   - p: number of predictor variables
#   - effect: effect size of influential predictors

# output: (niter x 4) matrix
#         1st column contains mstop values for original samples, 
#         2nd column contains mstop values for bootstrap samples (allowing training and test sets to overlap), 
#         3rd column contains mstop values for bootstrap samples (not allowing training and test sets to overlap),
#         4th column contains mstop values for subsamples

perform.study.CV <- function(seed, p = 200, effect = 1){
  
  require(mboost)
  
  # set up effect sizes
  beta <- rep(0, p)
  beta[((1:p) * 200/p) %in% c(1, 3, 5, 7, 9)]  <- effect
  beta[((1:p) * 200/p) %in% c(2, 4, 6, 8, 10)] <- -effect
  
  set.seed(seed)
  
  # predictors drawn from N(0,1)
  X <- matrix(rnorm(n = 100 * p), ncol = p)
  
  # response
  eta <- X %*% beta
  prob <- exp(eta)/(1 + exp(eta))
  y <- apply(prob, 1, function(z) rbinom(prob = z, n = 1, size = 1))
  
  # original sample
  data <- data.frame(id = 1:100, y = y, X = X)
  
  # bootstrap sample
  index     <- sort(sample(1:100, replace = TRUE, size = 100))
  boot.data <- data[index,]
  
  # subsample
  index.sub <- sort(sample(1:100, replace = FALSE, size = 63))
  sub.data  <- data[index.sub,]
  
  # fit models (on original sample, bootstrap sample and subsample) 
  model      <- glmboost(as.factor(y) ~ . - id, data = data,      control = boost_control(mstop = 1000), family = Binomial())
  model.boot <- glmboost(as.factor(y) ~ . - id, data = boot.data, control = boost_control(mstop = 1000), family = Binomial())
  model.sub  <- glmboost(as.factor(y) ~ . - id, data = sub.data,  control = boost_control(mstop = 1000), family = Binomial())
  
  # stratify according to person id
  cvindex <- sample(rep(1:5, length = length(unique(boot.data$id))))
  cvindex <- rep(cvindex, table(boot.data$id))
  cvindex <- sapply(1:5, function(z) as.numeric(cvindex != z))
  
  # perform 5 fold cross-validation and extract optimal number of boosting steps
  boost.steps.orig  <- mstop(cvrisk(model,      cv(model.weights(model),      type = "kfold", B = 5)))
  boost.steps.std   <- mstop(cvrisk(model.boot, cv(model.weights(model.boot), type = "kfold", B = 5)))
  boost.steps.strat <- mstop(cvrisk(model.boot, folds = cvindex))
  boost.steps.sub   <- mstop(cvrisk(model.sub,  cv(model.weights(model.sub),  type = "kfold", B = 5)))
    
  
  return(as.matrix(data.frame(original           = boost.steps.orig, 
                              bootstrap.standard = boost.steps.std, 
                              bootstrap.strat    = boost.steps.strat, 
                              subsample          = boost.steps.sub)))
}



sfExport("perform.study.CV")

# settings with weak effects
CV_setting_uncor_weak_p200 <- t(sfSapply(1:1000, function(z) perform.study.CV(seed = z, p = 200, effect = 1)))
save(CV_setting_uncor_weak_p200, file = "CV_setting_uncor_weak_p200.Rda")
CV_setting_uncor_weak_p1000 <- t(sfSapply(1:1000, function(z) perform.study.CV(seed = z, p = 1000, effect = 1)))
save(CV_setting_uncor_weak_p1000, file = "CV_setting_uncor_weak_p1000.Rda")
CV_setting_uncor_weak_p5000 <- t(sfSapply(1:1000, function(z) perform.study.CV(seed = z, p = 5000, effect = 1)))
save(CV_setting_uncor_weak_p5000, file = "CV_setting_uncor_weak_p5000.Rda")  

# settings with medium effects
CV_setting_uncor_medium_p200 <- t(sfSapply(1:1000, function(z) perform.study.CV(seed = z, p = 200, effect = 2)))
save(CV_setting_uncor_medium_p200, file = "CV_setting_uncor_medium_p200.Rda")
CV_setting_uncor_medium_p1000 <- t(sfSapply(1:1000, function(z) perform.study.CV(seed = z, p = 1000, effect = 2)))
save(CV_setting_uncor_medium_p1000, file = "CV_setting_uncor_medium_p1000.Rda")
CV_setting_uncor_medium_p5000 <- t(sfSapply(1:1000, function(z) perform.study.CV(seed = z, p = 5000, effect = 2)))
save(CV_setting_uncor_medium_p5000, file = "CV_setting_uncor_medium_p5000.Rda")  


sfStop()
