---
title: "Axonal Length"
author: "Sergio D?ez Hermano"
date: '`r format(Sys.Date(),"%e de %B, %Y")`'
output:
  html_document:
      highlight: tango
      toc: yes
      toc_depth: 4
      toc_float:
        collapsed: no
---

```{r setup, include=FALSE}
require(knitr)
# include this code chunk as-is to set options
opts_chunk$set(comment = NA, prompt = FALSE, fig.height=5, fig.width=5, dpi=300, fig.align = "center", echo = TRUE, message = FALSE, warning = FALSE, cache=FALSE)
Sys.setlocale("LC_TIME", "C")
```


##Load libraries and functions

```{r}
library(R.matlab)
library(lattice)
library(dplyr)
library(ggplot2)
library(ggridges)
library(hrbrthemes)
library(viridis)
library(ggpubr)
library(mltools)
library(data.table)
library(caret)
library(interactions)
# library(performance)
library(MASS)
library(geepack)
library(pstools)
library(MXM)
library(rmcorr)
library(multcomp)
library(Hmisc)


options(digits = 10)

# https://rstudio-pubs-static.s3.amazonaws.com/297778_5fce298898d64c81a4127cf811a9d486.html


################
# Bootstrap PI #
################

the.replication <- function(reg, s, x_Np1 = 0){
  # Make bootstrap residuals
  ep.star <- sample(s, size=length(reg$residuals),replace=TRUE)

  # Make bootstrap Y
  y.star <- fitted(reg) + ep.star

  # Do bootstrap regression
  lp <- model.frame(reg)[,2]
  nt <- model.frame(reg)[,3]
  bs.reg <- lm(y.star ~ lp:nt - 1)

  # Create bootstrapped adjusted residuals
  bs.lev <- influence(bs.reg)$hat
  bs.s   <- residuals(bs.reg)/sqrt(1-bs.lev)
  bs.s   <- bs.s - mean(bs.s)

  # Calculate draw on prediction error
  # xb.xb <- coef(my.reg)["(Intercept)"] - coef(bs.reg)["(Intercept)"] 
  xb.xb <- (coef(my.reg)[1] - coef(bs.reg)[1])*x_Np1
  return(unname(xb.xb + sample(bs.s, size=1)))
  
}


############################
# Bootstrap PI generalized #
############################

the.replication.gen <- function(reg, s, axType, x_Np1 = 0){
  # Make bootstrap residuals
  ep.star <- sample(s, size=length(reg$residuals), replace=TRUE)

  # Make bootstrap Y
  y.star <- fitted(reg) + ep.star

  # Do bootstrap regression
  lp <- model.frame(reg)[,2]
  nt <- model.frame(reg)[,3]
  bs.reg <- lm(y.star ~ lp:nt - 1)

  # Create bootstrapped adjusted residuals
  bs.lev <- influence(bs.reg)$hat
  bs.s   <- residuals(bs.reg)/sqrt(1-bs.lev)
  bs.s   <- bs.s - mean(bs.s)

  # Calculate draw on prediction error
  # xb.xb <- coef(my.reg)["(Intercept)"] - coef(bs.reg)["(Intercept)"] 
  xb.xb <- (coef(my.reg)[axType] - coef(bs.reg)[axType])*x_Np1
  return(unname(xb.xb + sample(bs.s, size=1)))
  
}


##############################
# Bootstrap PI per axon type #
##############################

the.replication.type <- function(reg, s, x_Np1 = 0){
  # Make bootstrap residuals
  ep.star <- sample(s, size=length(reg$residuals),replace=TRUE)

  # Make bootstrap Y
  y.star <- fitted(reg) + ep.star

  # Do bootstrap regression
  lp <- model.frame(reg)[,2]
  bs.reg <- lm(y.star ~ lp - 1)

  # Create bootstrapped adjusted residuals
  bs.lev <- influence(bs.reg)$hat
  bs.s   <- residuals(bs.reg)/sqrt(1-bs.lev)
  bs.s   <- bs.s - mean(bs.s)

  # Calculate draw on prediction error
  # xb.xb <- coef(my.reg)["(Intercept)"] - coef(bs.reg)["(Intercept)"] 
  xb.xb <- (coef(my.reg) - coef(bs.reg))*x_Np1
  return(unname(xb.xb + sample(bs.s, size=1)))
  
}

##############################
# Stratified random sampling #
##############################

stratified <- function(df, group, size, select = NULL, 
                       replace = FALSE, bothSets = FALSE) {
  if (is.null(select)) {
    df <- df
  } else {
    if (is.null(names(select))) stop("'select' must be a named list")
    if (!all(names(select) %in% names(df)))
      stop("Please verify your 'select' argument")
    temp <- sapply(names(select),
                   function(x) df[[x]] %in% select[[x]])
    df <- df[rowSums(temp) == length(select), ]
  }
  df.interaction <- interaction(df[group], drop = TRUE)
  df.table <- table(df.interaction)
  df.split <- split(df, df.interaction)
  if (length(size) > 1) {
    if (length(size) != length(df.split))
      stop("Number of groups is ", length(df.split),
           " but number of sizes supplied is ", length(size))
    if (is.null(names(size))) {
      n <- setNames(size, names(df.split))
      message(sQuote("size"), " vector entered as:\n\nsize = structure(c(",
              paste(n, collapse = ", "), "),\n.Names = c(",
              paste(shQuote(names(n)), collapse = ", "), ")) \n\n")
    } else {
      ifelse(all(names(size) %in% names(df.split)),
             n <- size[names(df.split)],
             stop("Named vector supplied with names ",
                  paste(names(size), collapse = ", "),
                  "\n but the names for the group levels are ",
                  paste(names(df.split), collapse = ", ")))
    }
  } else if (size < 1) {
    n <- round(df.table * size, digits = 0)
  } else if (size >= 1) {
    if (all(df.table >= size) || isTRUE(replace)) {
      n <- setNames(rep(size, length.out = length(df.split)),
                    names(df.split))
    } else {
      message(
        "Some groups\n---",
        paste(names(df.table[df.table < size]), collapse = ", "),
        "---\ncontain fewer observations",
        " than desired number of samples.\n",
        "All observations have been returned from those groups.")
      n <- c(sapply(df.table[df.table >= size], function(x) x = size),
             df.table[df.table < size])
    }
  }
  temp <- lapply(
    names(df.split),
    function(x) df.split[[x]][sample(df.table[x],
                                     n[x], replace = replace), ])
  set1 <- do.call("rbind", temp)
  
  if (isTRUE(bothSets)) {
    set2 <- df[!rownames(df) %in% rownames(set1), ]
    list(SET1 = set1, SET2 = set2)
  } else {
    set1
  }
}

```


