R/GiR_PP_Plots.R

In matthewjdenny/CCAS: Communication Content and Network Structure

#' @title Make Nice GiR PP-Plots
#' @description Makes Nice PP Plots with Getting it Right Results.
#'
#' @param forward_samples A dataframe with statistic generated by
#' test_internal_functions() with Getting_It_Right = TRUE.
#' @param backward_samples A dataframe with statistic generated by
#' test_internal_functions() with Getting_It_Right = TRUE.
#' @param generate_PP_plot Defaults to TRUE. If FALSE, then QQ plots are
#' generated.
#' @return A bunch of plots.
#' @export
GiR_PP_Plots <- function(forward_samples,
                         backward_samples,
                         generate_PP_plot = TRUE) {
    UMASS_BLUE <- rgb(51,51,153,155,maxColorValue = 255)
    UMASS_RED <- rgb(153,0,51,155,maxColorValue = 255)
    UMASS_GREEN <- rgb(0,102,102,255,maxColorValue = 255)
    UMASS_YELLOW <- rgb(255,255,102,255,maxColorValue = 255)
    UMASS_ORANGE <- rgb(255,204,51,255,maxColorValue = 255)

    nms <- colnames(forward_samples)

    for (i in 1:ncol(forward_samples)) {
        # don't plot MH accept rates
        if (colnames(forward_samples)[i] != "Mean_MH_Accept_Rate") {

            # start by generating reduced samples for t-tests
            if (nrow(forward_samples) > 20000) {
                thin <- seq(from = floor(nrow(forward_samples)/10), to = nrow(forward_samples),length.out = 10000)
                forward_test <- forward_samples[thin,i]
            } else {
                forward_test <- forward_samples[,i]
            }

            # start by generating reduced samples for t-tests
            if (nrow(backward_samples) > 20000) {
                thin <- seq(from = floor(nrow(backward_samples)/10), to = nrow(backward_samples),length.out = 10000)
                backward_test <- backward_samples[thin,i]
            } else {
                backward_test <- backward_samples[,i]
            }


            cat("Generating Plot",i,"of",ncol(forward_samples),"\n")

            if (generate_PP_plot) {
                forward <- forward_samples[,i]
                backward <- backward_samples[,i]
                # find all unique values across the two vectors
                uniqueValues <- sort(unique(c(forward,backward)))

                if (length(uniqueValues) > 1) {
                    # thinning
                    thin = 500
                    if (length(uniqueValues) > thin) {
                        uniqueValues <- uniqueValues[round(seq(1,length(uniqueValues),length.out = thin))]
                    }


                    # create vectors in which to store empirical quantiles
                    qx1 <- numeric(length(uniqueValues))
                    qx2 <- numeric(length(uniqueValues))

                    # loop (could be paralellized) to calculate empirical quantile position
                    # of each unique value in each vector
                    for(k in 1:length(uniqueValues)){
                        qx1[k] <- mean(forward <= uniqueValues[k])
                        qx2[k] <- mean(backward <= uniqueValues[k])
                    }

                    # calcualte autocorrelation in forward and backward chains
                    if (length(unique(forward)) > 1) {
                        ar1 <- stats::cor(forward[2:length(forward)],forward[1:(length(forward)-1)])
                    } else {
                        ar1 <- 1
                    }

                    if (length(unique(backward)) > 1) {
                        ar2 <- stats::cor(backward[2:length(backward)],backward[1:(length(backward)-1)])
                    } else {
                        ar2 <- 1
                    }

                    plot(x = qx1, y = qx2,
                         ylim = c(0,1),
                         xlim = c(0,1),
                         ylab = "Backward",
                         xlab = "Forward",
                         col = UMASS_BLUE,
                         pch = 4,
                         main = nms[i],
                         cex.lab=2,
                         cex.axis=1.4,
                         cex.main=2,cex = .5)
                    lines(x = c(0,1), y= c(0,1),
                          col = UMASS_RED,lwd = 3)
                    text(paste( "Backward Mean:", round(mean(backward_samples[,i]),4),
                                "\nForward Mean:", round(mean(forward_samples[,i]),4),
                                "\nt-test p-value:",
                                round(t.test(backward_test,
                                             forward_test)$p.value,4),
                                "\nMann-Whitney p-value:",
                                round(wilcox.test(backward_test,
                                                  forward_test)$p.value,4)),
                         x = 0.7,
                         y = 0.2,
                         cex = 1.5)
                    text(paste( "Backward Autocorr:", round(ar2,4),
                                "\nForward Autocorr:", round(ar1,4)),
                         x = 0.3,
                         y = 0.9,
                         cex = 1.5)
                }
            } else {
                all <- c(backward_samples[,i],forward_samples[,i])
                ylims <- c(min(all)-0.1*(max(abs(all))),
                           max(all)+0.1*(max(abs(all))))
                xlims <- ylims

                # set a higher number of quantiles for LSM parameters
                quantiles <- 50
                if (grepl("LSM",colnames(forward_samples)[i])) {
                    quantiles <- 1000
                }
                qqplot(x = quantile(forward_samples[,i],seq(0,1,length=quantiles)),
                       y = quantile(backward_samples[,i],seq(0,1,length=quantiles)),
                       ylim = ylims,
                       xlim = xlims,
                       ylab = "Backward",
                       xlab = "Forward",
                       col = UMASS_BLUE,
                       pch = 19,
                       main = nms[i],
                       cex.lab=2,
                       cex.axis=1.4,
                       cex.main=2)
                lines(x = xlims, y= ylims,
                      col = UMASS_RED,lwd = 3)
                text(paste( "Backward Mean:", round(mean(backward_samples[,i]),4),
                            "\nForward Mean:", round(mean(forward_samples[,i]),4),
                            "\nt-test p-value:",
                            round(t.test(backward_test,
                                         forward_test)$p.value,4),
                            "\nMann-Whitney p-value:",
                            round(wilcox.test(backward_test,
                                              forward_test)$p.value,4)),
                     x = xlims[2] - 0.3*abs(xlims[2] - xlims[1]),
                     y = ylims[1] + 0.2*abs(ylims[2] - ylims[1]),
                     cex = 1.5)
            }

        }
    }
}