# add_perm_p() is a function to replace ggpairs' (part of ggally) built-in 
# parametric significance-determination method with a permutations-based one.
#
# For example code used in our paper, see commented code below the function

add_perm_p <- function(data, mapping, n_perm = 1000, multiple_n = 1, multiple_method = "bonferroni", alpha = 0.05, textsize = 5, color = "black", ...) {
    # Credit to https://stackoverflow.com/a/61687822
    x = eval_data_col(data, mapping$x)
    y = eval_data_col(data, mapping$y)
    
    p_max <- 1/n_perm
    p_max_corr <- p.adjust(p_max, method = multiple_method, n=multiple_n)
    r_observed = cor(x, y)
    r_text <- paste("r =", round(r_observed, digits = 3))
    r_perm <- sapply(1:n_perm, FUN = function(i) cor(x=x, y=sample(y)))
    p_perm <- sum(abs(r_perm) >= abs(r_observed))/(n_perm + 1)
    p_corr <- p.adjust(p_perm, method = multiple_method, n=multiple_n)
    
    if(p_corr < p_max_corr) {
        p_text <- paste("p <", round(p_max_corr, digits = 4))
    } else if(p_corr > alpha) {
        p_text <- paste("p >", alpha)
        color = 'darkgray'
    } else {
        p_text <- paste("p =", round(p_corr, digits = 4))
    }

    text = paste(r_text, "\n", p_text)
    
    textplot = ggplot(data = data, mapping = mapping) + annotate(x = 0.5, y = 0.8, label = text, geom="Text", size=textsize, color = color) + theme_void()
    
    return(textplot)
}

## Load GGally library
# library(GGally)
#
## Load merge CSV of HCP's restricted and unrestricted data and subject IDs
# merged <- read.csv("./merged.csv")
# ids <- scan("./subject.ids")
#
## Extract relevant columns
# corrdata <- merged[which(merged$Subject %in% ids), 
#                    c('BMI', 'Flanker_AgeAdj', 'CardSort_AgeAdj', 'DDisc_AUC_200', 'DDisc_AUC_40K', 'PMAT24_A_CR', 'WM_Task_2bk_Acc')]
# corrdata <- corrdata[complete.cases(corrdata), ]
# colnames(corrdata) <- c('BMI', 'Flanker', 'CardSort', 'DDisc_200', 'DDisc_40K', 'PMAT24', 'WM_2bk') # Shorter column names for a neater plot
#
## Determine the number of comparisons to correct for
## We can either concentrate on one measure of interest:
##multiple_n = ncol(corrdata) - 1
##n_perm = 1000 * multiple_n # We want to estimate a corrected p of <= 0.001
## Or be more conservative and correct for all comparisons (which is what we did in our paper):
#n = ncol(corrdata)
#multiple_n = (n * (n-1)) /2 # Calculate the number of cells in the n*n matrix without the diagonal
#n_perm = 1000 * multiple_n # We want to estimate a corrected p of <= 0.001
#
## Produce plot
#splom <- ggpairs(corrdata,
#                 lower = list(continuous = wrap("smooth", method = "loess", alpha=0.2, size=1.25, color="dodgerblue3")),
#                 upper = list(continuous = wrap(add_perm_p, n_perm=n_perm, multiple_n=multiple_n)),
#                 diag = list(continuous = wrap("densityDiag", fill="chartreuse4", alpha=0.3))
#) + theme(axis.text = element_text(size=12), strip.text.x = element_text(size = 14), strip.text.y = element_text(size = 12))
#splom
#
## Save plot
#ggsave("/tmp/corr_measures_splom.png", plot = splom, width=35, height=23, units = "cm", dpi=600)