## processing figures for simulations
require(reshape2)

load("C:/Users/OllyC/Desktop/2019-DiffLoc/diffloc_sim_wild.rda") # cluster specific bootstrap


auc <- matrix(unlist(lapply(c(sim_results$ROCresBayes,
                              sim_results$ROCresDir, 
                              sim_results$ROCresmr), function(x)x$auc)), ncol = 3, byrow = FALSE)
colnames(auc) <- c("PG", "Dir", "MR")
boxplot(auc[,c(3,2,1)], ylim = c(0.5,1), col = c("blue", "orange", "darkgreen"))

auc_long <- as.data.frame(melt(auc[,c(3,2,1)]))
colnames(auc_long) <- c("Iteration", "Method", "AUC")

gg <- ggplot(data = auc_long, aes(x = Method, y = AUC, fill = Method)) + geom_boxplot(outlier.colour="black", outlier.shape=16,
                                                                                      outlier.size=2, notch=FALSE)
gg <- gg + theme_classic() + scale_fill_manual(values=c("blue", "orange", "darkgreen", "darkorchid4"))
gg <- gg + theme(legend.position="none") + scale_y_continuous(limits = c(0.5,1))
gg <- gg + ggtitle("Cluster specific noise distribution") + xlab("Method") + ylab("AUC")
gg <- gg + scale_x_discrete(labels = c("MR", "Dirichlet", "Polya-Gamma"))
gg


plot(ROC_bayes, type = "b", col = "darkgreen", lwd = 1.5, pch = 19, asp = NA, legacy.axes = T, 
     xlab = "False Positive Rate",
     ylab = "True Positive Rate",
     main = "ROC curve demonstrating proposed methods outperform current approach")
lines(ROC_mr, type = "b", col = "blue", lwd = 1.5, pch = 19, asp = NA)
lines(ROC_dir,type = "b", col = "orange", lwd = 1.5, pch = 19, asp = NA)
legend("bottomright", c("MR", "Dirichlet", "Polya-Gamma"), col = c("blue", "orange", "darkgreen"), pch = 19)

dirnumtrue <- dirnumfalse <- pgnumtrue <- pgnumfalse <- mrnumtrue <- mrnumfalse <- vector(length = 10)


for (l in seq.int(10)){
  
  probdir <- sim_results$probdirres[[l]]
  probpg <- sim_results$probpgres[[l]]
  mrres <- sim_results$mrreslist[[l]]
  perm1_names <- sim_results$perm1_nameslist[[l]]
  
  dirnumtrue[l] <- sum(names(probdir[order(probdir)])[(1 - probdir[order(probdir)]/ncol(resdyndir$allocOut[[1]])) > 0.999] %in% perm1_names)
  dirnumfalse[l] <- sum(!names(probdir[order(probdir)])[(1 - probdir[order(probdir)]/ncol(resdyndir$allocOut[[1]])) > 0.999] %in% perm1_names)
  pgnumtrue[l] <- sum(names(probpg[order(probpg)])[(1 - probpg[order(probpg)]/ncol(resdynpg$allocOut[[1]])) > 0.999] %in% perm1_names)
  pgnumfalse[l] <- sum(!names(probpg[order(probpg)])[(1 - probpg[order(probpg)]/ncol(resdynpg$allocOut[[1]])) > 0.999] %in% perm1_names)
  mrnumtrue[l] <- sum(names(mrres$Mscore[(mrres$Mscore > 2) & mrres$Rscore > 0.9 ] ) %in% perm1_names)
  mrnumfalse[l] <- sum(!names(mrres$Mscore[(mrres$Mscore > 2) & mrres$Rscore > 0.9 ] ) %in% perm1_names)
}

numtrue <- cbind(dirnumtrue, pgnumtrue, mrnumtrue, pg1numtrue)

true_long <- as.data.frame(melt(numtrue[,c(3,1,2)]))
colnames(true_long) <- c("Iteration", "Method", "Relocalisations")

