# Load results out of the four jobs used to obtain the results
# and combine to corresponding four lists to a single list:


load("./ICSVCodeData/Results/ResultsLargeDataStudy_1.Rda")
res1 <- Results

load("./ICSVCodeData/Results/ResultsLargeDataStudy_2.Rda")
res2 <- Results

load("./ICSVCodeData/Results/ResultsLargeDataStudy_3.Rda")
res3 <- Results

load("./ICSVCodeData/Results/ResultsLargeDataStudy_4.Rda")
res4 <- Results

Results <- c(res1, res2, res3, res4)




# Store results in a data frame:


setting <- sapply(Results, function(x) x$setting)
acctrain <- sapply(Results, function(x) x$acctrain)
acctest <- sapply(Results, function(x) x$acctest)
batchremmethod <- sapply(Results, function(x) x$batchremmethod)

ntrain <- sapply(Results, function(x) x$ntrain)

true1pred1plslda <- sapply(Results, function(x) x$true1pred1plslda)
true1pred2plslda <- sapply(Results, function(x) x$true1pred2plslda)
true2pred1plslda <- sapply(Results, function(x) x$true2pred1plslda)
true2pred2plslda <- sapply(Results, function(x) x$true2pred2plslda)

true1pred1plsldavarsel <- sapply(Results, function(x) x$true1pred1plsldavarsel)
true1pred2plsldavarsel <- sapply(Results, function(x) x$true1pred2plsldavarsel)
true2pred1plsldavarsel <- sapply(Results, function(x) x$true2pred1plsldavarsel)
true2pred2plsldavarsel <- sapply(Results, function(x) x$true2pred2plsldavarsel)

true1pred1boosting <- sapply(Results, function(x) x$true1pred1boosting)
true1pred2boosting <- sapply(Results, function(x) x$true1pred2boosting)
true2pred1boosting <- sapply(Results, function(x) x$true2pred1boosting)
true2pred2boosting <- sapply(Results, function(x) x$true2pred2boosting)

true1pred1boostingvarsel <- sapply(Results, function(x) x$true1pred1boostingvarsel)
true1pred2boostingvarsel <- sapply(Results, function(x) x$true1pred2boostingvarsel)
true2pred1boostingvarsel <- sapply(Results, function(x) x$true2pred1boostingvarsel)
true2pred2boostingvarsel <- sapply(Results, function(x) x$true2pred2boostingvarsel)

true1pred1nsc <- sapply(Results, function(x) x$true1pred1nsc)
true1pred2nsc <- sapply(Results, function(x) x$true1pred2nsc)
true2pred1nsc <- sapply(Results, function(x) x$true2pred1nsc)
true2pred2nsc <- sapply(Results, function(x) x$true2pred2nsc)

true1pred1rf <- sapply(Results, function(x) x$true1pred1rf)
true1pred2rf <- sapply(Results, function(x) x$true1pred2rf)
true2pred1rf <- sapply(Results, function(x) x$true2pred1rf)
true2pred2rf <- sapply(Results, function(x) x$true2pred2rf)

true1pred1knnvarsel <- sapply(Results, function(x) x$true1pred1knnvarsel)
true1pred2knnvarsel <- sapply(Results, function(x) x$true1pred2knnvarsel)
true2pred1knnvarsel <- sapply(Results, function(x) x$true2pred1knnvarsel)
true2pred2knnvarsel <- sapply(Results, function(x) x$true2pred2knnvarsel)

traintestsetsize <- sapply(Results, function(x) x$traintestsetsize)
sizetrain <- rep("large", length(Results))
sizetrain[traintestsetsize %in% c("trainsmalltestlarge", "trainsmalltestsmall",
                                  "trainsmalltestverysmall")] <- "small"
sizetest <- rep("large", length(Results))
sizetest[traintestsetsize %in% c("trainlargetestsmall", "trainsmalltestsmall")] <- "small"
sizetest[traintestsetsize %in% c("trainlargetestverysmall", "trainsmalltestverysmall")] <- "verysmall"

propclass2 <- sapply(Results, function(x) x$propclass2)


resultswide <- data.frame(setting=setting, acctrain=acctrain, acctest=acctest,
                          batchremmethod=batchremmethod, sizetrain=sizetrain, sizetest=sizetest, propclass2=propclass2, ntrain=ntrain, 
                          true1pred1plslda = true1pred1plslda, true1pred2plslda = true1pred2plslda, 
                          true2pred1plslda = true2pred1plslda, true2pred2plslda = true2pred2plslda, 
                          true1pred1plsldavarsel = true1pred1plsldavarsel, true1pred2plsldavarsel = true1pred2plsldavarsel, 
                          true2pred1plsldavarsel = true2pred1plsldavarsel, true2pred2plsldavarsel = true2pred2plsldavarsel,
                          true1pred1boosting = true1pred1boosting, true1pred2boosting = true1pred2boosting, 
                          true2pred1boosting = true2pred1boosting, true2pred2boosting = true2pred2boosting, 
                          true1pred1boostingvarsel = true1pred1boostingvarsel, true1pred2boostingvarsel = true1pred2boostingvarsel, 
                          true2pred1boostingvarsel = true2pred1boostingvarsel, true2pred2boostingvarsel = true2pred2boostingvarsel, 
                          true1pred1nsc = true1pred1nsc, true1pred2nsc = true1pred2nsc, 
                          true2pred1nsc = true2pred1nsc, true2pred2nsc = true2pred2nsc, 
                          true1pred1rf = true1pred1rf, true1pred2rf = true1pred2rf, 
                          true2pred1rf = true2pred1rf, true2pred2rf = true2pred2rf, 
                          true1pred1knnvarsel = true1pred1knnvarsel, true1pred2knnvarsel = true1pred2knnvarsel, 
                          true2pred1knnvarsel = true2pred1knnvarsel, true2pred2knnvarsel = true2pred2knnvarsel)

