###############################################################################
#                                                                             #
#   Reproducible files to the paper 'Pitfalls of performing 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>                     #
#                                                                             #
#   File for reproducing Figure 1                                             #
#                                                                             #
###############################################################################

# function to perform the study to obtain test statistics for original samples, bootstrap samples and subsamples if H0 is true in 'real world'

# input parameters: 
#   - niter: number of iterations
#   - n: number of observations
#   - df: degrees of freedom the test should have

# output: (niter x 3) matrix
#         1st column contains values for test statistic in original samples, 
#         2nd column contains values for bootstrap samples, 
#         3rd column contains values for subsamples 

perform.LR.test <- function(niter = 100, n = 1000, df = 1){
  
  teststat.orig <- vector()
  teststat.boot  <- vector()
  teststat.sub   <- vector()
  
  for(i in 1:niter){
  
    set.seed(i)
  
    # generate data
    x    <- matrix(rnorm(n = df * n, 0, 1), nrow = n)
    beta <- rep(0, df) # under H0
    
    data <- data.frame(y = x %*% beta + rnorm(n, 0, 1),
                       x = x)

    # perform test for original sample
    teststat.orig[i] <- obtain.teststat(dat = data)
    
    # perform test for bootstrap sample    
    data.boot <- data[sample(1:n, replace = TRUE, size = n),]
    teststat.boot[i] <- obtain.teststat(dat = data.boot)

    # perform test for subsample (results not shown in the paper)
    data.sub <- data[sample(1:n, replace = FALSE, size = floor(0.632 * n)),]
    teststat.sub[i] <- obtain.teststat(dat = data.sub)
    

  }
  
  teststat <- data.frame(original  = teststat.orig,
                         bootstrap = teststat.boot,
                         subsample = teststat.sub)
}


# function to compute test statistic from LR-test (is called by function perform.study())

obtain.teststat <- function(dat){
  
  model_full    <- lm(y ~ ., data = dat)
  model_reduced <- lm(y ~ 1, data = dat)
  
  teststat <- as.numeric(unlist(anova(model_full, model_reduced)["F"]))[2]
  
}



# compute test statistics from LR-test under null hypothesis (this may take some days)

teststat.LR.test.H0 <- perform.LR.test(niter = 500000, n = 100000, df = 1)  # test with 1 degrees of freedom
save(teststat.LR.test.H0, file = "teststat_LR_test_H0.Rda") # save result in the current working directory



# create Figure 1: plot of densities T and T* (is saved in the current working directory)

graphics.off()
pdf(height = 4, width = 4.5, file = "testatistic_LRtest.pdf")

plot(density(teststat.LR.test.H0[,1], from = 0, adjust = 0.9), type = "l", xlim = c(0, 8), ylim = c(0, 1.3), xlab = "Test statistic", main = "Empirical density functions for LR-test", cex.lab = 1.1, cex.main = 1.2, lwd = 2)
lines(seq(0.000001, 20, length = 100000), dchisq(seq(0.000001, 20, length = 100000), df = 1), lty = 1, col = "red", lwd = 0.8)
lines(density(teststat.LR.test.H0[,2], from = 0, adjust = 0.9), lty = 2, lwd = 2)
legend("topright", col = c("black", "black", "red"), lty = c(1, 2, 1), lwd = c(2, 2, 1), legend = c("T", "T*", expression(paste(chi^2, "(1)", sep = ""))))

graphics.off()
