# Load results:
load("./FAbatchPaper/Results/SimulationResults.Rda")


# Combine results to a data frame:

indscen <- sapply(SimulationResults, function(x) x$indscen)
metric <- sapply(SimulationResults, function(x) x$metric)
batchremmethod <- sapply(SimulationResults, function(x) x$batchremmethod)
iteration <- sapply(SimulationResults, function(x) x$iteration)
metrval <- sapply(SimulationResults, function(x) x$metrval)

results <- data.frame(indscen=indscen, metric=metric,
  batchremmethod=batchremmethod, iteration=iteration, metrval=metrval)

results$batchremmethod <- as.character(results$batchremmethod)
results$batchremmethod[results$batchremmethod=="standardize"] <- "stand"
results$batchremmethod[results$batchremmethod=="meancenter"] <- "meanc"
results$batchremmethod <- factor(results$batchremmethod, levels=c("none", 
  "fabatch", "combat", "sva", "meanc", "stand", "ratiog", "ratioa")) 

results$indscen <- factor(as.character(results$indscen), 
  levels=c("b0scen", "bjscen", "btildescen", "CommonCor", "BatchspecificCor", "BatchClassspecificCor"))

results$metric <- as.character(results$metric)
results$metric[results$metric=="sep"] <- "sepscore"
results$metric[results$metric=="kldist"] <- "klmetr"
results$metric[results$metric=="skew"] <- "skewdiv"
results$metric[results$metric=="cor"] <- "corbeaf"
results$metric <- factor(results$metric, levels=c("sepscore", "avedist", 
  "klmetr", "pvca", "diffexpr", "skewdiv", "corbeaf"))


  
  
  
# Boxplots of metric values:
##############################

library("ggplot2")

rows <- c("Design A - ComCor",
"Design A - BatchCor",
"Design A - BatchClassCor",
"Design B - ComCor",
"Design B - BatchCor",
"Design B - BatchClassCor")

for(metind in 1:length(levels(results$metric))) {

data1 <- subset(results, metric==levels(metric)[metind])
levels(data1$indscen)[levels(data1$indscen)=="b0scen"] <- rows[1]
levels(data1$indscen)[levels(data1$indscen)=="bjscen"] <- rows[2]
levels(data1$indscen)[levels(data1$indscen)=="btildescen"] <- rows[3]
levels(data1$indscen)[levels(data1$indscen)=="CommonCor"] <- rows[4]
levels(data1$indscen)[levels(data1$indscen)=="BatchspecificCor"] <- rows[5]
levels(data1$indscen)[levels(data1$indscen)=="BatchClassspecificCor"] <- rows[6]

data1$indscen <- factor(data1$indscen, levels=levels(data1$indscen)[c(1,4,2,5,3,6)])

p1 <- ggplot(data=data1, aes(x=batchremmethod, y=metrval)) + 
  geom_boxplot() + facet_wrap(~ indscen, nrow=3, ncol=2) + 
  scale_x_discrete(labels=gsub("Standardize", "Stand.", gsub("Meancenter", "Meanc.", levels(data1$batchremmethod)))) +
  labs(x="Method", y=levels(data1$metric)[metind]) +
  theme(axis.title=element_text(size=15), legend.position="none",
  strip.text.x=element_text(size=12), axis.text=element_text(size=12),
  axis.text.x=element_text(angle=45, vjust=1, hjust=1)) +
  theme(axis.title=element_blank()) # + scale_fill_manual(values=c("lightgrey", rep("white", length(levels(data1$batchremmethod))-1)))

p1

ggsave(file=paste("./FAbatchPaper/Results/SupplementaryFigure", metind + 7, ".pdf", sep=""), width=9*(7/8), height=10*(7/8))

}






# Calculate mean performances of the methods for the individual metrics,
# with respect to the original metric values and also with respect to
# the ranks
###########################################################################


# Calculate ranks of the methods with respect to the individual metrics:

results$simdataset <- factor(paste(results$indscen, results$iteration, sep="_"))


resultsranks <- results
isdecreasing <- c(FALSE, FALSE, FALSE, TRUE, TRUE, FALSE, TRUE)