results <- reshape(resultswide, varying=list(c("true1pred1plslda", "true1pred1plsldavarsel", "true1pred1boosting", "true1pred1boostingvarsel", "true1pred1nsc", "true1pred1rf", "true1pred1knnvarsel"),
                                             c("true1pred2plslda", "true1pred2plsldavarsel", "true1pred2boosting", "true1pred2boostingvarsel", "true1pred2nsc", "true1pred2rf", "true1pred2knnvarsel"),
                                             c("true2pred1plslda", "true2pred1plsldavarsel", "true2pred1boosting", "true2pred1boostingvarsel", "true2pred1nsc", "true2pred1rf", "true2pred1knnvarsel"),
                                             c("true2pred2plslda", "true2pred2plsldavarsel", "true2pred2boosting", "true2pred2boostingvarsel", "true2pred2nsc", "true2pred2rf", "true2pred2knnvarsel")), 
                   v.names=c("true1pred1", "true1pred2", "true2pred1", "true2pred2"), 
                   timevar="classifier", times=c("plslda", "plsldavarsel", "boosting", "boostingvarsel", "nsc", "rf", "knnvarsel"),
                   direction="long")
results$classifier <- as.factor(results$classifier)


results$true1pred1 <- as.numeric(results$true1pred1)
results$true1pred2 <- as.numeric(results$true1pred2)
results$true2pred1 <- as.numeric(results$true2pred1)
results$true2pred2 <- as.numeric(results$true2pred2)


results$classifier <- factor(results$classifier, levels=c("plslda", "plsldavarsel", 
                                                          "boosting", "boostingvarsel",
                                                          "nsc", "rf", "knnvarsel"))







# Separately sum up the TP-, the TN-, the FP- and the FN-values and calculate the 
# MCC_rule-values:


library("plyr")
resultsum <- ddply(results, .variables=c("setting", "batchremmethod", "classifier", 
                                         "sizetrain", "sizetest", "acctrain"), .fun=summarise, true1pred1=sum(true1pred1),
                   true1pred2=sum(true1pred2), true2pred1=sum(true2pred1), true2pred2=sum(true2pred2))

resultsum$true1pred1 <- as.numeric(resultsum$true1pred1)
resultsum$true1pred2 <- as.numeric(resultsum$true1pred2)
resultsum$true2pred1 <- as.numeric(resultsum$true2pred1)
resultsum$true2pred2 <- as.numeric(resultsum$true2pred2)


MCC <- function(tp, tn, fp, fn)
  (tp*tn - fp*fn)/sqrt((tp+fp)*(tp+fn)*(tn+fp)*(tn+fn))

resultsum$mcc <- MCC(resultsum$true2pred2, resultsum$true1pred1,
                     resultsum$true1pred2, resultsum$true2pred1)


resultsum$mcc[is.na(resultsum$mcc)] <- 0 


trainingdatasize <- expand.grid(sizetrain=unique(results$sizetrain), 
    acctrain=unique(results$acctrain), stringsAsFactors=FALSE)[,2:1]
trainingdatasize$ntrain <- apply(trainingdatasize, 1, function(x) 
  results$ntrain[(results$acctrain==x[1]) & (results$sizetrain==x[2])][1])

resultsum$ntrain <- apply(cbind(as.character(resultsum$acctrain), as.character(resultsum$sizetrain)), 1, function(x) 
  trainingdatasize$ntrain[(trainingdatasize$acctrain==x[1]) & (trainingdatasize$sizetrain==x[2])])



# Assign different labels to the setting parameters:

levels(resultsum$sizetrain)
levels(resultsum$sizetest)

levels(resultsum$sizetrain) <- c("trainlarge", "trainsmall")
levels(resultsum$sizetest) <- c("testlarge", "testsmall", "testverysmall")

levels(resultsum$classifier)
levels(resultsum$classifier) <- c("PLS-LDA", "PLS-LDAvarsel", "Boosting", "Boostingvarsel", "NSC", "RF", "kNNvarsel")

levels(resultsum$setting)
levels(resultsum$setting) <- c("Addon norm.", "Sep. norm.")

resultsum$batchremmethod <- as.character(resultsum$batchremmethod)
resultsum$batchremmethod[resultsum$batchremmethod=="meancenter"] <- "meanc"
resultsum$batchremmethod[resultsum$batchremmethod=="standardize"] <- "stand"
resultsum$batchremmethod[resultsum$batchremmethod=="fsvafast"] <- "fsva_f"
resultsum$batchremmethod[resultsum$batchremmethod=="fsvaexact"] <- "fsva_e"
resultsum$batchremmethod <- factor(resultsum$batchremmethod, 
                                   levels=c("none", "combat", "meanc", "stand", "ratiog", "ratioa", "fsva_f", "fsva_e", "fabatch"))










