# Plots for CEEA presentation, June 2015

# Run the following code once to make sure you have the right R packages
#install.packages(c('gdata', 'data.table', 'car', 'boot'))

# Please note, results might differ slightly on your computer, however, they will
# not be much different to those presented in the paper.
# That's because of the random number generator used for bootstrapping.

# Just plot the raw data
bitmap('raw-data-tests.png', type="png256", width=14, height=14*3/4, res=300, pointsize=14)
par(mar=c(4.2, 4.2, 1.5, 0.5))  # (bottom, left, top, right); defaults are par(mar=c(5, 4, 4, 2) + 0.1)
par(cex.lab=1.5, cex.main=1.2, cex.sub=1.8, cex.axis=1.8 )
layout(matrix(c(1,2,3,4), 2, 2))

for (i in c(1,2,3,4)) {
  require(gdata)
  library(gdata)
  
  # Download the data from http://yint.org/unlimited-time-tests
  # and save as an Excel file in the local directory. 
  # Use an R library to read in this data.
  data <- read.xls('Unlimited-time--Raw-data.xlsx', sheet=i)
  colnames(data) <- c("Time", "Grade", "Details")
  
  # Re-order the rows by "Time" order
  library(data.table)
  data <- data.table(data)
  data[order(data$Time)]
  
  plot(data$Time, data$Grade, 
       xlab="Time duration [minutes]",        
       ylim=c(31, 101),
       ylab="Student's percentage",
       main=""
  )
  
  minx = min(data$Time)
  #text(minx+15, 50, paste("r =", round(cor(data$Time, data$Grade),2)), cex=2)
  text(minx+15, 95, paste("Test", i), cex=3)  
}
dev.off()

# Add the smoother and regression model
bitmap('raw-data-tests-with-model.png', type="png256", width=14, height=14*3/4, res=300, pointsize=14)
par(mar=c(4.2, 4.2, 1.5, 0.5))  # (bottom, left, top, right); defaults are par(mar=c(5, 4, 4, 2) + 0.1)
par(cex.lab=1.5, cex.main=1.2, cex.sub=1.8, cex.axis=1.8 )
layout(matrix(c(1,2,3,4), 2, 2))

for (i in c(1,2,3,4)) {
  require(gdata)
  library(gdata)
  # Download the data from http://yint.org/unlimited-time-tests
  # and save as an Excel file. Use an R library to read in this data.
  data <- read.xls('Unlimited-time--Raw-data.xlsx', sheet=i)
  colnames(data) <- c("Time", "Grade", "Details")
  
  # Re-order the rows by "Time" order
  library(data.table)
  data <- data.table(data)
  data[order(data$Time)]
  
  plot(data$Time, data$Grade, 
       xlab="Time duration [minutes]",        
       ylim=c(31, 101),
       ylab="Student's percentage",
       main=""
  )
  
  minx = min(data$Time)
  text(minx+15, 50, paste("r =", round(cor(data$Time, data$Grade),2)), cex=2)
  text(minx+15, 95, paste("Test", i), cex=3)
  
  model <- lm(data$Grade ~ data$Time)  
  abline(model, col="blue", lwd=2)
  
  lines(lowess(data$Grade ~ data$Time), lwd=2, col="red", lty=2)
}
dev.off()

bitmap('histograms.png', type="png256", width=14, height=14*3/4, res=300, pointsize=14)
par(mar=c(4.2, 4.2, 1.5, 0.5))  # (bottom, left, top, right); defaults are par(mar=c(5, 4, 4, 2) + 0.1)
par(cex.lab=1.5, cex.main=1.2, cex.sub=1.8, cex.axis=1.8 )
layout(matrix(c(1,2,3,4), 2, 2))
for (i in 1:4) {
  require(gdata)
  library(gdata)
  # Download the data from http://yint.org/unlimited-time-tests
  # and save as an Excel file. Use an R library to read in this data.
  data <- read.xls('Unlimited-time--Raw-data.xlsx', sheet=i)
  colnames(data) <- c("Time", "Grade", "Details")
  print(i)
  # Re-order the rows by "Time" order
  library(data.table)
  data <- data.table(data)
  data[order(data$Time)]
  
  # Code based on that from: http://www.statmethods.net/advstats/bootstrapping.html
  library(boot)
  coefficient <- function(formula, data, indices) {
    d <- data[indices,]
    fit <- lm(formula, data=d)
    return(as.numeric(fit$coefficients[2]))
  } 
  
  SE <- function(formula, data, indices) {
    d <- data[indices,]
    fit <- lm(formula, data=d)
    return(summary(fit)$sigma)
  } 
  
  # Bootstrap with 10000 replications.
  # Change the "statistic" from "coefficient" if you'd like to generate 
  # different bootstrap statistics (e.g. for standard error)
  results <- boot(data=data, statistic=coefficient, 
                  R=10000, formula=Grade~Time)
  
  hist(results$t, xlab=expression('Slope coefficient, b'[T]), 
       main=paste("Bootstrap slope coefficient for test", i)
  )
   
  abline(v=boot.ci(results, type="bca")$bca[4], lwd=3, col="blue", lty=3)
  abline(v=boot.ci(results, type="bca")$bca[5], lwd=3, col="blue", lty=3)
  abline(v=0, col="red", lwd=4)
  abline(v=mean(c(boot.ci(results, type="bca")$bca[4], boot.ci(results, type="bca")$bca[5])), col="blue", lwd=6)
}
dev.off()

bitmap('confidence-intervals.png', type="png256", width=14, height=14*3/4, res=300, pointsize=14)
par(mar=c(4.2, 4.2, 1.5, 0.5))  # (bottom, left, top, right); defaults are par(mar=c(5, 4, 4, 2) + 0.1)
par(cex.lab=1.5, cex.main=1.2, cex.sub=1.8, cex.axis=1.8 )
plot(c(0,0), c(0.5,8.5), type="l", col="red", lwd=4, ylab="Test number", 
     xlab=expression('Confidence interval for slope coefficient, b'[T]),
     xlim=c(-0.5, 0.5))
for (i in 1:8) {
  require(gdata)
  library(gdata)
  # Download the data from http://yint.org/unlimited-time-tests
  # and save as an Excel file. Use an R library to read in this data.
  data <- read.xls('Unlimited-time--Raw-data.xlsx', sheet=i)
  colnames(data) <- c("Time", "Grade", "Details")
  print(i)
  # Re-order the rows by "Time" order
  library(data.table)
  data <- data.table(data)
  data[order(data$Time)]
  
  # Code based on that from: http://www.statmethods.net/advstats/bootstrapping.html
  library(boot)
  coefficient <- function(formula, data, indices) {
    d <- data[indices,]
    fit <- lm(formula, data=d)
    return(as.numeric(fit$coefficients[2]))
  } 
  
  SE <- function(formula, data, indices) {
    d <- data[indices,]
    fit <- lm(formula, data=d)
    return(summary(fit)$sigma)
  } 
  
  # Bootstrap with 10000 replications.
  # Change the "statistic" from "coefficient" if you'd like to generate 
  # different bootstrap statistics (e.g. for standard error)
  results <- boot(data=data, statistic=coefficient, 
                  R=10000, formula=Grade~Time)
  lines(c(boot.ci(results, type="bca")$bca[4], boot.ci(results, type="bca")$bca[5]), 
       c(i,i), lwd=5, type="b")
  text(boot.ci(results, type="bca")$bca[5]+0.1, i, paste("Test", i), cex=1.5)
  
}
dev.off()