#COMMANDS USED TO OBTAIN GRAPHS AND TABLES

setwd("UNIT:/FOLDER") #set your own folder
source("functions.R")
# Uncomment and execute sentences reading data of interest
data1 <- read.table(file="prodCICLOPA1.txt", header=T)
#data2 <- read.table(file="prodCICLOPA2.txt", header=T)
#data3 <- read.table(file="prodCICLOPB1.txt", header=T)
#data4 <- read.table(file="prodCICLOPB2.txt", header=T)

# View data:
fix(data1)
produc.vec <- data1[,1]
#produc.vec <- data2[,1]
#produc.vec <- data3[,1]
#produc.vec <- data4[,1]

# COMMAND: generation of 10000 production bootstrap replicates (CPI values are randomly
# generated. according to a normal distribution with mean (120+180)/2 = 150
# and standard deviation (180 - 150) / qnorm(0.975), the default option)
boot.prods <- bootReplicates.prod(produc.vec)
boot.prods[,1:10]
rownames(boot.prods)

# Percentile confidence interval from these resamples
bootPercCI.prod(boot.replicates = boot.prods)
# BCa confidence interval (original data vector is required)
bootBCaCI.prod(produc.vec, boot.replicates = boot.prods)                                


# COMMAND: Comparing two or more datasets:
# By means of percentile interval and histogram:
compareProds(data1[,1], data2[,1], name1 = "A1", name2 = "A2")
# By means of percentile interval and box-plot:
compareProds(data1[,1], data2[,1], name1 = "A1", name2 = "A2", graphic = "boxplot")
# By means of BCa interval and histogram:
compareProds(data1[,1], data2[,1], name1 = "A1", name2 = "A2", confInt = "BCa")
# By means of BCa interval and box-plot:
compareProds(data1[,1], data2[,1], data3[,1], data4[,1], name1 = "C-Bef", name2 = "C-Aft", name3 = "I-Bef", name4 = "I-Aft", confInt = "BCa", graphic = "boxplot") #*********************************THIS GENERATES FIG. 3 IN THE MANUSCRIPT AND VALUES OF TABLE 4******************************************************************************************************************************************************************************

# Similar comparisons but from a previously generated set of resamples:       
boot.prods1 <- bootReplicates.prod(data1[,1])
boot.prods2 <- bootReplicates.prod(data2[,1])
boot.prods3 <- bootReplicates.prod(data3[,1])
boot.prods4 <- bootReplicates.prod(data4[,1])
compareProds(boot1 = boot.prods1, boot2 = boot.prods2, name1 = "A1", name2 = "A2")
compareProds(boot1 = boot.prods1, boot2 = boot.prods2, name1 = "A1", name2 = "A2", graphic = "boxplot")
compareProds(exemple1[,1], exemple2[,1], boot.prods1, boot.prods2, "A1", "A2", confInt = "BCa")
compareProds(data1[,1], data2[,1],data3[,1], data4[,1],
  boot.prods1, boot.prods2, boot.prods3, boot.prods4,
  "C-Bef", "C-Aft", "I-Bef", "I-Aft", 
  confInt = "BCa", graphic = "boxplot")   #FIG. 3 IN THE MANUSCRIPT ALSO.
  
compareProds(data1[,1], data2[,1],data3[,1], data4[,1],
  boot.prods1, boot.prods2, boot.prods3, boot.prods4,
  "C-Bef", "C-Aft", "I-Bef", "I-Aft", 
  confInt = "BCa")      # Histogram