# Figures S1 to S7 in the Supplementary Materials: Boxplots of the MCC_rule-values


library("ggplot2")

classtemp <- "Boosting"

margins <- c(0.5,0.3,0.3,0.3)

p1 <- ggplot(data=resultsum[resultsum$classifier==classtemp & resultsum$sizetest=="testlarge",], aes(x=batchremmethod, y=mcc)) + 
  facet_wrap(~ sizetrain, nrow=1, ncol=2, labeller=labeller(sizetrain=c('trainlarge'="Large training dataset", 'trainsmall'="Small training dataset"))) + geom_boxplot(aes(colour=setting)) + theme_bw() + 
  theme(legend.position="none", axis.text.x=element_text(angle=45, vjust=1, hjust=1, size=rel(1.1), colour="black"), 
        plot.title=element_text(size=rel(1.1)), plot.margin = unit(c(0.3,0.3,0.3,0.3), "cm"), strip.text.x=element_text(size=rel(1.2))) + labs(x="batch effect adjustment method", y="MCC") +
  scale_y_continuous(limits = range(resultsum$mcc[resultsum$classifier==classtemp])) + ylab(expression(MCC["rule"])) +
  ggtitle("Large test dataset")

p2 <- ggplot(data=resultsum[resultsum$classifier==classtemp & resultsum$sizetest=="testsmall",], aes(x=batchremmethod, y=mcc)) + 
  facet_wrap(~ sizetrain, nrow=1, ncol=2, labeller=labeller(sizetrain=c('trainlarge'="Large training dataset", 'trainsmall'="Small training dataset"))) + geom_boxplot(aes(colour=setting)) + theme_bw() + 
  theme(legend.position="none", axis.text.x=element_text(angle=45, vjust=1, hjust=1, size=rel(1.1), colour="black"), plot.title=element_text(size=rel(1.1)), plot.margin = unit(margins, "cm"), strip.text.x=element_text(size=rel(1.2))) + labs(x="batch effect adjustment method", y="MCC") +
  scale_y_continuous(limits = range(resultsum$mcc[resultsum$classifier==classtemp])) + ylab(expression(MCC["rule"])) +
  ggtitle("Small test dataset")

p3 <- ggplot(data=resultsum[resultsum$classifier==classtemp & resultsum$sizetest=="testverysmall",], aes(x=batchremmethod, y=mcc)) + 
  facet_wrap(~ sizetrain, nrow=1, ncol=2, labeller=labeller(sizetrain=c('trainlarge'="Large training dataset", 'trainsmall'="Small training dataset"))) + geom_boxplot(aes(colour=setting)) + theme_bw() + 
  theme(legend.position="none", axis.text.x=element_text(angle=45, vjust=1, hjust=1, size=rel(1.1), colour="black"), plot.title=element_text(size=rel(1.1)), plot.margin = unit(margins, "cm"), strip.text.x=element_text(size=rel(1.2))) + labs(x="batch effect adjustment method", y="MCC") +
  scale_y_continuous(limits = range(resultsum$mcc[resultsum$classifier==classtemp])) + ylab(expression(MCC["rule"])) +
  ggtitle("Very small test dataset")

library("grid")
library("gridExtra")
pall <- grid.arrange(p1, p2, p3, ncol = 1)

grid.draw(pall)

ggsave(file=paste("./ICSVCodeData/Results/", classtemp, ".eps", sep=""), plot=pall, width=7, height=10)





classtemp <- "Boostingvarsel"

margins <- c(0.5,0.3,0.3,0.3)

p1 <- ggplot(data=resultsum[resultsum$classifier==classtemp & resultsum$sizetest=="testlarge",], aes(x=batchremmethod, y=mcc)) + 
  facet_wrap(~ sizetrain, nrow=1, ncol=2, labeller=labeller(sizetrain=c('trainlarge'="Large training dataset", 'trainsmall'="Small training dataset"))) + geom_boxplot(aes(colour=setting)) + theme_bw() + 
  theme(legend.position="none", axis.text.x=element_text(angle=45, vjust=1, hjust=1, size=rel(1.1), colour="black"), 
        plot.title=element_text(size=rel(1.1)), plot.margin = unit(c(0.3,0.3,0.3,0.3), "cm"), strip.text.x=element_text(size=rel(1.2))) + labs(x="batch effect adjustment method", y="MCC") +
  scale_y_continuous(limits = range(resultsum$mcc[resultsum$classifier==classtemp])) + ylab(expression(MCC["rule"])) +
  ggtitle("Large test dataset")

p2 <- ggplot(data=resultsum[resultsum$classifier==classtemp & resultsum$sizetest=="testsmall",], aes(x=batchremmethod, y=mcc)) + 
  facet_wrap(~ sizetrain, nrow=1, ncol=2, labeller=labeller(sizetrain=c('trainlarge'="Large training dataset", 'trainsmall'="Small training dataset"))) + geom_boxplot(aes(colour=setting)) + theme_bw() + 
  theme(legend.position="none", axis.text.x=element_text(angle=45, vjust=1, hjust=1, size=rel(1.1), colour="black"), plot.title=element_text(size=rel(1.1)), plot.margin = unit(margins, "cm"), strip.text.x=element_text(size=rel(1.2))) + labs(x="batch effect adjustment method", y="MCC") +
  scale_y_continuous(limits = range(resultsum$mcc[resultsum$classifier==classtemp])) + ylab(expression(MCC["rule"])) +
  ggtitle("Small test dataset")