##Load data

```{r}
# Get file paths
axonNames <- list.files(paste("ProyeccionData/", sep=""), full.names=T)

# Load matlab file
axonFiles <- lapply(axonNames, function(x) readMat(x))
names(axonFiles) <- c("Type1", "Type2", "Type3", "Type4")

# VAR NAMES
# REAL_LENGTH, PROYECTION_LENGTH_XY, PROYECTION_LENGTH_XZ, PROYECTION_LENGTH_YZ

# Extract data
realLength   <- lapply(axonFiles, function(x) x$REAL.LENGTH)
planeXY      <- lapply(axonFiles, function(x) x$PROYECTION.LENGTH.XY)
planeXZ      <- lapply(axonFiles, function(x) x$PROYECTION.LENGTH.XZ)
planeYZ      <- lapply(axonFiles, function(x) x$PROYECTION.LENGTH.YZ)

# Get number of neurons per type
numberTypes  <- unlist(lapply(realLength, function(x) length(x)))

# Store in data frames by plane
xyData <- data.frame(id = 1:length(unlist(realLength)), 
                      LR = unlist(realLength), LP = unlist(planeXY), Plane = rep("XY", sum(numberTypes)),
                      neurType = rep(c("Type1", "Type2", "Type3", "Type4"), numberTypes))

xzData <- data.frame(id = 1:length(unlist(realLength)), 
                      LR = unlist(realLength), LP = unlist(planeXZ), Plane = rep("XZ", sum(numberTypes)),
                      neurType = rep(c("Type1", "Type2", "Type3", "Type4"), numberTypes))

yzData <- data.frame(id = 1:length(unlist(realLength)), 
                      LR = unlist(realLength), LP = unlist(planeYZ), Plane = rep("YZ", sum(numberTypes)),
                      neurType = rep(c("Type1", "Type2", "Type3", "Type4"), numberTypes))

# Merge into a single data frame and sort by id
allData <- do.call("rbind", list(xyData, xzData, yzData))
allData <- allData[order(allData$id),]
# Add different number for every neuron-plane combination
allData$NeurPlane <- c(rep(c(1,2,3), numberTypes[1]), rep(c(4,5,6), numberTypes[2]), 
                       rep(c(7,8,9), numberTypes[3]), rep(c(10,11,12), numberTypes[4]))
allData$interact <-interaction(allData$neurType, allData$Plane, lex.order=T)
allData$idR <- 1:dim(allData)[1]

# Create data frame w/o Type4
dataNo4 <- allData[allData$neurType != "Type4", ]; dataNo4$neurType <- factor(dataNo4$neurType)

# Data in short format
dataShort <- data.frame(id = 1:length(unlist(realLength)), 
                      LR = unlist(realLength), XY = unlist(planeXY), XZ = unlist(planeXZ), YZ = unlist(planeYZ),
                      neurType = rep(c("Type1", "Type2", "Type3", "Type4"), numberTypes))
```


