FR_before_adj_comparison.R

FR_before_adj_comparison.R

# produce plots and CSV files with looping rather than pasting repeated commands

library(moments)
setwd("~/Documents/TBP_Plan_1/FR_before_adj/")

# Load data files
load("Study1a_M_F0_100_vf1_FR_before_adj.RData")
load("Study1b_M_F0_100_vf1_FR_before_adj.RData")
load("Study2b_M_F0_100_vf1_FR_before_adj_1000scen.RData")
load("Study3b_M_F0_100_vf1_FR_before_adj_1000scen.RData")
load("Study4b_M_F0_100_vf1_FR_before_adj_1000scen.RData")
load("Study5b_M_F0_100_vf1_FR_before_adj_1000scen.RData")
load("Study6b_M_F0_100_vf1_FR_before_adj_1000scen.RData")
load("Study7b_M_F0_100_vf1_trigger_100_140_FR_before_adj.RData")
load("Study8b_M_F0_100_vf1_trigger_100_140_FR_before_adj.RData")
load("Study9b_M_F0_100_vf1_trigger_100_140_FR_before_adj.RData")

#' @title obj The name of the objects to include in the plot
#' @description Create a density plot comparing the funding ratio
#'    under different policy types
#' @param obj A vector of obj names to be included in the plot
#' @param t The time to produce the plot at 
#' @param plotTitle optional overwrite for plot title
#' @param ylim y limits for plots (optional)
#' @param legendNames Titles for plots
densityPlot <- function(obj="FR.before.adj.1a", t=25,
                        plotTitle=NA, ylim=NA,
                        legendNames=substr(obj,nchar(obj)-2+1,nchar(obj))) {
    
    if(is.na(plotTitle)) {
        plotTitle <- paste0("Kernel density plot of funded ratio at t=", t)
    }
    
    numObjs <- length(obj)
    tmp <- get(obj[1]) 
    if(!is.na(ylim[1])) {
        plot(density(tmp[t,]), main=plotTitle, xlab="Funded Ratio", 
             xlim=c(50,180), ylim=ylim)
    } else {
        plot(density(tmp[t,]), main=plotTitle, xlab="Funded Ratio", 
             xlim=c(50,180))
    }
    
    legColors <- c(1)
    legLTY <- c(1)
    
    if(numObjs > 1) { # lines
        for(i in 2:numObjs) {
            tmp <- get(obj[i])
            lines(density(tmp[t,]), lty=(i+1)%%2+1, col=floor(i/2))
            
            legColors[i] <- i%%2 + 1
            legLTY[i] <- ceiling(i/2)
        }
    }
    
    # legend
    legend("topright", legendNames, col=legColors, lty=legLTY)
}

csvSummary <- function(obj="FR.before.adj.1a", t=25, fname=paste0("FR_table3_t=",t,".csv")) {
    numObjects <- length(obj)
    statTable <- matrix(0, 16, numObjects)
    cnames <- substr(obj,nchar(obj)-2+1,nchar(obj))
    
    for(i in 1:numObjects) {
        tmpstr <- gsub(cnames[i], paste0("q.",cnames[i]), obj[i])
        tmp <- get(tmpstr) # quantiles
        statTable[1:5,i] <- tmp[t,1:5] # percentiles
        
        tmp <- get(obj[i]) # normal
        statTable[6,i] <- mean(tmp[t,]) # mean
        statTable[7,i] <- sd(tmp[t,]) # variance
        statTable[8,i] <- skewness(tmp[t,]) # skewness
        statTable[9,i] <- kurtosis(tmp[t,]) # kurtosis
        tmp.ecdf <- ecdf(tmp[t,]) # ecdf function can return percentile
        statTable[10,i] <- tmp.ecdf(80) # < 80
        statTable[11,i] <- tmp.ecdf(90) - tmp.ecdf(80) # (80,90)
        statTable[12,i] <- tmp.ecdf(100) - tmp.ecdf(90) # (90,100)
        statTable[13,i] <- tmp.ecdf(110) - tmp.ecdf(100) # (100,110)
        statTable[14,i] <- tmp.ecdf(120) - tmp.ecdf(110) # (110, 120)
        statTable[15,i] <- tmp.ecdf(140) - tmp.ecdf(120) # (120, 140)
        statTable[16,i] <- 1 - tmp.ecdf(140) # > 140
    }
    statTable <- data.frame(statTable)
    colnames(statTable) <- cnames
    rownames(statTable) <- c("5th percentile","25th percentile","50th percentile",
                             "75th percentile","95th percentile","mean","SD",
                             "skewness","kurtosis","Pr(FR<=80)","Pr(80<FR<=90)",
                             "Pr(90<FR<=100)","Pr(100<FR<=110)","Pr(110<FR<=120)",
                             "Pr(120<FR<=140)","Pr(FR>140)")
    write.csv(x=statTable,file=fname)
}

# plots
for(time in c(10,25,50,100)) {
    densityPlot(obj=c("FR.before.adj.1a", "FR.before.adj.1b"), t=time, 
                legendNames=c("study 1a", "study 1b"))
    
    densityPlot(obj=c("FR.before.adj.1b", "FR.before.adj.2b",
                      "FR.before.adj.3b", "FR.before.adj.4b", "FR.before.adj.5b"), t=time,
                legendNames=c("1b (100)", "2b (90-110)", "3b (80-120)", "4b (100-120)",
                              "5b (100-140)"))
}

densityPlot(obj=c("FR.before.adj.2b", "FR.before.adj.6b"), t=25,
            legendNames=c("2b (90-110)", "6b (90-110)"))
densityPlot(obj=c("FR.before.adj.5b", "FR.before.adj.7b",
                  "FR.before.adj.8b", "FR.before.adj.9b"), t=25, 
            legendNames=c("5b (100-140)", "7b (100-140)", "8b (100-140)", "9b (100-140)"),
            ylim=c(0,0.0325))

# CSV files
csvSummary(obj=c("FR.before.adj.1b", "FR.before.adj.2b",
    "FR.before.adj.3b", "FR.before.adj.4b", "FR.before.adj.5b"), t=5)
csvSummary(obj=c("FR.before.adj.1b", "FR.before.adj.2b",
    "FR.before.adj.3b", "FR.before.adj.4b", "FR.before.adj.5b",
    "FR.before.adj.7b", "FR.before.adj.8b", "FR.before.adj.9b"), t=25)