p3 <- ggplot(data=resultsum[resultsum$classifier==classtemp & resultsum$sizetest=="testverysmall",], aes(x=batchremmethod, y=mcc)) + 
  facet_wrap(~ sizetrain, nrow=1, ncol=2, labeller=labeller(sizetrain=c('trainlarge'="Large training dataset", 'trainsmall'="Small training dataset"))) + geom_boxplot(aes(colour=setting)) + theme_bw() + 
  theme(legend.position="none", axis.text.x=element_text(angle=45, vjust=1, hjust=1, size=rel(1.1), colour="black"), plot.title=element_text(size=rel(1.1)), plot.margin = unit(margins, "cm"), strip.text.x=element_text(size=rel(1.2))) + labs(x="batch effect adjustment method", y="MCC") +
  scale_y_continuous(limits = range(resultsum$mcc[resultsum$classifier==classtemp])) + ylab(expression(MCC["rule"])) +
  ggtitle("Very small test dataset")

pall <- grid.arrange(p1, p2, p3, ncol = 1)

grid.draw(pall)

ggsave(file=paste("./ICSVCodeData/Results/", classtemp, ".eps", sep=""), plot=pall, width=7, height=10)





classtemp <- "NSC"

margins <- c(0.5,0.3,0.3,0.3)

p1 <- ggplot(data=resultsum[resultsum$classifier==classtemp & resultsum$sizetest=="testlarge",], aes(x=batchremmethod, y=mcc)) + 
  facet_wrap(~ sizetrain, nrow=1, ncol=2, labeller=labeller(sizetrain=c('trainlarge'="Large training dataset", 'trainsmall'="Small training dataset"))) + geom_boxplot(aes(colour=setting)) + theme_bw() + 
  theme(legend.position="none", axis.text.x=element_text(angle=45, vjust=1, hjust=1, size=rel(1.1), colour="black"), 
        plot.title=element_text(size=rel(1.1)), plot.margin = unit(c(0.3,0.3,0.3,0.3), "cm"), strip.text.x=element_text(size=rel(1.2))) + labs(x="batch effect adjustment method", y="MCC") +
  scale_y_continuous(limits = range(resultsum$mcc[resultsum$classifier==classtemp])) + ylab(expression(MCC["rule"])) +
  ggtitle("Large test dataset")

p2 <- ggplot(data=resultsum[resultsum$classifier==classtemp & resultsum$sizetest=="testsmall",], aes(x=batchremmethod, y=mcc)) + 
  facet_wrap(~ sizetrain, nrow=1, ncol=2, labeller=labeller(sizetrain=c('trainlarge'="Large training dataset", 'trainsmall'="Small training dataset"))) + geom_boxplot(aes(colour=setting)) + theme_bw() + 
  theme(legend.position="none", axis.text.x=element_text(angle=45, vjust=1, hjust=1, size=rel(1.1), colour="black"), plot.title=element_text(size=rel(1.1)), plot.margin = unit(margins, "cm"), strip.text.x=element_text(size=rel(1.2))) + labs(x="batch effect adjustment method", y="MCC") +
  scale_y_continuous(limits = range(resultsum$mcc[resultsum$classifier==classtemp])) + ylab(expression(MCC["rule"])) +
  ggtitle("Small test dataset")

p3 <- ggplot(data=resultsum[resultsum$classifier==classtemp & resultsum$sizetest=="testverysmall",], aes(x=batchremmethod, y=mcc)) + 
  facet_wrap(~ sizetrain, nrow=1, ncol=2, labeller=labeller(sizetrain=c('trainlarge'="Large training dataset", 'trainsmall'="Small training dataset"))) + geom_boxplot(aes(colour=setting)) + theme_bw() + 
  theme(legend.position="none", axis.text.x=element_text(angle=45, vjust=1, hjust=1, size=rel(1.1), colour="black"), plot.title=element_text(size=rel(1.1)), plot.margin = unit(margins, "cm"), strip.text.x=element_text(size=rel(1.2))) + labs(x="batch effect adjustment method", y="MCC") +
  scale_y_continuous(limits = range(resultsum$mcc[resultsum$classifier==classtemp])) + ylab(expression(MCC["rule"])) +
  ggtitle("Very small test dataset")

pall <- grid.arrange(p1, p2, p3, ncol = 1)

grid.draw(pall)

ggsave(file=paste("./ICSVCodeData/Results/", classtemp, ".eps", sep=""), plot=pall, width=7, height=10)





classtemp <- "PLS-LDA"

margins <- c(0.5,0.3,0.3,0.3)

p1 <- ggplot(data=resultsum[resultsum$classifier==classtemp & resultsum$sizetest=="testlarge",], aes(x=batchremmethod, y=mcc)) + 
  facet_wrap(~ sizetrain, nrow=1, ncol=2, labeller=labeller(sizetrain=c('trainlarge'="Large training dataset", 'trainsmall'="Small training dataset"))) + geom_boxplot(aes(colour=setting)) + theme_bw() + 
  theme(legend.position="none", axis.text.x=element_text(angle=45, vjust=1, hjust=1, size=rel(1.1), colour="black"), 
        plot.title=element_text(size=rel(1.1)), plot.margin = unit(c(0.3,0.3,0.3,0.3), "cm"), strip.text.x=element_text(size=rel(1.2))) + labs(x="batch effect adjustment method", y="MCC") +
  scale_y_continuous(limits = range(resultsum$mcc[resultsum$classifier==classtemp])) + ylab(expression(MCC["rule"])) +
  ggtitle("Large test dataset")