##CV-Bootstrap Pred Interval

RUN ONLY ONCE, AFTERWARDS LOAD THE FILE IN MEAN ERRORS CHUNK

```{r}
# Ensure reproducibility
set.seed(123456)

axList <- list()

for (h in unique(allData$neurType)) {
  
  # Subset axon type
  dataCV <- allData[allData$neurType == h, ]
  
  # List to store CV repetitions
  cvList <- list ()
  cvRand <- list()
  
  # Stratified random 5-fold sets
  setList <- list()
  sampSize <- 0.2*length(unique(dataCV$id))
  
  groupDat1 <- dataCV %>% group_by(Plane)
  set1<- sample_n(groupDat1, sampSize)
  
  groupDat2 <- groupDat1[-set1$idR, ]
  set2 <- sample_n(groupDat2, sampSize)
  
  groupDat3 <- groupDat2[-set2$idR, ]
  set3 <- sample_n(groupDat3, sampSize)
  
  groupDat4 <- groupDat3[-set3$idR, ]
  set4 <- sample_n(groupDat4, sampSize)
  
  groupDat5 <- groupDat4[-set4$idR, ]
  set5 <- sample_n(groupDat5, sampSize)
  
  # Store sets (ungrouped)
  setList <- lapply(list(set1, set2, set3, set4, set5), ungroup)
  axList[[h]]$sets <- setList
  
  # List to store CV
  cvList <- list()
  
  for (i in 1:5) {
    
    # Stratified random training/test sets
    trainSet <- do.call(rbind, setList[-i])
    testSet <- setList[[i]]
    
    # OLS with training set
    my.reg <- lm(LR ~ LP - 1, trainSet)
    # Create adjusted residuals
    leverage <- influence(my.reg)$hat
    my.s.resid <- residuals(my.reg)/sqrt(1-leverage)
    my.s.resid <- my.s.resid - mean(my.s.resid)
    
    # List to store planes
    planeList <- list()
    
    # Bootstrapping on test set
    for (j in c("XY", "XZ", "YZ")) {
      
      # Subset plane
      dataSub <- testSet[testSet$Plane == j, ]
      
      # List to store errors
      epList <- list()
      epCount <- 1
      
      for (k in dataSub[ , "LP"]) {
        
        # Do bootstrap with 100 replications
        ep.draws <- replicate(n=100, the.replication.type(reg = my.reg, s = my.s.resid, x_Np1 = k))
        
        # Store in list
        epList[[epCount]] <- ep.draws
        epCount <- epCount + 1
        
      }
      
      # Store in list
      planeList[[j]] <- epList
      
      print(paste("CV", i, "Plane", j, "finished"))
      
    }
    
    # Store in list
    cvList[[i]] <- planeList
    
  }
  
  # Store in list
  axList[[h]]$cv <- cvList
  
  print(paste("Axon", h, "finished"))
  
}



saveRDS(axList, paste("projection_5CV_allAxons.rds", sep = ""))
```