for(i in 1:length(levels(results$metric))) {
  for(j in 1:length(levels(results$simdataset))) {
    if(isdecreasing[i])
      ranktemp <- rank(-results[(results$metric == levels(results$metric)[i]) & (results$simdataset == levels(results$simdataset)[j]),]$metrval)
    else
      ranktemp <- rank(results[(results$metric == levels(results$metric)[i]) & (results$simdataset == levels(results$simdataset)[j]),]$metrval)
    resultsranks[(results$metric == levels(results$metric)[i]) & (results$simdataset == levels(results$simdataset)[j]),]$metrval <- ranktemp
  }
  cat(paste("Iteration ", i, " of ", length(levels(results$metric)), sep=""), "\n")
}


# Because the above step is very time consuming it is best to perform it only once
# and store the results to spare time. Once it has been performed one, the above
# part can be commented out and the code below used:

# save(results, resultsranks, file="./FAbatchPaper/Results/resultswithranktable.Rda")
# load("./FAbatchPaper/Results/resultswithranktable.Rda")



# Take means over datasets: 
  
library("plyr")
resultsum <- ddply(results, .variables=c("indscen", "metric", "batchremmethod"), .fun=summarise, mean.metrval=mean(metrval))


# Take means of the ranks over datasets: 

library("plyr")
resultsumranks <- ddply(resultsranks, .variables=c("indscen", "metric", "batchremmethod"), .fun=summarise, mean.metrval=mean(metrval))



# Generate list containing the means of the metrics for the individual methods
# sorted by mean performance:

isdecreasing <- c(FALSE, FALSE, FALSE, TRUE, TRUE, FALSE, TRUE)
ranklist <- list()

for(metind in 1:length(levels(resultsum$metric))) {

ranklist[[metind]] <- list()

subsetbool <- (resultsum$indscen==levels(resultsum$indscen)[1]) & 
  (resultsum$metric==levels(resultsum$metric)[metind])
ordering <- order(resultsum$mean.metrval[subsetbool], decreasing=isdecreasing[metind])
namesmet <- as.character(resultsum$batchremmethod[subsetbool][ordering])
valuesmet <- resultsum$mean.metrval[subsetbool][ordering]
names(valuesmet) <- namesmet
ranklist[[metind]][[1]] <- valuesmet

for(i in 2:length(unique(resultsum$indscen))) {
  subsetbool <- (resultsum$indscen==levels(resultsum$indscen)[i]) &  
    (resultsum$metric==levels(resultsum$metric)[metind])
  ordering <- order(resultsum$mean.metrval[subsetbool], decreasing=isdecreasing[metind])
  namesmet <- as.character(resultsum$batchremmethod[subsetbool][ordering])
  valuesmet <- resultsum$mean.metrval[subsetbool][ordering]
  names(valuesmet) <- namesmet
  ranklist[[metind]][[i]] <- valuesmet
}

names(ranklist[[metind]]) <- paste(levels(resultsum$indscen), levels(resultsum$metric)[metind], sep="_")

}

names(ranklist) <- levels(resultsum$metric)




# Generate list containing the means of the ranks of the metrics for the individual methods
# sorted by mean performance according to the ranks:

rankranklist <- list()

for(metind in 1:length(levels(resultsumranks$metric))) {

rankranklist[[metind]] <- list()

subsetbool <- (resultsumranks$indscen==levels(resultsumranks$indscen)[1]) & 
  (resultsumranks$metric==levels(resultsumranks$metric)[metind])
ordering <- order(resultsumranks$mean.metrval[subsetbool])
namesmet <- as.character(resultsumranks$batchremmethod[subsetbool][ordering])
valuesmet <- resultsumranks$mean.metrval[subsetbool][ordering]
names(valuesmet) <- namesmet
rankranklist[[metind]][[1]] <- valuesmet

for(i in 2:length(unique(resultsumranks$indscen))) {
  subsetbool <- (resultsumranks$indscen==levels(resultsumranks$indscen)[i]) &  
    (resultsumranks$metric==levels(resultsumranks$metric)[metind])
  ordering <- order(resultsumranks$mean.metrval[subsetbool])
  namesmet <- as.character(resultsumranks$batchremmethod[subsetbool][ordering])
  valuesmet <- resultsumranks$mean.metrval[subsetbool][ordering]
  names(valuesmet) <- namesmet
  rankranklist[[metind]][[i]] <- valuesmet
}

names(rankranklist[[metind]]) <- paste(levels(resultsumranks$indscen), levels(resultsumranks$metric)[metind], sep="_")

}