p2 <- ggplot(data=resultsum[resultsum$classifier==classtemp & resultsum$sizetest=="testsmall",], aes(x=batchremmethod, y=mcc)) + 
  facet_wrap(~ sizetrain, nrow=1, ncol=2, labeller=labeller(sizetrain=c('trainlarge'="Large training dataset", 'trainsmall'="Small training dataset"))) + geom_boxplot(aes(colour=setting)) + theme_bw() + 
  theme(legend.position="none", axis.text.x=element_text(angle=45, vjust=1, hjust=1, size=rel(1.1), colour="black"), plot.title=element_text(size=rel(1.1)), plot.margin = unit(margins, "cm"), strip.text.x=element_text(size=rel(1.2))) + labs(x="batch effect adjustment method", y="MCC") +
  scale_y_continuous(limits = range(resultsum$mcc[resultsum$classifier==classtemp])) + ylab(expression(MCC["rule"])) +
  ggtitle("Small test dataset")

p3 <- ggplot(data=resultsum[resultsum$classifier==classtemp & resultsum$sizetest=="testverysmall",], aes(x=batchremmethod, y=mcc)) + 
  facet_wrap(~ sizetrain, nrow=1, ncol=2, labeller=labeller(sizetrain=c('trainlarge'="Large training dataset", 'trainsmall'="Small training dataset"))) + geom_boxplot(aes(colour=setting)) + theme_bw() + 
  theme(legend.position="none", axis.text.x=element_text(angle=45, vjust=1, hjust=1, size=rel(1.1), colour="black"), plot.title=element_text(size=rel(1.1)), plot.margin = unit(margins, "cm"), strip.text.x=element_text(size=rel(1.2))) + labs(x="batch effect adjustment method", y="MCC") +
  scale_y_continuous(limits = range(resultsum$mcc[resultsum$classifier==classtemp])) + ylab(expression(MCC["rule"])) +
  ggtitle("Very small test dataset")

pall <- grid.arrange(p1, p2, p3, ncol = 1)

grid.draw(pall)

ggsave(file=paste("./ICSVCodeData/Results/", classtemp, ".eps", sep=""), plot=pall, width=7, height=10)





classtemp <- "PLS-LDAvarsel"

margins <- c(0.5,0.3,0.3,0.3)

p1 <- ggplot(data=resultsum[resultsum$classifier==classtemp & resultsum$sizetest=="testlarge",], aes(x=batchremmethod, y=mcc)) + 
  facet_wrap(~ sizetrain, nrow=1, ncol=2, labeller=labeller(sizetrain=c('trainlarge'="Large training dataset", 'trainsmall'="Small training dataset"))) + geom_boxplot(aes(colour=setting)) + theme_bw() + 
  theme(legend.position="none", axis.text.x=element_text(angle=45, vjust=1, hjust=1, size=rel(1.1), colour="black"), 
        plot.title=element_text(size=rel(1.1)), plot.margin = unit(c(0.3,0.3,0.3,0.3), "cm"), strip.text.x=element_text(size=rel(1.2))) + labs(x="batch effect adjustment method", y="MCC") +
  scale_y_continuous(limits = range(resultsum$mcc[resultsum$classifier==classtemp])) + ylab(expression(MCC["rule"])) +
  ggtitle("Large test dataset")

p2 <- ggplot(data=resultsum[resultsum$classifier==classtemp & resultsum$sizetest=="testsmall",], aes(x=batchremmethod, y=mcc)) + 
  facet_wrap(~ sizetrain, nrow=1, ncol=2, labeller=labeller(sizetrain=c('trainlarge'="Large training dataset", 'trainsmall'="Small training dataset"))) + geom_boxplot(aes(colour=setting)) + theme_bw() + 
  theme(legend.position="none", axis.text.x=element_text(angle=45, vjust=1, hjust=1, size=rel(1.1), colour="black"), plot.title=element_text(size=rel(1.1)), plot.margin = unit(margins, "cm"), strip.text.x=element_text(size=rel(1.2))) + labs(x="batch effect adjustment method", y="MCC") +
  scale_y_continuous(limits = range(resultsum$mcc[resultsum$classifier==classtemp])) + ylab(expression(MCC["rule"])) +
  ggtitle("Small test dataset")

p3 <- ggplot(data=resultsum[resultsum$classifier==classtemp & resultsum$sizetest=="testverysmall",], aes(x=batchremmethod, y=mcc)) + 
  facet_wrap(~ sizetrain, nrow=1, ncol=2, labeller=labeller(sizetrain=c('trainlarge'="Large training dataset", 'trainsmall'="Small training dataset"))) + geom_boxplot(aes(colour=setting)) + theme_bw() + 
  theme(legend.position="none", axis.text.x=element_text(angle=45, vjust=1, hjust=1, size=rel(1.1), colour="black"), plot.title=element_text(size=rel(1.1)), plot.margin = unit(margins, "cm"), strip.text.x=element_text(size=rel(1.2))) + labs(x="batch effect adjustment method", y="MCC") +
  scale_y_continuous(limits = range(resultsum$mcc[resultsum$classifier==classtemp])) + ylab(expression(MCC["rule"])) +
  ggtitle("Very small test dataset")