gg2 <- ggplot(data = true_long, aes(x = Method, y = Relocalisations, fill = Method)) + geom_boxplot(outlier.colour="black", outlier.shape=16,
                                                                                                    outlier.size=2, notch=FALSE)
gg2 <- gg2 + theme_classic() + scale_fill_manual(values=c("blue", "orange", "darkgreen", "darkorchid4"))
gg2 <- gg2 + theme(legend.position="none")
gg2 <- gg2 + ggtitle("Cluster specific noise distribution") + xlab("Method") + ylab("Number of Relocalisations")
gg2 <- gg2 + scale_x_discrete(labels = c("MR", "Dirichlet", "Polya-Gamma","Polya-Gamma unshrunk"))
gg2




require(patchwork)

gg / gg2

t.test(auc[,1], auc[,2])
t.test(auc[,1], auc[,3])
t.test(auc[,2], auc[,3])

t.test(numtrue[,1], numtrue[,3])
t.test(numtrue[,1], numtrue[,2])
t.test(numtrue[,2], numtrue[,3])

load("C:/Users/OllyC/Desktop/2019-DiffLoc/diffloc_sim_wild_batch_rand.rda") # random batch effect

auc <- matrix(unlist(lapply(c(sim_results$ROCresBayes,
                              sim_results$ROCresDir, 
                              sim_results$ROCresmr), function(x)x$auc)), ncol = 3, byrow = FALSE)
colnames(auc) <- c("PG", "Dir", "MR")
boxplot(auc[,c(3,2,1)], ylim = c(0.5,1), col = c("blue", "orange", "darkgreen"))

auc_long <- as.data.frame(melt(auc[,c(3,2,1)]))
colnames(auc_long) <- c("Iteration", "Method", "AUC")

gg1 <- ggplot(data = auc_long, aes(x = Method, y = AUC, fill = Method)) + geom_boxplot(outlier.colour="black", outlier.shape=16,
                                                                                      outlier.size=2, notch=FALSE)
gg1 <- gg1 + theme_classic() + scale_fill_manual(values=c("blue", "orange", "darkgreen", "darkorchid4"))
gg1 <- gg1 + theme(legend.position="none") + scale_y_continuous(limits = c(0.3,1))
gg1 <- gg1 + ggtitle("Random batch effects") + xlab("Method") + ylab("AUC")
gg1 <- gg1 + scale_x_discrete(labels = c("MR", "Dirichlet", "Polya-Gamma"))
gg1


plot(ROC_bayes, type = "b", col = "darkgreen", lwd = 1.5, pch = 19, asp = NA, legacy.axes = T, 
     xlab = "False Positive Rate",
     ylab = "True Positive Rate",
     main = "ROC curve demonstrating proposed methods outperform current approach")
lines(ROC_mr, type = "b", col = "blue", lwd = 1.5, pch = 19, asp = NA)
lines(ROC_dir,type = "b", col = "orange", lwd = 1.5, pch = 19, asp = NA)
legend("bottomright", c("MR", "Dirichlet", "Polya-Gamma"), col = c("blue", "orange", "darkgreen"), pch = 19)

dirnumtrue <- dirnumfalse <- pgnumtrue <- pgnumfalse <- mrnumtrue <- mrnumfalse <- vector(length = 10)


for (l in seq.int(10)){
  
  probdir <- sim_results$probdirres[[l]]
  probpg <- sim_results$probpgres[[l]]
  mrres <- sim_results$mrreslist[[l]]
  perm1_names <- sim_results$perm1_nameslist[[l]]
  
  dirnumtrue[l] <- sum(names(probdir[order(probdir)])[(1 - probdir[order(probdir)]/ncol(resdyndir$allocOut[[1]])) > 0.999] %in% perm1_names)
  dirnumfalse[l] <- sum(!names(probdir[order(probdir)])[(1 - probdir[order(probdir)]/ncol(resdyndir$allocOut[[1]])) > 0.999] %in% perm1_names)
  pgnumtrue[l] <- sum(names(probpg[order(probpg)])[(1 - probpg[order(probpg)]/ncol(resdynpg$allocOut[[1]])) > 0.999] %in% perm1_names)
  pgnumfalse[l] <- sum(!names(probpg[order(probpg)])[(1 - probpg[order(probpg)]/ncol(resdynpg$allocOut[[1]])) > 0.999] %in% perm1_names)
  mrnumtrue[l] <- sum(names(mrres$Mscore[(mrres$Mscore > 2) & mrres$Rscore > 0.9 ] ) %in% perm1_names)
  mrnumfalse[l] <- sum(!names(mrres$Mscore[(mrres$Mscore > 2) & mrres$Rscore > 0.9 ] ) %in% perm1_names)
}