##Estimate porcentual errors

```{r, fig.width=8, fig.height=8}
# Load 5-fold CV
axList <- readRDS("ProyeccionAnalysis/projection_5CV_allAxons.rds")

# Estimate and store in list
peList <- list()
peCount <- 1
lrLength <- list()

for (h in c("Type1", "Type2", "Type3", "Type4")) {
  
  for (i in 1:5) {
    
    for (j in c("XY", "XZ", "YZ")) {
      
      # Extract LR values
      lr <- axList[[h]]$sets[[i]]
      lr <- lr[lr$Plane == j, "LR"]
      
      # Extract LP values
      lp <- axList[[h]]$cv[[i]][[j]][[1]]
      
      for (k in 1:dim(lr)[1]) {
        
        # Divide errors by LR to get porcentual errors
        pe <- lp[k, ] / lr$LR[k]
        
        # Store in list
        peList[[peCount]] <- pe*100
        peCount <-peCount + 1
      }
      
    }
    
  }
  
  lrLength[h] <- dim(lr)[1]
  
}

# Store in data frame
lrLength <- unlist(lrLength)
peFrame <- data.frame(pe = unlist(peList),
                      neurType = rep(c("Type1", "Type2", "Type3", "Type4"), times = lrLength*5*3*100),
                      Plane = c(rep(rep(c("XY", "XZ", "YZ"), each = lrLength[1]*100), 5),
                                rep(rep(c("XY", "XZ", "YZ"), each = lrLength[2]*100), 5),
                                rep(rep(c("XY", "XZ", "YZ"), each = lrLength[3]*100), 5),
                                rep(rep(c("XY", "XZ", "YZ"), each = lrLength[4]*100), 5)))

peFrame$interact <- interaction(peFrame$neurType, peFrame$Plane, lex.order=T)
```


####Plot density and histogram

```{r}
# Plot
# peGroup <- peFrame %>% group_by(interact)
# peSample <- sample_n(peGroup, 100)

setwd("ProyeccionFigures/")
png(filename=paste("prederrorDistributions_5CV.png", sep=""), type="cairo",
    units="in", width=8, height=10, pointsize=12,
    res=300)

(p <- ggplot(peFrame, aes(x = `pe`, y = `interact`, fill = ifelse(..x.. >=-5 & ..x.. <=5, ">= -5", "< 5"))) +
  scale_fill_manual(values = c("gray70", "coral2"), name = NULL) +
  geom_density_ridges_gradient(alpha = 0.8) +
  theme_ridges() + 
  theme(legend.position = "none") + 
  xlab("% Error")) + 
  xlim(c(-40,40))

dev.off()

png(filename=paste("prederrorHistogram_5CV.png", sep=""), type="cairo",
    units="in", width=8, height=10, pointsize=12,
    res=300)

# Histogram
# nint 300 bootPI
# nint 500 bootPI_full
# nint 400 bootPI_full_replace
(h <- histogram( ~ pe | Plane + neurType, peFrame, nint = 300, type ="density", xlim = c(-50,50)))

dev.off()
```

#####Insets