pall <- grid.arrange(p1, p2, p3, ncol = 1)

grid.draw(pall)

ggsave(file=paste("./ICSVCodeData/Results/", classtemp, ".eps", sep=""), plot=pall, width=7, height=10)





classtemp <- "RF"

margins <- c(0.5,0.3,0.3,0.3)

p1 <- ggplot(data=resultsum[resultsum$classifier==classtemp & resultsum$sizetest=="testlarge",], aes(x=batchremmethod, y=mcc)) + 
  facet_wrap(~ sizetrain, nrow=1, ncol=2, labeller=labeller(sizetrain=c('trainlarge'="Large training dataset", 'trainsmall'="Small training dataset"))) + geom_boxplot(aes(colour=setting)) + theme_bw() + 
  theme(legend.position="none", axis.text.x=element_text(angle=45, vjust=1, hjust=1, size=rel(1.1), colour="black"), 
        plot.title=element_text(size=rel(1.1)), plot.margin = unit(c(0.3,0.3,0.3,0.3), "cm"), strip.text.x=element_text(size=rel(1.2))) + labs(x="batch effect adjustment method", y="MCC") +
  scale_y_continuous(limits = range(resultsum$mcc[resultsum$classifier==classtemp])) + ylab(expression(MCC["rule"])) +
  ggtitle("Large test dataset")

p2 <- ggplot(data=resultsum[resultsum$classifier==classtemp & resultsum$sizetest=="testsmall",], aes(x=batchremmethod, y=mcc)) + 
  facet_wrap(~ sizetrain, nrow=1, ncol=2, labeller=labeller(sizetrain=c('trainlarge'="Large training dataset", 'trainsmall'="Small training dataset"))) + geom_boxplot(aes(colour=setting)) + theme_bw() + 
  theme(legend.position="none", axis.text.x=element_text(angle=45, vjust=1, hjust=1, size=rel(1.1), colour="black"), plot.title=element_text(size=rel(1.1)), plot.margin = unit(margins, "cm"), strip.text.x=element_text(size=rel(1.2))) + labs(x="batch effect adjustment method", y="MCC") +
  scale_y_continuous(limits = range(resultsum$mcc[resultsum$classifier==classtemp])) + ylab(expression(MCC["rule"])) +
  ggtitle("Small test dataset")

p3 <- ggplot(data=resultsum[resultsum$classifier==classtemp & resultsum$sizetest=="testverysmall",], aes(x=batchremmethod, y=mcc)) + 
  facet_wrap(~ sizetrain, nrow=1, ncol=2, labeller=labeller(sizetrain=c('trainlarge'="Large training dataset", 'trainsmall'="Small training dataset"))) + geom_boxplot(aes(colour=setting)) + theme_bw() + 
  theme(legend.position="none", axis.text.x=element_text(angle=45, vjust=1, hjust=1, size=rel(1.1), colour="black"), plot.title=element_text(size=rel(1.1)), plot.margin = unit(margins, "cm"), strip.text.x=element_text(size=rel(1.2))) + labs(x="batch effect adjustment method", y="MCC") +
  scale_y_continuous(limits = range(resultsum$mcc[resultsum$classifier==classtemp])) + ylab(expression(MCC["rule"])) +
  ggtitle("Very small test dataset")

pall <- grid.arrange(p1, p2, p3, ncol = 1)

grid.draw(pall)

ggsave(file=paste("./ICSVCodeData/Results/", classtemp, ".eps", sep=""), plot=pall, width=7, height=10)




classtemp <- "kNNvarsel"

margins <- c(0.5,0.3,0.3,0.3)

p1 <- ggplot(data=resultsum[resultsum$classifier==classtemp & resultsum$sizetest=="testlarge",], aes(x=batchremmethod, y=mcc)) + 
  facet_wrap(~ sizetrain, nrow=1, ncol=2, labeller=labeller(sizetrain=c('trainlarge'="Large training dataset", 'trainsmall'="Small training dataset"))) + geom_boxplot(aes(colour=setting)) + theme_bw() + 
  theme(legend.position="none", axis.text.x=element_text(angle=45, vjust=1, hjust=1, size=rel(1.1), colour="black"), 
        plot.title=element_text(size=rel(1.1)), plot.margin = unit(c(0.3,0.3,0.3,0.3), "cm"), strip.text.x=element_text(size=rel(1.2))) + labs(x="batch effect adjustment method", y="MCC") +
  scale_y_continuous(limits = range(resultsum$mcc[resultsum$classifier==classtemp])) + ylab(expression(MCC["rule"])) +
  ggtitle("Large test dataset")

p2 <- ggplot(data=resultsum[resultsum$classifier==classtemp & resultsum$sizetest=="testsmall",], aes(x=batchremmethod, y=mcc)) + 
  facet_wrap(~ sizetrain, nrow=1, ncol=2, labeller=labeller(sizetrain=c('trainlarge'="Large training dataset", 'trainsmall'="Small training dataset"))) + geom_boxplot(aes(colour=setting)) + theme_bw() + 
  theme(legend.position="none", axis.text.x=element_text(angle=45, vjust=1, hjust=1, size=rel(1.1), colour="black"), plot.title=element_text(size=rel(1.1)), plot.margin = unit(margins, "cm"), strip.text.x=element_text(size=rel(1.2))) + labs(x="batch effect adjustment method", y="MCC") +
  scale_y_continuous(limits = range(resultsum$mcc[resultsum$classifier==classtemp])) + ylab(expression(MCC["rule"])) +
  ggtitle("Small test dataset")