numtrue <- cbind(dirnumtrue, pgnumtrue, mrnumtrue, pg1numtrue)

true_long <- as.data.frame(melt(numtrue[,c(3,1,2)]))
colnames(true_long) <- c("Iteration", "Method", "Relocalisations")

gg3 <- ggplot(data = true_long, aes(x = Method, y = Relocalisations, fill = Method)) + geom_boxplot(outlier.colour="black", outlier.shape=16,
                                                                                                    outlier.size=2, notch=FALSE)
gg3 <- gg3 + theme_classic() + scale_fill_manual(values=c("blue", "orange", "darkgreen", "darkorchid4"))
gg3 <- gg3 + theme(legend.position="none")
gg3 <- gg3 + ggtitle("Random batch effects") + xlab("Method") + ylab("Number of Relocalisations")
gg3 <- gg3 + scale_x_discrete(labels = c("MR", "Dirichlet", "Polya-Gamma","Polya-Gamma unshrunk"))
gg3

gg1/gg3

t.test(auc[,1], auc[,2])
t.test(auc[,1], auc[,3])
t.test(auc[,2], auc[,3])

t.test(numtrue[,1], numtrue[,3])
t.test(numtrue[,1], numtrue[,2])
t.test(numtrue[,2], numtrue[,3])


load("C:/Users/OllyC/Desktop/2019-DiffLoc/diffloc_sim_wild_batch_sys.rda") # systematic batch effect

auc <- matrix(unlist(lapply(c(sim_results$ROCresBayes,
                              sim_results$ROCresDir, 
                              sim_results$ROCresmr), function(x)x$auc)), ncol = 3, byrow = FALSE)
colnames(auc) <- c("PG", "Dir", "MR")
boxplot(auc[,c(3,2,1)], ylim = c(0.5,1), col = c("blue", "orange", "darkgreen"))

auc_long <- as.data.frame(melt(auc[,c(3,2,1)]))
colnames(auc_long) <- c("Iteration", "Method", "AUC")

gg4 <- ggplot(data = auc_long, aes(x = Method, y = AUC, fill = Method)) + geom_boxplot(outlier.colour="black", outlier.shape=16,
                                                                                       outlier.size=2, notch=FALSE)
gg4 <- gg4 + theme_classic() + scale_fill_manual(values=c("blue", "orange", "darkgreen", "darkorchid4"))
gg4 <- gg4 + theme(legend.position="none") + scale_y_continuous(limits = c(0.3,1))
gg4 <- gg4 + ggtitle("Systematic batch effects") + xlab("Method") + ylab("AUC")
gg4 <- gg4 + scale_x_discrete(labels = c("MR", "Dirichlet", "Polya-Gamma"))
gg4


plot(ROC_bayes, type = "b", col = "darkgreen", lwd = 1.5, pch = 19, asp = NA, legacy.axes = T, 
     xlab = "False Positive Rate",
     ylab = "True Positive Rate",
     main = "ROC curve demonstrating proposed methods outperform current approach")
lines(ROC_mr, type = "b", col = "blue", lwd = 1.5, pch = 19, asp = NA)
lines(ROC_dir,type = "b", col = "orange", lwd = 1.5, pch = 19, asp = NA)
legend("bottomright", c("MR", "Dirichlet", "Polya-Gamma"), col = c("blue", "orange", "darkgreen"), pch = 19)

dirnumtrue <- dirnumfalse <- pgnumtrue <- pgnumfalse <- mrnumtrue <- mrnumfalse <- vector(length = 10)