```{r}
# Find mean bandwidth
bwFrame <- aggregate(pe ~ neurType + Plane, peFrame, function(x) density(x)$bw)

# Named color vector
colVec <- c("limegreen", "blue", "red")
names(colVec) <- c("XY", "XZ", "YZ")

setwd("ProyeccionFigures/")

for (i in unique(peFrame$neurType)) {
  
  png(filename=paste("insetProb_", i, ".png", sep=""), type="cairo",
    units="in", width=5, height=5, pointsize=12,
    res=300)

  par(mfrow = c(3,2),
      oma = c(5,4,0,0) + 0.1,
      mar = c(0,0,1,1) + 0.1)
  
  for (j in unique(peFrame$Plane)) {
    
    # Estimate density
    dens <- density(peFrame[peFrame$neurType == i & peFrame$Plane == j, "pe"], bw = mean(bwFrame$pe))

    # 5%
    x1 <- min(which(dens$x >= -5))  
    x2 <- max(which(dens$x <  5))
    
    plot(dens, ylim = c(0, 0.15), xlim = c(-20,20), zero.line = FALSE, main = "", ylab = "", xlab = "",  axes = FALSE)
    with(dens, polygon(x=c(x[c(x1,x1:x2,x2)]), y= c(0, y[x1:x2], 0), col=colVec[[j]], border = NA))
    lines(dens)
    
    # Add y axis
    axis(side = 2, at=seq(-0.05,0.15,0.05), labels = seq(-0.05,0.15,0.05), cex.axis = 1.2)

    # Different x axis for bottom plot
    if (j == "YZ") {
      axis(side = 1, at=seq(-30,30,10), labels = seq(-30,30,10), cex.axis = 1.2)    
    } else {
      axis(side = 1, at=seq(-30,30,10), labels = FALSE)    
    }
   
    
    # 10%
    x1 <- min(which(dens$x >= -10))  
    x2 <- max(which(dens$x <  10))
    
    plot(dens, ylim = c(0, 0.15), xlim = c(-20,20), zero.line = FALSE, main = "", ylab = "", xlab = "",  axes = FALSE)
    with(dens, polygon(x=c(x[c(x1,x1:x2,x2)]), y= c(0, y[x1:x2], 0), col=colVec[[j]], border = NA))
    lines(dens)
    
    # Add y axis
    axis(side = 2, at=seq(-0.05,0.15,0.05), labels = FALSE)

    # Different x axis for bottom plot
    if (j == "YZ") {
      axis(side = 1, at=seq(-30,30,10), labels = seq(-30,30,10), cex.axis = 1.2)    
    } else {
      axis(side = 1, at=seq(-30,30,10), labels = FALSE)    
    }
    
  }
  
  
  title(xlab = "% Error",
        ylab = "Density",
        outer = TRUE, line = 3,
        cex.lab = 1.5)
  
  dev.off()
  
}


```


####Errors vs cumul probability

RUN ONLY ONCE, AFTERWARDS LOAD THE FILE IN PLOT CHUNKS

```{r, fig.height=5, fig.width=20}
# Inverse quantile
quantInv <- function(distr, value) ecdf(distr)(value)

#ecdf plot example
# plot(ecdf(peFrame[peFrame$interact == "Type1.XY", "pe"]))
# lines(ecdf(peFrame[peFrame$interact == "Type1.XZ", "pe"]), col = "red")
# lines(ecdf(peFrame[peFrame$interact == "Type1.YZ", "pe"]), col = "blue")

probList <- list()
binProb <- 0.5
probSeq <- seq(binProb, 100, binProb)

for (i in unique(peFrame$interact)) {
  
  dataProb <- peFrame[peFrame$interact == i, "pe"]
  probVec <- numeric()
  
  for (j in probSeq) {
    
    errProb <- quantInv(dataProb, j) - quantInv(dataProb, -j)
    probVec <- c(probVec, errProb)
    
  }
  
  probList[[i]] <- probVec
  
}

setwd("ProyeccionAnalysis/")
saveRDS(probList, "cumProb_5CV.rds")
```