names(rankranklist) <- levels(resultsumranks$metric)




# Generate LaTeX-Code for the tables in the Supplementary Materials of the paper:

rowslong <- c("Factor induced correlations - Common Correlations",
"Factor induced correlations - Batch-specific Correlations",
"Factor induced correlations - Batch-class-specific Correlations",
"Correlations estimated on real data - Common Correlations",
"Correlations estimated on real data - Batch-specific Correlations",
"Correlations estimated on real data - Batch-class-specific Correlation")


tablesrankranklist <- mapply(function(xmeans, xranks, indmetr) {contents <- unlist(lapply(1:6, function(y) {
    contentxmeans <- round(xmeans[[y]], 5); 
    namescontentxmeans <- names(xmeans[[y]]);
    whichfaxmeans <- which(names(xmeans[[y]])=="fabatch");
    contentxmeans[whichfaxmeans] <- paste("\\textbf{", contentxmeans[whichfaxmeans], "}", sep="");
    namescontentxmeans[whichfaxmeans] <- paste("\\textbf{", namescontentxmeans[whichfaxmeans], "}", sep="");
	contentxranks <- round(xranks[[y]], 5); 
    namescontentxranks <- names(xranks[[y]]);
    whichfaxranks <- which(names(xranks[[y]])=="fabatch");
    contentxranks[whichfaxranks] <- paste("\\textbf{", contentxranks[whichfaxranks], "}", sep="");
    namescontentxranks[whichfaxranks] <- paste("\\textbf{", namescontentxranks[whichfaxranks], "}", sep="");
c(paste("\\multicolumn{9}{c}{", rowslong[y], "}\\\\", sep=""),
"\\hline",
paste("\\multirow{2}{*}{mean values} & ", paste(namescontentxmeans, collapse=" & "), 
"\\\\", sep=""),
paste(" & ", paste(contentxmeans, collapse=" & "), "\\\\", sep=""),
"\\hline\\hline",
paste("\\multirow{2}{*}{mean ranks} & ", paste(namescontentxranks, collapse=" & "), 
"\\\\", sep=""),
paste(" & ", paste(contentxranks, collapse=" & "), "\\\\", sep=""),
"\\hline",
"\\multicolumn{9}{c}{}\\\\[-5.5pt]")}));
contents <- contents[-length(contents)];
c("\\newpage",
"",
paste("\\refstepcounter{", gsub("b0scen_", "", names(xmeans)[1]), "s}", sep=""),
paste("\\label{t:", gsub("b0scen_", "", names(xmeans)[1]), "}", sep=""),
"",
"\\vspace*{2cm}",
"\\hbox{\\hspace{-3cm}\\hvFloat[%",
"onlyText=true,",
"nonFloat=true,",
"objectAngle=90,%",
"capPos=l,%",
"capAngle=90,",
"capWidth=h]{table}{",
"\\small",
"%",
"\\bgroup",
"\\def\\arraystretch{0.75}",
"\\begin{tabular}{| l | c | c | c | c | c | c | c | c |}",
contents, 
"\\end{tabular}",
"\\egroup",
paste("}{\\fontsize{9pt}{9pt}\\selectfont Supplementary Table ", indmetr, ": ", gsub("b0scen_", "", names(xmeans)[1]), "}{lab", indmetr, "}}", sep=""))
}, ranklist, rankranklist, 1:length(ranklist), SIMPLIFY=FALSE)


for(i in 1:length(tablesrankranklist)) {

cat("\n")
cat("\n")
cat("\n")
cat("\n")

for(j in 1:length(tablesrankranklist[[i]])) {
  cat(tablesrankranklist[[i]][j],"\n")
}

}