for (l in seq.int(10)){
  
  probdir <- sim_results$probdirres[[l]]
  probpg <- sim_results$probpgres[[l]]
  mrres <- sim_results$mrreslist[[l]]
  perm1_names <- sim_results$perm1_nameslist[[l]]
  
  dirnumtrue[l] <- sum(names(probdir[order(probdir)])[(1 - probdir[order(probdir)]/ncol(resdyndir$allocOut[[1]])) > 0.999] %in% perm1_names)
  dirnumfalse[l] <- sum(!names(probdir[order(probdir)])[(1 - probdir[order(probdir)]/ncol(resdyndir$allocOut[[1]])) > 0.999] %in% perm1_names)
  pgnumtrue[l] <- sum(names(probpg[order(probpg)])[(1 - probpg[order(probpg)]/ncol(resdynpg$allocOut[[1]])) > 0.999] %in% perm1_names)
  pgnumfalse[l] <- sum(!names(probpg[order(probpg)])[(1 - probpg[order(probpg)]/ncol(resdynpg$allocOut[[1]])) > 0.999] %in% perm1_names)
  mrnumtrue[l] <- sum(names(mrres$Mscore[(mrres$Mscore > 2) & mrres$Rscore > 0.9 ] ) %in% perm1_names)
  mrnumfalse[l] <- sum(!names(mrres$Mscore[(mrres$Mscore > 2) & mrres$Rscore > 0.9 ] ) %in% perm1_names)
}

numtrue <- cbind(dirnumtrue, pgnumtrue, mrnumtrue, pg1numtrue)

true_long <- as.data.frame(melt(numtrue[,c(3,1,2)]))
colnames(true_long) <- c("Iteration", "Method", "Relocalisations")

gg5 <- ggplot(data = true_long, aes(x = Method, y = Relocalisations, fill = Method)) + geom_boxplot(outlier.colour="black", outlier.shape=16,
                                                                                                    outlier.size=2, notch=FALSE)
gg5 <- gg5 + theme_classic() + scale_fill_manual(values=c("blue", "orange", "darkgreen", "darkorchid4"))
gg5 <- gg5 + theme(legend.position="none")
gg5 <- gg5 + ggtitle("Systematic batch effects") + xlab("Method") + ylab("Number of Relocalisations")
gg5 <- gg5 + scale_x_discrete(labels = c("MR", "Dirichlet", "Polya-Gamma","Polya-Gamma unshrunk"))
gg5

gg4/gg5




load("C:/Users/OllyC/Desktop/2019-DiffLoc/diffloc_sim_wild_perm.rda") # cluster specific bootstrap + permutation

auc <- matrix(unlist(lapply(c(sim_results$ROCresBayes,
                              sim_results$ROCresDir, 
                              sim_results$ROCresmr), function(x)x$auc)), ncol = 3, byrow = FALSE)
colnames(auc) <- c("PG", "Dir", "MR")
boxplot(auc[,c(3,2,1)], ylim = c(0.5,1), col = c("blue", "orange", "darkgreen"))

auc_long <- as.data.frame(melt(auc[,c(3,2,1)]))
colnames(auc_long) <- c("Iteration", "Method", "AUC")

gg6 <- ggplot(data = auc_long, aes(x = Method, y = AUC, fill = Method)) + geom_boxplot(outlier.colour="black", outlier.shape=16,
                                                                                       outlier.size=2, notch=FALSE)
gg6 <- gg6 + theme_classic() + scale_fill_manual(values=c("blue", "orange", "darkgreen", "darkorchid4"))
gg6 <- gg6 + theme(legend.position="none") + scale_y_continuous(limits = c(0.3,1))
gg6 <- gg6 + ggtitle("Fraction Swapping") + xlab("Method") + ylab("AUC")
gg6 <- gg6 + scale_x_discrete(labels = c("MR", "Dirichlet", "Polya-Gamma"))
gg6