#####Plot

```{r, fig.height=5, fig.width=20}
# Plot with limited axis
setwd("ProyeccionFigures/")
probList <- readRDS("cumProb_5CV.rds")
png(filename=paste("errorProbabilityCum_bin", binProb, "_5CV.png", sep=""), type="cairo",
    units="in", width=20, height=5, pointsize=12,
    res=300)

par(mfrow = c(1,4))
hlist <- list()
axCount <- 1


for (i in seq(1,10,3)) {
  # Extract height for 5% error
  h5 <- sapply(probList[i:(i+2)], function(x) x[5/binProb])
  max5 <- max(h5)
  
  # Extract height for 10% error
  h10 <- sapply(probList[i:(i+2)], function(x) x[10/binProb])
  max10 <- max(h10)
  
  # Plot cumulative curves
  plot(probSeq, probList[[i]], type = "l", col = "limegreen", lwd = 1.5, ylim = c(0,1), xlim = c(0.5, 40),
       ylab = "Cumulative probability", xlab = "% Error", main = paste("Axon Type", axCount))
  lines(probSeq, probList[[i + 1]], col = "red", lwd = 1.5)
  lines(probSeq, probList[[i + 2]], col = "blue", lwd = 1.5)
  
  # Plot vertical segments at 5% and 10%
  segments(5, -0.1, 5, max5, lty = "dashed", col = "gray")
  segments(10, -0.1, 10, max10, lty = "dashed", col = "gray")
  
   # Plot horizontal segments at cumul prob for 5% and 10%
  transp <- 0.8
  segments(-2, h5[1], 5, h5[1], lty = "dashed", col = alpha("limegreen", transp))
  segments(-2, h5[2], 5, h5[2], lty = "dashed", col = alpha("red", transp))
  segments(-2, h5[3], 5, h5[3], lty = "dashed", col = alpha("blue", transp))
  
  segments(-2, h10[1], 10, h10[1], lty = "dashed", col = alpha("limegreen", transp))
  segments(-2, h10[2], 10, h10[2], lty = "dashed", col = alpha("red", transp))
  segments(-2, h10[3], 10, h10[3], lty = "dashed", col = alpha("blue", transp))
  
  # Add legend
  # legend("bottomright", c("XY", "XZ", "YZ"), lwd = 1.5, col = c("limegreen", "red", "blue"))
  
  # Store to check later
  hlist[[paste("Type", axCount, sep="")]] <- data.frame(h5 = h5, h10 = h10)
  axCount <- axCount + 1
}

dev.off()

# Plot with free axis
png(filename=paste("errorProbabilityCum_bin", binProb, "_freeAxis_5CV.png", sep=""), type="cairo",
    units="in", width=20, height=5, pointsize=12,
    res=300)

par(mfrow = c(1,4))
hlist <- list()
axCount <- 1


for (i in seq(1,10,3)) {
  # Extract height for 5% error
  h5 <- sapply(probList[i:(i+2)], function(x) x[5/binProb])
  max5 <- max(h5)
  
  # Extract height for 10% error
  h10 <- sapply(probList[i:(i+2)], function(x) x[10/binProb])
  max10 <- max(h10)
  
  # Plot cumulative curves
  plot(probSeq, probList[[i]], type = "l", col = "green", lwd = 1.5,
       ylab = "Cumulative probability", xlab = "% Error", main = paste("Axon Type", axCount))
  lines(probSeq, probList[[i + 1]], col = "red", lwd = 1.5)
  lines(probSeq, probList[[i + 2]], col = "blue", lwd = 1.5)
  
  # Plot vertical segments at 5% and 10%
  segments(5, -0.1, 5, max5, lty = "dashed", col = "gray")
  segments(10, -0.1, 10, max10, lty = "dashed", col = "gray")
  
   # Plot horizontal segments at cumul prob for 5% and 10%
  segments(-2, h5[1], 5, h5[1], lty = "dashed", col = "yellowgreen")
  segments(-2, h5[2], 5, h5[2], lty = "dashed", col = "tomato")
  segments(-2, h5[3], 5, h5[3], lty = "dashed", col = "lightskyblue")
  
  segments(-2, h10[1], 10, h10[1], lty = "dashed", col = "yellowgreen")
  segments(-2, h10[2], 10, h10[2], lty = "dashed", col = "tomato")
  segments(-2, h10[3], 10, h10[3], lty = "dashed", col = "lightskyblue")
  
  # Add legend
  legend("bottomright", c("XY", "XZ", "YZ"), lwd = 1.5, col = c("green", "red", "blue"))
  
  # Store to check later
  hlist[[paste("Type", axCount, sep="")]] <- data.frame(h5 = h5, h10 = h10)
  axCount <- axCount + 1
}

dev.off()
```