p3 <- ggplot(data=resultsum[resultsum$classifier==classtemp & resultsum$sizetest=="testverysmall",], aes(x=batchremmethod, y=mcc)) + 
  facet_wrap(~ sizetrain, nrow=1, ncol=2, labeller=labeller(sizetrain=c('trainlarge'="Large training dataset", 'trainsmall'="Small training dataset"))) + geom_boxplot(aes(colour=setting)) + theme_bw() + 
  theme(legend.position="none", axis.text.x=element_text(angle=45, vjust=1, hjust=1, size=rel(1.1), colour="black"), plot.title=element_text(size=rel(1.1)), plot.margin = unit(margins, "cm"), strip.text.x=element_text(size=rel(1.2))) + labs(x="batch effect adjustment method", y="MCC") +
  scale_y_continuous(limits = range(resultsum$mcc[resultsum$classifier==classtemp])) + ylab(expression(MCC["rule"])) +
  ggtitle("Very small test dataset")

pall <- grid.arrange(p1, p2, p3, ncol = 1)

grid.draw(pall)

ggsave(file=paste("./ICSVCodeData/Results/", classtemp, ".eps", sep=""), plot=pall, width=7, height=10)





# Figure 1: MCC_rule-values for the 34 datasets for each setting without addon batch effect adjustment.

p <- ggplot(data=resultsum[resultsum$batchremmethod=="none",], aes(x=classifier, y=mcc)) + 
  facet_wrap(~ sizetrain + sizetest, nrow=2, ncol=3) + geom_boxplot(aes(fill=setting)) + theme_bw() + 
  theme(legend.position="none", axis.text.x=element_text(angle=45, vjust=1, hjust=1, colour="black"))  + 
  labs(x="classification method", y=expression(MCC["rule"])) + scale_fill_grey(start = 0.8, end = 1)

p

ggsave(file="./ICSVCodeData/Results/Figure1.eps", width=10, height=7)





# Figure 2: Medians of the MCC_rule-values over the 34 datasets for each setting 
# with large training dataset and separate normalization.

resultsmedians <- ddply(resultsum, .variables=c("setting", "batchremmethod", "classifier", "sizetrain", "sizetest"), 
                     .fun=summarise, mcc=median(mcc))
resultsmedianssub <- resultsmedians[(resultsmedians$setting == "Sep. norm.") & (resultsmedians$sizetrain == "trainlarge"),]

p <- ggplot(data=resultsmedianssub, aes(x=batchremmethod, y=mcc, colour=sizetest, group=sizetest)) + 
  facet_wrap(~ classifier, nrow=4, ncol=2) + geom_point(aes(shape=sizetest)) + 
  geom_line(aes(linetype=sizetest)) + 
  scale_colour_grey(start = 0, end = 0.7)+ theme_bw() +
  theme(legend.position="none", axis.text.x=element_text(angle=45, vjust=1, hjust=1, size=rel(1.05), colour="black"), strip.text.x=element_text(size=rel(1.1))) +
  labs(x="batch effect adjustment method", y=expression(paste("Median of ", MCC["rule"], "-values", sep="")))
p  

ggsave(file="./ICSVCodeData/Results/Figure2.eps", width=6, height=8.5)





# In the Results section we state that in the case of Boosting two to three datasets worked
# very bad as training and test datasets, while the other datasets exhibihit an almost 
# perfect performance.
# To investigate this we plot the MCC-values calculated for each training test set pair
# for all settings without batch effect adjustment, Boosting as classifier and 
# addon normalization as normalization procedure.

results$mcc <- MCC(results$true2pred2, results$true1pred1,
                   results$true1pred2, results$true2pred1)

results2 <- results[results$batchremmethod=="none" & results$classifier=="boosting" & results$setting=="addonnorm",]

results2$batchremmethod <- NULL
results2$classifier <- NULL
results2$setting <- NULL

levels(results2$sizetrain)
levels(results2$sizetrain) <- c("trainlarge", "trainsmall")

levels(results2$sizetest)
levels(results2$sizetest) <- c("testlarge", "testsmall", "testverysmall")

p <- ggplot(data=results2, aes(x=acctrain, y=mcc)) + facet_wrap(~sizetrain + sizetest, nrow=3, ncol=2) + geom_boxplot() + geom_line(aes(group=acctest))

p

# --> For two to three test datasets, the MCC-values are consistently and considerably smaller than for the other
#     test datasets. For the latter, the MCC-values are close to one.
#     Similarly, for two training datasets, the MCC-values are considerably smaller than for the other training
#     datasets, where for the latter they are again equal to one.













# Again in the Results section we state that the small 25%-quartiles
# of the MCC_rule-values in the case of Random Forest after batch effect
# adjustment can be attributed to the results of five datasets with
# considerably smaller MCC-values after batch effect adjustment.
# In the following we identify these five datasets:


# Consider only ComBat as batch effect adjustment method:
resultsum2 <- resultsum[resultsum$classifier=="RF" & 
                          resultsum$sizetrain=="trainlarge" & 
                          resultsum$sizetest=="testlarge" &
                          resultsum$setting=="Sep. norm." &
                          resultsum$batchremmethod %in% c("none", "combat"),]
resultsum2$classifier <- NULL
resultsum2$sizetrain <- NULL
resultsum2$sizetest <- NULL
resultsum2$setting <- NULL

# Accession numbers of datasets:
acctrains <- as.character(unique(resultsum2$acctrain))

# Calculate difference between MCC-value before and after batch effect adjustment 
# for each dataset:
mccdiff <- sapply(acctrains, function(x) resultsum2$mcc[resultsum2$batchremmethod=="none" & resultsum2$acctrain==x] -
                    resultsum2$mcc[resultsum2$batchremmethod=="combat" & resultsum2$acctrain==x])
names(mccdiff) <- acctrains

sort(mccdiff)

# --> Identify problematic datasets, where this difference is larger than 0.4:
acctrainsproblem <- acctrains[mccdiff > 0.4]
acctrainsok <- setdiff(unique(resultsum2$acctrain), acctrainsproblem)

# Add variable to dataframe, which indicates whether or not the corresponding training dataset
# belongs to the problematic datasets:
resultsum2$bigdiff <- ifelse(resultsum2$acctrain %in% acctrainsproblem, "Yes", "No")

# Boxplot of MCC_rule-value before and after batch effect adjustment highlighting the
# problematic datasets:
p <- ggplot(data=resultsum2, aes(x=batchremmethod, y=mcc)) + geom_boxplot() + geom_line(aes(group=acctrain, colour=bigdiff)) + geom_point(aes(colour=bigdiff))

p

# --> The small 25%-quartiles of the boxplots after batch effect adjustment
#     can be attributed to these five datsets, where the MCC-values are strongly
#     decreased after batch effect adjustment.



# Investigate properties of these datasets which could be responsible for their bad performance when
# used as training datasets with addon batch effect adjustment:


# Load problematic and unproblematic datasets:

# NOTE: This requires that the raw datasets have been downloaded (by using the script 
# "DownloadData.R" in the subfolder "DataRetrieval") and normalized (by using the script 
# "NormalizationTrainingScript.R" in the subfolder "JobScripts"):

files <- list.files("./ICSVCodeData/Datasets/DatanormfilesLarge", full.names=TRUE)
accnums <- gsub("_datanorm.Rda", "", gsub("./ICSVCodeData/Datasets/DatanormfilesLarge/", "", files))

Xbad <- ybad <- list()

for(i in seq(along=acctrainsproblem)) {
  
  load(files[accnums==acctrainsproblem[i]])
  
  Xbad[[i]] <- X
  ybad[[i]] <- y
  
}

names(Xbad) <- names(ybad) <- acctrainsproblem


Xok <- yok <- list()

for(i in seq(along=acctrainsok)) {
  
  load(files[accnums==acctrainsok[i]])
  
  Xok[[i]] <- X
  yok[[i]] <- y
  
}

names(Xok) <- names(yok) <- acctrainsok



# The five datasets which perform bad when used as trainng datasets in combination with ComBat
# feature a significantly stronger imbalance:

propsmallerclassbad <- sapply(ybad, function(x) min(c(mean(x==1), mean(x==2))))
propsmallerclassok <- sapply(yok, function(x) min(c(mean(x==1), mean(x==2))))

p <- ggplot(data=data.frame(props=c(propsmallerclassbad, propsmallerclassok), 
                            baddataset=c(rep("yes", length(propsmallerclassbad)), 
                                         rep("no", length(propsmallerclassok)))), aes(x=baddataset, y=props)) + 
  geom_boxplot() +  geom_jitter(position=position_jitter(w=.2, h=0))
p

# Wilcoxon test (reported in paper):

wilcox.test(propsmallerclassbad, propsmallerclassok)

# --> Significant difference in proportions of smaller class between problematic
#     and unproblematic datasets.



# In the Discussion section we state that, as expected, for the five problematic
# datasets ComBat lead to classifying almost all test observations as the class
# overrepresented in the training data. This is verified in the following:

resultsum2[resultsum2$acctrain %in% acctrainsproblem,]

unlist(lapply(ybad, function(x) names(table(x))[which.max(table(x))]))




# In the Results section, we state that after excluding the results corresponding
# to training datasets with imbalanced class frequencies, batch effect adjustment
# by combat, meanc, stand, ratiog and ratioa lead to a quite strong improvement 
# of the prediction accuracy of Random Forest. This is verified in the following:

resultsum2 <- resultsum[resultsum$classifier=="RF",]

propall <- c(propsmallerclassbad, propsmallerclassok)
resultsum2$propsmall <- sapply(resultsum2$acctrain, function(x) propall[names(propall)==x])


p <- ggplot(data=resultsum2[resultsum2$propsmall > 0.35,], aes(x=batchremmethod, y=mcc)) + 
  facet_wrap(~ sizetest + sizetrain, nrow=3, ncol=2) + geom_boxplot(aes(colour=setting))

p

# --> Except in the settings with very small test datasets and with addon normalization -
#     where addon batch effect adjustment is not recommended in general - we do see
#     quite a strong improvement.