plot(ROC_bayes, type = "b", col = "darkgreen", lwd = 1.5, pch = 19, asp = NA, legacy.axes = T, 
     xlab = "False Positive Rate",
     ylab = "True Positive Rate",
     main = "ROC curve demonstrating proposed methods outperform current approach")
lines(ROC_mr, type = "b", col = "blue", lwd = 1.5, pch = 19, asp = NA)
lines(ROC_dir,type = "b", col = "orange", lwd = 1.5, pch = 19, asp = NA)
legend("bottomright", c("MR", "Dirichlet", "Polya-Gamma"), col = c("blue", "orange", "darkgreen"), pch = 19)

dirnumtrue <- dirnumfalse <- pgnumtrue <- pgnumfalse <- mrnumtrue <- mrnumfalse <- vector(length = 10)


for (l in seq.int(10)){
  
  probdir <- sim_results$probdirres[[l]]
  probpg <- sim_results$probpgres[[l]]
  mrres <- sim_results$mrreslist[[l]]
  perm1_names <- sim_results$perm1_nameslist[[l]]
  
  dirnumtrue[l] <- sum(names(probdir[order(probdir)])[(1 - probdir[order(probdir)]/ncol(resdyndir$allocOut[[1]])) > 0.999] %in% perm1_names)
  dirnumfalse[l] <- sum(!names(probdir[order(probdir)])[(1 - probdir[order(probdir)]/ncol(resdyndir$allocOut[[1]])) > 0.999] %in% perm1_names)
  pgnumtrue[l] <- sum(names(probpg[order(probpg)])[(1 - probpg[order(probpg)]/ncol(resdynpg$allocOut[[1]])) > 0.999] %in% perm1_names)
  pgnumfalse[l] <- sum(!names(probpg[order(probpg)])[(1 - probpg[order(probpg)]/ncol(resdynpg$allocOut[[1]])) > 0.999] %in% perm1_names)
  mrnumtrue[l] <- sum(names(mrres$Mscore[(mrres$Mscore > 2) & mrres$Rscore > 0.9 ] ) %in% perm1_names)
  mrnumfalse[l] <- sum(!names(mrres$Mscore[(mrres$Mscore > 2) & mrres$Rscore > 0.9 ] ) %in% perm1_names)
}

numtrue <- cbind(dirnumtrue, pgnumtrue, mrnumtrue, pg1numtrue)

true_long <- as.data.frame(melt(numtrue[,c(3,1,2)]))
colnames(true_long) <- c("Iteration", "Method", "Relocalisations")

gg7 <- ggplot(data = true_long, aes(x = Method, y = Relocalisations, fill = Method)) + geom_boxplot(outlier.colour="black", outlier.shape=16,
                                                                                                    outlier.size=2, notch=FALSE)
gg7 <- gg7 + theme_classic() + scale_fill_manual(values=c("blue", "orange", "darkgreen", "darkorchid4"))
gg7 <- gg7 + theme(legend.position="none")
gg7 <- gg7 + ggtitle("Fraction Swapping") + xlab("Method") + ylab("Number of Relocalisations")
gg7 <- gg7 + scale_x_discrete(labels = c("MR", "Dirichlet", "Polya-Gamma","Polya-Gamma unshrunk"))
gg7

gg6/gg7

load("C:/Users/OllyC/Desktop/2019-DiffLoc/diffloc_sim_wild_batch_perm.rda") # all effects, systematic batch


auc <- matrix(unlist(lapply(c(sim_results$ROCresBayes,
                              sim_results$ROCresDir, 
                              sim_results$ROCresmr), function(x)x$auc)), ncol = 3, byrow = FALSE)
colnames(auc) <- c("PG", "Dir", "MR")
boxplot(auc[,c(3,2,1)], ylim = c(0.5,1), col = c("blue", "orange", "darkgreen"))

auc_long <- as.data.frame(melt(auc[,c(3,2,1)]))
colnames(auc_long) <- c("Iteration", "Method", "AUC")

gg8 <- ggplot(data = auc_long, aes(x = Method, y = AUC, fill = Method)) + geom_boxplot(outlier.colour="black", outlier.shape=16,
                                                                                       outlier.size=2, notch=FALSE)