##Full figure

```{r, fig.width=20}
## show the regions that have been allocated to each plot
# layout.show(8)


#############
# Load data #
#############

probList <- readRDS("ProyeccionAnalysis/cumProb_5CV.rds")


##############
# MAIN PLOTS #
##############

setwd("ProyeccionFigures/")
png(filename=paste("errorProbabilityCum_bin", binProb, "_5CV_FULL.png", sep=""), type="cairo",
    units="in", width=20, height=5, pointsize=12,
    res=300)

# Graphical parameters
par(oma = c(1,4,1,1),
    mar = c(4,0,3,0) + 0.5,
    cex.axis = 1.5,
    cex.lab = 1.5,
    cex.main = 2)

layout(matrix(c(1,1,2,2,3,3,0,4,4,
                1,5,2,7,3,9,0,4,11,
                1,6,2,8,3,10,0,4,12,
                1,0,2,0,3,0,0,4,0), 4, 9, byrow = TRUE),
        widths = c(2.5,2.5,2.5,2.5,2.5,2.5,1,3,2.5),
        heights = c(1.6,1.6,1.6,1.2))
# layout.show(12)

# Plot loop
hlist <- list()
axCount <- 1


for (i in seq(1,10,3)) {
  
  # Extract height for 5% error
  h5 <- sapply(probList[i:(i+2)], function(x) x[5/binProb])
  
  # Extract height for 10% error
  h10 <- sapply(probList[i:(i+2)], function(x) x[10/binProb])
  
  # Plot cumulative error and Axis labels
  
  if (axCount == 1) {
    
    plot(probSeq, probList[[i]], type = "l", col = "limegreen", lwd = 1.5, ylim = c(0,1), xlim = c(0.5, 40),
         main = paste("Axon Type", axCount), yaxt = "n", xaxt = "n",  xlab = "", ylab = "")
    lines(probSeq, probList[[i + 1]], col = "red", lwd = 1.5)
    lines(probSeq, probList[[i + 2]], col = "blue", lwd = 1.5)
    title(ylab = "Cumulative probability",
          outer = TRUE, line=3)
    
    
  } else if (axCount == 2) {
    
    plot(probSeq, probList[[i]], type = "l", col = "limegreen", lwd = 1.5, ylim = c(0,1), xlim = c(0.5, 40),
         main = paste("Axon Type", axCount), yaxt = "n", xaxt = "n",  xlab = "", ylab = "")
    lines(probSeq, probList[[i + 1]], col = "red", lwd = 1.5)
    lines(probSeq, probList[[i + 2]], col = "blue", lwd = 1.5)
    title(xlab="% Error", line=3)
    
    
  } else if (axCount == 3) {
    
    plot(probSeq, probList[[i]], type = "l", col = "limegreen", lwd = 1.5, ylim = c(0,1), xlim = c(0.5, 40),
         main = paste("Axon Type", axCount), yaxt = "n", xaxt = "n",  xlab = "", ylab = "")
    lines(probSeq, probList[[i + 1]], col = "red", lwd = 1.5)
    lines(probSeq, probList[[i + 2]], col = "blue", lwd = 1.5)
    
  } else if (axCount == 4) {
    
    par(mar = c(4,5,3,0) + 0.5)
    plot(probSeq, probList[[i]], type = "l", col = "limegreen", lwd = 1.5, ylim = c(0,1), xlim = c(0.5, 40),
         main = paste("Axon Type", axCount), yaxt = "n", xaxt = "n",  xlab = "", ylab = "")
    lines(probSeq, probList[[i + 1]], col = "red", lwd = 1.5)
    lines(probSeq, probList[[i + 2]], col = "blue", lwd = 1.5)
    title(ylab = "Cumulative probability", xlab="% Error", line=3)
    
  }
 
  # Add x axis
  axis(side = 1, at=seq(-10,40,10), labels = seq(-10,40,10))
  
  # Different y axis for Type1 and Type4
  if (axCount == 1) {
    
    axis(side = 2, at=seq(-0.2,1,0.2), labels = seq(-0.2,1,0.2))
    
  } else if (axCount == 4) {
    
    axis(side = 1, at=seq(-10,40,10), labels = seq(-10,40,10))
    axis(side = 2, at=seq(-0.2,1,0.2), labels = seq(-0.2,1,0.2))
    
  }
  
  # Plot vertical segments at 5% and 10% for XY plane
  segments(5, -0.1, 5, h5[1], lty = "dashed", col = "darkgray")
  segments(10, -0.1, 10, h10[1], lty = "dashed", col = "darkgray")
  
   # Plot horizontal segments at cumul prob for 5% and 10% for XY plane
  transp <- 0.8
  segments(-2, h5[1], 5, h5[1], lty = "dashed", col = alpha("darkgray", transp))
  segments(-2, h10[1], 10, h10[1], lty = "dashed", col = alpha("darkgray", transp))
  
  # Store to check later
  hlist[[paste("Type", axCount, sep="")]] <- data.frame(h5 = h5, h10 = h10)
  axCount <- axCount + 1
}


###############
# PLOT INSETS #
###############

# Find mean bandwidth
bwFrame <- aggregate(pe ~ neurType + Plane, peFrame, function(x) density(x)$bw)

for (i in unique(peFrame$neurType)) {
  
  # Estimate density
  dens <- density(peFrame[peFrame$neurType == i & peFrame$Plane == "XY", "pe"], bw = mean(bwFrame$pe))

  # 5%
  x1 <- min(which(dens$x >= -5))  
  x2 <- max(which(dens$x <  5))
  
  par(mar = c(0.2,5,0.2,2))
  plot(dens, ylim = c(0, 0.15), xlim = c(-20,20), zero.line = FALSE, main = "", ylab = "", xlab = "",  axes = FALSE)
  with(dens, polygon(x=c(x[c(x1,x1:x2,x2)]), y= c(0, y[x1:x2], 0), col="limegreen", border = NA))
  lines(dens)
  
  # Add y axis
  axis(side = 2, at=seq(-0.05,0.15,0.05), labels = seq(-0.05,0.15,0.05), cex.axis = 1.2)
    
  
  # 10%
  x1 <- min(which(dens$x >= -10))  
  x2 <- max(which(dens$x <  10))
  
  plot(dens, ylim = c(0, 0.15), xlim = c(-20,20), zero.line = FALSE, main = "", ylab = "", xlab = "",  axes = FALSE)
  with(dens, polygon(x=c(x[c(x1,x1:x2,x2)]), y= c(0, y[x1:x2], 0), col="limegreen", border = NA))
  lines(dens)
  
  # Add y axis
  axis(side = 2, at=seq(-0.05,0.15,0.05), labels = seq(-0.05,0.15,0.05), cex.axis = 1.2)
  # Add x axis
  axis(side = 1, at=seq(-30,30,10), labels = seq(-30,30,10), cex.axis = 1.2)    
      
  
}

dev.off()
```