gg8 <- gg8 + theme_classic() + scale_fill_manual(values=c("blue", "orange", "darkgreen", "darkorchid4"))
gg8 <- gg8 + theme(legend.position="none") + scale_y_continuous(limits = c(0.3,1))
gg8 <- gg8 + ggtitle("Batch effects and fraction swapping") + xlab("Method") + ylab("AUC")
gg8 <- gg8 + scale_x_discrete(labels = c("MR", "Dirichlet", "Polya-Gamma"))
gg8


plot(ROC_bayes, type = "b", col = "darkgreen", lwd = 1.5, pch = 19, asp = NA, legacy.axes = T, 
     xlab = "False Positive Rate",
     ylab = "True Positive Rate",
     main = "ROC curve demonstrating proposed methods outperform current approach")
lines(ROC_mr, type = "b", col = "blue", lwd = 1.5, pch = 19, asp = NA)
lines(ROC_dir,type = "b", col = "orange", lwd = 1.5, pch = 19, asp = NA)
legend("bottomright", c("MR", "Dirichlet", "Polya-Gamma"), col = c("blue", "orange", "darkgreen"), pch = 19)

dirnumtrue <- dirnumfalse <- pgnumtrue <- pgnumfalse <- mrnumtrue <- mrnumfalse <- vector(length = 10)


for (l in seq.int(10)){
  
  probdir <- sim_results$probdirres[[l]]
  probpg <- sim_results$probpgres[[l]]
  mrres <- sim_results$mrreslist[[l]]
  perm1_names <- sim_results$perm1_nameslist[[l]]
  
  dirnumtrue[l] <- sum(names(probdir[order(probdir)])[(1 - probdir[order(probdir)]/ncol(resdyndir$allocOut[[1]])) > 0.999] %in% perm1_names)
  dirnumfalse[l] <- sum(!names(probdir[order(probdir)])[(1 - probdir[order(probdir)]/ncol(resdyndir$allocOut[[1]])) > 0.999] %in% perm1_names)
  pgnumtrue[l] <- sum(names(probpg[order(probpg)])[(1 - probpg[order(probpg)]/ncol(resdynpg$allocOut[[1]])) > 0.999] %in% perm1_names)
  pgnumfalse[l] <- sum(!names(probpg[order(probpg)])[(1 - probpg[order(probpg)]/ncol(resdynpg$allocOut[[1]])) > 0.999] %in% perm1_names)
  mrnumtrue[l] <- sum(names(mrres$Mscore[(mrres$Mscore > 2) & mrres$Rscore > 0.9 ] ) %in% perm1_names)
  mrnumfalse[l] <- sum(!names(mrres$Mscore[(mrres$Mscore > 2) & mrres$Rscore > 0.9 ] ) %in% perm1_names)
}

numtrue <- cbind(dirnumtrue, pgnumtrue, mrnumtrue, pg1numtrue)

true_long <- as.data.frame(melt(numtrue[,c(3,1,2)]))
colnames(true_long) <- c("Iteration", "Method", "Relocalisations")

gg9 <- ggplot(data = true_long, aes(x = Method, y = Relocalisations, fill = Method)) + geom_boxplot(outlier.colour="black", outlier.shape=16,
                                                                                                   outlier.size=2, notch=FALSE)
gg9 <- gg9 + theme_classic() + scale_fill_manual(values=c("blue", "orange", "darkgreen", "darkorchid4"))
gg9 <- gg9 + theme(legend.position="none")
gg9 <- gg9 + ggtitle("Batch effects and fraction swapping") + xlab("Method") + ylab("Number of Relocalisations")
gg9 <- gg9 + scale_x_discrete(labels = c("MR", "Dirichlet", "Polya-Gamma","Polya-Gamma unshrunk"))
gg9

gg8/gg9



# cluster specific bootstrap
# random batch effect
# systematic batch effect
# fraction permutations
# all effects systematic batch effects
(gg + gg1 + gg4 + gg6 +  gg8)/(gg2 + gg3 + gg5 + gg7 + gg9) 
