R scripts used to statistically assess differences in performance between methods for axon length measurement

Code/Statistical Analysis/EsferasAnalysis/desktop.ini

@@ -0,0 +1,5 @@
+[.ShellClassInfo]
+InfoTip=Esta carpeta se ha compartido online.
+IconFile=C:\Program Files\Google\Drive\googledrivesync.exe
+IconIndex=16
+IconResource=C:\Program Files\Google\Drive File Stream\52.0.6.0\GoogleDriveFS.exe,23

Code/Statistical Analysis/EsferasData/desktop.ini

@@ -0,0 +1,5 @@
+[.ShellClassInfo]
+InfoTip=Esta carpeta se ha compartido online.
+IconFile=C:\Program Files\Google\Drive\googledrivesync.exe
+IconIndex=16
+IconResource=C:\Program Files\Google\Drive File Stream\52.0.6.0\GoogleDriveFS.exe,23

Code/Statistical Analysis/EsferasFigures/desktop.ini

@@ -0,0 +1,5 @@
+[.ShellClassInfo]
+InfoTip=Esta carpeta se ha compartido online.
+IconFile=C:\Program Files\Google\Drive\googledrivesync.exe
+IconIndex=16
+IconResource=C:\Program Files\Google\Drive File Stream\52.0.6.0\GoogleDriveFS.exe,23

Code/Statistical Analysis/EstereoAnalysis/desktop.ini

@@ -0,0 +1,5 @@
+[.ShellClassInfo]
+InfoTip=Esta carpeta se ha compartido online.
+IconFile=C:\Program Files\Google\Drive\googledrivesync.exe
+IconIndex=16
+IconResource=C:\Program Files\Google\Drive File Stream\52.0.6.0\GoogleDriveFS.exe,23

Code/Statistical Analysis/EstereoAnalysis/stereo_boot10CVPI_LRLP_pred_Type1.rds

@@ -0,0 +1 @@
+/annex/objects/MD5-s26153526--62faab3218cbe1c6de8ceedc70a38e4c

Code/Statistical Analysis/EstereoAnalysis/stereo_boot10CVPI_LRLP_pred_Type2.rds

@@ -0,0 +1 @@
+/annex/objects/MD5-s44538132--23b8d63af4f3ce1f0671dd3e84dbff5e

Code/Statistical Analysis/dataforPepo.Rmd

@@ -0,0 +1,484 @@
+---
+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")
+```
+
+##Libraries and functions
+
+```{r}
+suppressMessages(library(R.matlab))  #import mat files
+```
+
+
+##PROJECTIONS
+
+###Cumulative probability
+
+####Format and save
+
+```{r}
+# Load rds
+probList <- readRDS("ProyeccionAnalysis/cumProb_5CV.rds")
+binProb <- 0.5
+probSeq <- seq(binProb, 100, binProb)
+
+# Store as long format data frame
+probFrame <- data.frame(axonType = rep(c("Type1", "Type2", "Type3", "Type4"), each = 3*200),
+                        Plane = rep(rep(c("XY", "XZ", "YZ"), each = 200),4),
+                        error = rep(probSeq, 12),
+                        cumulProb = unlist(probList))
+
+# Save as csv
+write.csv(probFrame, "dataPepo/projection_error_cumulativeProbability.csv", sep=",", row.names = FALSE)
+```
+
+###Error density
+
+####Estimation
+
+```{r}
+# Load 5-fold CV
+axList <- readRDS("ProyeccionAnalysis/modelBased_projections_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:length(lr)) {
+        
+        # Divide errors by LR to get porcentual errors
+        pe <- lp[k, ] / lr[k]
+        
+        # Store in list
+        peList[[peCount]] <- pe*100
+        peCount <-peCount + 1
+      }
+      
+    }
+    
+  }
+  
+  lrLength[h] <- length(lr)
+  
+}
+
+# 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)
+
+# Find mean bandwidth
+bwFrame <- aggregate(pe ~ neurType + Plane, peFrame, function(x) density(x)$bw)
+
+# Estimate density
+densList <- by(peFrame$pe, peFrame$interact, function(x) density(x, bw = mean(bwFrame$pe)))
+```
+
+####Format and save
+
+```{r}
+# Extract from list
+densXY <- lapply(densList, function(x) cbind(error = x$x, density = x$y))
+densFrame <- data.frame(axonType = rep(c("Type1", "Type2", "Type3", "Type4"), each = 512*3),
+                        Plane = rep(rep(c("XY", "XZ", "YZ"), each = 512), 4))
+
+# Store in long format data frame
+densFrame <- cbind(densFrame, do.call(rbind, densXY))
+
+# Save density as csv
+write.csv(densFrame, "dataPepo/projection_error_Densities.csv", sep=",", row.names = FALSE)
+
+# Save original data as csv
+oriFrame <- cbind(axonType = peFrame[, "neurType"], Plane = peFrame[, "Plane"], error = peFrame[, "pe"], data.frame(abserror = abs(peFrame$pe)))
+write.csv(oriFrame, "dataPepo/projection_error_OriginalData.csv", sep=",", row.names = FALSE)
+```
+
+###Mean errors
+
+####Estimate, format and save
+
+```{r}
+# Mean absolute error
+meanErr <- aggregate(abs(pe) ~ Plane + neurType, peFrame, mean)
+
+# SD absolute error
+sdErr <- aggregate(abs(pe) ~ Plane + neurType, peFrame, sd)
+
+# SE absolute error
+seErr <- aggregate(abs(pe) ~ Plane + neurType, peFrame, sjstats::se)
+
+# Common data frame
+meansdErr <- data.frame(axonType = meanErr$neurType, 
+                        Plane = meanErr$Plane, 
+                        mean = meanErr$`abs(pe)`, 
+                        sd = sdErr$`abs(pe)`,
+                        se = seErr$`abs(pe)`,
+                        ci95 = 1.96*seErr$`abs(pe)`)
+
+# Save as csv
+write.csv(meansdErr, "dataPepo/projection_error_MeanSD.csv", sep=",", row.names = FALSE)
+```
+
+###Five-number summary
+
+####Format and save
+```{r}
+# Store boxplot data
+boxData <- boxplot(abs(pe) ~ Plane + neurType, peFrame)
+fiveNum <- boxData$stats
+
+# Reformat as data frame
+fiveNum.frame <- data.frame(axonType = rep(c("Type1", "Type2", "Type3", "Type4"), each = 15),
+                            Plane = rep(rep(c("XY", "XZ", "YZ"), each = 5), 4),
+                            quantileType = rep(c("lowlim", "p25", "p50", "p75", "uplim"), 12),
+                            quantileDat = as.vector(fiveNum))
+
+# Save
+write.csv(fiveNum.frame, "dataPepo/projection_error_boxplot.csv", sep=",", row.names = FALSE)
+```
+
+
+##STEREO PLANES
+
+###Load data
+
+```{r}
+# Get file paths
+axonNames <- list.files(paste("EstereoDataPlanos/", sep=""), full.names=T)
+
+# Load matlab file
+axonFiles <- lapply(axonNames, function(x) readMat(x))
+names(axonFiles) <- c("Type1", "Type2", "Type3", "Type4")
+
+# Extract data
+# errorMatrix contains 4 arrays, one per neuron type
+# within each array, the dimensions corresponds to step:dist:neuron
+# this means that each array has as many matrices as neuron of a certain type
+dist         <- lapply(axonFiles, function(x) x$Distancias.Entre.Planos)[[1]]
+step         <- lapply(axonFiles, function(x) x$Step.Lengths)[[1]]
+errorMatrix  <- lapply(axonFiles, function(x) x$MATRIZ.ERROR.LENGTH)
+lpMatrix     <- lapply(axonFiles, function(x) x$MATRIZ.ESTIMATED.AXON.LENGTH)
+lrVals       <- lapply(axonFiles, function(x) x$AXON.REAL.LENGTH)
+
+# Get number of neurons per type
+numberTypes  <- unlist(lapply(errorMatrix, function(x) dim(x)[3]))
+
+# Vectorizing the arrays goes as follows: neuron, dist, step
+errorVector <- unlist(lapply(errorMatrix, function(x) as.vector(x)))
+lpVector <- unlist(lapply(lpMatrix, function(x) as.vector(x)))
+lrVector <- unlist(lapply(lrVals, function(x) as.vector(x)))
+
+# Store in data frames
+allData <- data.frame(id = rep(1:sum(numberTypes), each = 90),
+                      neurType = rep(c("Type1", "Type2", "Type3", "Type4"), times = 90*numberTypes),
+                      dist = rep(rep(dist, each = 9), sum(numberTypes)), 
+                      step = rep(rep(step, 10), sum(numberTypes)),
+                      error = abs(errorVector),
+                      error2 = errorVector,
+                      LRe = lpVector,
+                      LR = rep(lrVector, each = 90))
+allData$errorRaw <- allData$LR - allData$LRe
+allData$interact <- interaction(allData$step, allData$dist, lex.order = T)
+allData$interact2 <- interaction(allData$neurType, allData$step, allData$dist, lex.order = T)
+```
+
+###Mean errors
+
+####Estimate, format and save
+
+```{r}
+# Estimate mean residual error
+meanFrame <- aggregate(error ~ dist + step + neurType, allData, mean)
+n <- 200
+grid <- expand.grid(dist = seq(3, 30, length = n), step = seq(70, 150, length = n))
+
+# Smooth lm
+for (i in unique(meanFrame$neurType)) {
+  
+  axlm <- lm(error ~ dist * step, data = meanFrame[meanFrame$neurType == i, ])
+  print(summary(axlm))
+  axfit <- cbind(grid, prediction = predict(axlm, grid))
+  
+  # Save as csv
+  write.csv(axfit, paste("dataPepo/stereoPlanes_meanError_smooth_", i, ".csv", sep=""))
+  
+}
+```
+
+###Error density
+
+####Estimation
+
+```{r}
+# Find mean bandwidth
+bwFrame <- aggregate(error2 ~ neurType + dist + step, allData, function(x) density(x)$bw)
+
+# Estimate density using mean bw per axonType
+densList <- list()
+
+for (i in c("Type1", "Type2", "Type3", "Type4")) {
+  
+  axData <- allData[allData$neurType == i, ]
+  densList[[i]] <- by(axData$error2, axData$interact, function(x) density(x, 
+                                                                           bw = ceiling(mean(bwFrame[bwFrame$neurType == i, "error2"]))))
+  
+}
+
+# Restore format
+densList <- unlist(densList, recursive = FALSE)
+```
+
+
+####Format and save
+
+```{r}
+# Extract from list
+densXY <- lapply(densList, function(x) cbind(error = x$x, density = x$y))
+densFrame <- data.frame(axonType = rep(c("Type1", "Type2", "Type3", "Type4"), each = 512*90),
+                        step = rep(rep(unique(allData$step), each = 512*10), 4), 
+                        step = rep(rep(unique(allData$dist), each = 512), 4*9))
+
+# Store in long format data frame
+densFrame <- cbind(densFrame, do.call(rbind, densXY))
+
+# Save density as csv
+write.csv(densFrame, "dataPepo/stereoPlanes_error_Densities.csv", sep=",", row.names = FALSE)
+
+# Save original data as csv
+oriFrame <- cbind(axonType = allData[, "neurType"], dist = allData[, "dist"], step = allData[, "step"],
+                  error = allData[, "error2"], abserror = allData[, "error"])
+write.csv(oriFrame, "dataPepo/stereoPlanes_error_OriginalData.csv", sep=",", row.names = FALSE)
+```
+
+###Cumul probability
+
+####Estimate, format and save
+
+```{r}
+# Reformat as dataframe
+binProb <- 0.5
+probSeq <- seq(binProb, 100, binProb)
+axProb <- readRDS("EstereoAnalysis/errorProbs_RMSE.rds")
+probFrame <- data.frame(error = rep(probSeq, 90*4), 
+                        neurType = rep(c("Type1", "Type2", "Type3", "Type4"), each = length(probSeq)*90),
+                        dist = rep(rep(unique(allData$dist), each = length(probSeq)*9), 4),
+                        step = rep(rep(unique(allData$step), each = length(probSeq)), 10*4),
+                        prob = unlist(axProb))
+
+
+n <- 200
+grid <- expand.grid(dist = seq(3, 30, length = n), step = seq(70, 150, length = n))
+
+# Smooth lm
+for (h in c(5, 10)) {
+  
+  for (i in unique(probFrame$neurType)) {
+    
+    axlm <- lm(prob ~ dist * step, data = probFrame[probFrame$error == h & probFrame$neurType == i, ])
+    print(summary(axlm))
+    axfit <- cbind(grid, prediction = predict(axlm, grid))
+    
+    # Save as csv
+    write.csv(axfit, paste("dataPepo/stereoPlanes_cumulProb_error_", h, "_smooth_", i, ".csv", sep=""))
+    
+  }
+  
+  
+}
+
+```
+
+##STEREO SPHERES
+
+###Load data
+
+```{r}
+# Get file paths
+axonNames <- list.files(paste("EsferasData/", sep=""), full.names=T)
+
+# Load matlab file
+axonFiles <- lapply(axonNames, function(x) readMat(x))
+names(axonFiles) <- c("Type1", "Type2", "Type3_1", "Type3_2", "Type4")
+
+# Extract data
+# errorMatrix contains 4 arrays, one per neuron type
+# within each array, the dimensions corresponds to step:diams:neuron
+# this means that each array has as many matrices as neuron of a certain type
+diams         <- lapply(axonFiles, function(x) x$Diams.Sonda)[[1]]
+step         <- lapply(axonFiles, function(x) x$Step.Lengths)[[1]]
+errorMatrix  <- lapply(axonFiles, function(x) x$MATRIZ.ERROR.LENGTH)
+lpMatrix     <- lapply(axonFiles, function(x) x$ESTIMATED.AXON.LENGTH)
+lrVals       <- lapply(axonFiles, function(x) x$AXON.REAL.LENGTH)
+
+# Get number of neurons per type
+numberTypes  <- unlist(lapply(errorMatrix, function(x) dim(x)[3]))
+
+# Vectorizing the arrays goes as follows: neuron, diams, step
+errorVector <- unlist(lapply(errorMatrix, function(x) as.vector(x)))
+lpVector <- unlist(lapply(lpMatrix, function(x) as.vector(x)))
+lrVector <- unlist(lapply(lrVals, function(x) as.vector(x)))
+
+# Store in data frames
+allData <- data.frame(id = rep(1:sum(numberTypes), each = 81),
+                      diams = rep(rep(diams, each = 9), sum(numberTypes)), 
+                      step = rep(rep(step, 9), sum(numberTypes)),
+                      error = abs(errorVector),
+                      error2 = errorVector,
+                      LRe = lpVector,
+                      LR = rep(lrVector, each = 81))
+
+# Remove all LR = 0 (duplicated Type3)
+allData <- allData[allData$LR != 0, ]
+
+# Get number of neurons per type
+numberTypes  <- numberTypes[-3]
+names(numberTypes) <- c("Type1", "Type2", "Type3", "Type4")
+
+# Reassign id and neurType
+allData$id <- rep(1:sum(numberTypes), each = 81)
+allData$neurType <- rep(c("Type1", "Type2", "Type3", "Type4"), times = 81*numberTypes)
+
+allData$errorRaw <- allData$LR - allData$LRe
+allData$interact <- interaction(allData$step, allData$diams, lex.order = T)
+allData$interact2 <- interaction(allData$neurType, allData$step, allData$diams, lex.order = T)
+```
+
+###Mean errors
+
+####Estimate, format and save
+
+```{r}
+# Estimate mean residual error
+meanFrame <- aggregate(error ~ diams + step + neurType, allData, mean)
+n <- 200
+grid <- expand.grid(diams = seq(10, 50, length = n), step = seq(70, 150, length = n))
+
+# Smooth lm
+for (i in unique(meanFrame$neurType)) {
+  
+  axlm <- lm(error ~ diams * step, data = meanFrame[meanFrame$neurType == i, ])
+  print(summary(axlm))
+  axfit <- cbind(grid, prediction = predict(axlm, grid))
+  
+  # Save as csv
+  write.csv(axfit, paste("dataPepo/stereoSpheres_meanError_smooth_", i, ".csv", sep=""))
+  
+}
+```
+
+###Error density
+
+####Estimation
+
+```{r}
+# Find mean bandwidth
+bwFrame <- aggregate(error2 ~ neurType + diams + step, allData, function(x) density(x)$bw)
+
+# Estimate density using mean bw per axonType
+densList <- list()
+
+for (i in c("Type1", "Type2", "Type3", "Type4")) {
+  
+  axData <- allData[allData$neurType == i, ]
+  densList[[i]] <- by(axData$error2, axData$interact, function(x) density(x, 
+                                                                           bw = ceiling(mean(bwFrame[bwFrame$neurType == i, "error2"]))))
+  
+}
+
+# Restore format
+densList <- unlist(densList, recursive = FALSE)
+```
+
+
+####Format and save
+
+```{r}
+# Extract from list
+densXY <- lapply(densList, function(x) cbind(error = x$x, density = x$y))
+densFrame <- data.frame(axonType = rep(c("Type1", "Type2", "Type3", "Type4"), each = 512*81),
+                        step = rep(rep(unique(allData$step), each = 512*9), 4), 
+                        step = rep(rep(unique(allData$diams), each = 512), 4*9))
+
+# Store in long format data frame
+densFrame <- cbind(densFrame, do.call(rbind, densXY))
+
+# Save density as csv
+write.csv(densFrame, "dataPepo/stereoSpheres_error_Densities.csv", sep=",", row.names = FALSE)
+
+# Save original data as csv
+oriFrame <- cbind(axonType = allData[, "neurType"], diams = allData[, "diams"], step = allData[, "step"],
+                  error = allData[, "error2"], abserror = allData[, "error"])
+write.csv(oriFrame, "dataPepo/stereoSpheres_error_OriginalData.csv", sep=",", row.names = FALSE)
+```
+
+###Cumul probability
+
+####Estimate, format and save
+
+```{r}
+# Reformat as dataframe
+binProb <- 0.5
+probSeq <- seq(binProb, 100, binProb)
+axProb <- readRDS("EsferasAnalysis/errorProbs_RMSE.rds")
+probFrame <- data.frame(error = rep(probSeq, 81*4), 
+                        neurType = rep(c("Type1", "Type2", "Type3", "Type4"), each = length(probSeq)*81),
+                        diams = rep(rep(unique(allData$diams), each = length(probSeq)*9), 4),
+                        step = rep(rep(unique(allData$step), each = length(probSeq)), 9*4),
+                        prob = unlist(axProb))
+
+
+n <- 200
+grid <- expand.grid(diams = seq(10, 50, length = n), step = seq(70, 150, length = n))
+
+# Smooth lm
+for (h in c(5, 10)) {
+  
+  for (i in unique(probFrame$neurType)) {
+    
+    axlm <- lm(prob ~ diams * step, data = probFrame[probFrame$error == h & probFrame$neurType == i, ])
+    axfit <- cbind(grid, prediction = predict(axlm, grid))
+    
+    # Save as csv
+    write.csv(axfit, paste("dataPepo/stereoSpheres_cumulProb_error_", h, "_smooth_", i, ".csv", sep=""))
+    
+  }
+  
+  
+}
+
+```

Code/Statistical Analysis/desktop.ini

@@ -0,0 +1,5 @@
+[.ShellClassInfo]
+InfoTip=Esta carpeta se ha compartido online.
+IconFile=C:\Program Files\Google\Drive\googledrivesync.exe
+IconIndex=16
+IconResource=C:\Program Files\Google\Drive File Stream\52.0.6.0\GoogleDriveFS.exe,23

Code/Statistical Analysis/modelBased_projections_errors_5foldCV.Rmd

@@ -0,0 +1,843 @@
+---
+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()
+```
+
+

Code/Statistical Analysis/modelBased_projections_errors_CV.Rmd

@@ -0,0 +1,846 @@
+---
+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)
+# 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
+
+This code is meant to be run independently for every neurType and in parallel. RUN AT CREMOTO
+
+```{r}
+# Ensure reproducibility
+set.seed(123456)
+
+# Subset axon type
+dataCV <- allData[allData$neurType == "Type1", ]
+
+# List to store CV repetitions
+cvList <- list ()
+cvRand <- list()
+
+for (i in 1:10) {
+  
+  # Stratified random training/test sets
+  stratDat <- stratified(dataCV, "Plane", 0.7, bothSets = TRUE, replace = FALSE)
+  trainSet <- stratDat$SET1
+  testSet <- stratDat$SET2
+  # Store train and test set
+  cvRand[[i]] <- list(trainSet, testSet)
+  
+  # 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()
+  planeRand <- list()
+  # Assign test data to bootstrapping
+  dataToBoot <- testSet
+    
+  # Bootstrapping on test set
+    for (j in c("XY", "XZ", "YZ")) {
+      
+      # Subset plane
+      dataSub <- dataToBoot[dataToBoot$Plane == j, ]
+      
+      # List to store errors
+      epList <- list()
+      epCount <- 1
+      
+      for (k in dataSub[ , "LP"]) {
+        
+        # Predict LR
+        y.p <- coef(my.reg)*k
+        
+        # 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
+        
+        print(epCount)
+        
+      }
+      
+      # Store in list
+      planeList[[j]] <- epList
+      
+      print(paste("CV", i, "Plane", j, "finished"))
+      
+    }
+    
+    # Store in list
+    cvList[[i]] <- planeList
+    print(paste("CV", i, "finished"))
+    
+  }
+  
+  # names(axList) <- c("Type1", "Type2", "Type3", "Type4")
+  # names(axRand) <- c("Type1", "Type2", "Type3", "Type4")
+  # 
+
+
+
+# saveRDS(list(distList, distRand), paste("stereo_bootPI_full_Type", axCount, "_boot20.rds", sep = ""))
+```
+
+
+####Check convergence 
+
+```{r}
+# Load object
+bootCheck <- readRDS("bootPI_full.rds")
+# First list is prediction error, second is random subset
+axList <- bootCheck[[1]]
+axRand <- bootCheck[[2]]
+# One example
+
+# OLS
+my.reg <- lm(LR ~ LP:neurType - 1, allData)
+
+# Estimate and store in list
+lrpList <- list()
+lrpCount <- 1
+axCount <- 1
+
+for (i in c("Type1", "Type2", "Type3", "Type4")) {
+
+  for (j in c("XY", "XZ", "YZ")) {
+    
+     # Extract LP and LR values
+    lp <- axRand[[i]][[j]]$LP
+    lr <- axRand[[i]][[j]]$LR
+    
+    for (k in 1:length(lp)) {
+      
+      # Predicted LR
+      lr.p <- coef(my.reg)[axCount]*lp[k]
+      
+      # Extract absolute errors
+      ae <- axList[[i]][[j]][[k]]
+      
+      # Empty vectors for storage
+      lrp.Low <- numeric()
+      lrp.Up <- numeric()
+      
+      # Estimate PI95 for boot n=10...100
+      for (l in 4:100) {
+        
+        lrp <- lr.p + quantile(ae[1:l], probs = c(0.025,0.975))
+        
+        # Store in vectors
+        lrp.Low <- c(lrp.Low, lrp[1])
+        lrp.Up <- c(lrp.Up, lrp[2])
+        
+      }
+      
+      # Store in data frame
+      lrpFrame <- data.frame(nboot = 4:100, lowLim = lrp.Low, upLim = lrp.Up, LR = rep(lr[k], length(4:100)))
+      
+      # Store in list
+      lrpList[[lrpCount]] <- lrpFrame
+      lrpCount <- lrpCount + 1
+    }
+    
+  }
+  
+  axCount <- axCount + 1
+  
+}
+
+```
+
+#####Reformat
+
+```{r}
+# 1/3 random observations per plane, no replacement
+numExtract <- floor(1/3*numberTypes)
+
+checkFrame <- do.call(rbind, lrpList)
+checkFrame$neurType <- rep(c("Type1", "Type2", "Type3", "Type4"), times = numExtract*(length(4:100))*3)
+checkFrame$Plane <- c(rep(c("XY", "XZ", "YZ"), each = numExtract[1]*length(4:100)),
+                     rep(c("XY", "XZ", "YZ"), each = numExtract[2]*length(4:100)),
+                     rep(c("XY", "XZ", "YZ"), each = numExtract[3]*length(4:100)),
+                     rep(c("XY", "XZ", "YZ"), each = numExtract[4]*length(4:100)))
+
+
+lowMean <- aggregate(lowLim ~ neurType + Plane + nboot, checkFrame, mean)
+upMean <- aggregate(upLim ~ neurType + Plane + nboot, checkFrame, mean)
+LRmean <- aggregate(LR ~ neurType + Plane + nboot, checkFrame, mean)
+
+cumError <- with(lowMean[lowMean$neurType == "Type1" & lowMean$Plane == "XY", ], 
+             abs((lowLim[-1] - lowLim[-length(lowLim)])/lowLim[-length(lowLim)])*100)
+```
+
+#####Plot
+
+```{r, fig.height=20, fig.width=30}
+# Separate plots
+png(filename=paste("convergence_bootPI_separate.png", sep=""), type="cairo",
+    units="in", width=30, height=20, pointsize=12,
+    res=300)
+
+par(mfrow=c(4,6), cex = 1.1)
+
+for (i in c("Type1", "Type2", "Type3", "Type4")) {
+  
+  for (j in c("XY", "XZ", "YZ")) {
+    
+    plot(lowLim ~ nboot, lowMean[lowMean$neurType == i & lowMean$Plane == j, ], type = "l", lwd = 2)
+    plot(upLim ~ nboot, upMean[upMean$neurType == i & upMean$Plane == j, ], type = "l", lwd = 2)
+    
+  }
+  
+}
+
+dev.off()
+
+# Same plot
+png(filename=paste("convergence_bootPI_together.png", sep=""), type="cairo",
+    units="in", width=15, height=20, pointsize=12,
+    res=300)
+
+par(mfrow=c(4,3), cex = 1.1)
+
+for (i in c("Type1", "Type2", "Type3", "Type4")) {
+  
+  for (j in c("XY", "XZ", "YZ")) {
+    
+    minY <- min(lowMean[lowMean$neurType == i & lowMean$Plane == j, "lowLim"])
+    maxY <- max(upMean[upMean$neurType == i & upMean$Plane == j, "upLim"])
+    plot(lowLim ~ nboot, lowMean[lowMean$neurType == i & lowMean$Plane == j, ], type = "l", lwd = 2, ylim = c(minY, maxY))
+    lines(upLim ~ nboot, upMean[upMean$neurType == i & upMean$Plane == j, ], type = "l", lwd = 2)
+    # lines(LR ~ nboot, LRmean[LRmean$neurType == i & LRmean$Plane == j, ], type = "l", lwd = 2, col = "red")
+    
+  }
+  
+}
+
+dev.off()
+
+# Diff plot
+png(filename=paste("convergence_bootPI_Diff.png", sep=""), type="cairo",
+    units="in", width=15, height=20, pointsize=12,
+    res=300)
+
+par(mfrow=c(4,3), cex = 1.1)
+
+for (i in c("Type1", "Type2", "Type3", "Type4")) {
+  
+  for (j in c("XY", "XZ", "YZ")) {
+    
+    minY <- min(lowMean[lowMean$neurType == i & lowMean$Plane == j, "lowLim"])
+    maxY <- max(upMean[upMean$neurType == i & upMean$Plane == j, "upLim"])
+    diflu <- upMean[upMean$neurType == i & upMean$Plane == j, "upLim"] - lowMean[lowMean$neurType == i & lowMean$Plane == j, "lowLim"]
+    plot(4:100, diflu, type = "l", col = "red", lwd = 2)
+    
+  }
+  
+}
+
+dev.off()
+
+# Same plot random intervals
+
+for (m in 1:4) {
+  
+  png(filename=paste("convergence_bootPI_together_example_rand", m, ".png", sep=""), type="cairo",
+      units="in", width=15, height=20, pointsize=12,
+      res=300)
+  
+  par(mfrow=c(4,3), cex = 1.1)
+  
+  for (i in c("Type1", "Type2", "Type3", "Type4")) {
+    
+    for (j in c("XY", "XZ", "YZ")) {
+      
+      subData <- checkFrame[checkFrame$neurType == i & checkFrame$Plane == j, ]
+      randLR <- sample(subData[, "LR"], 1)
+      # randLR <- sample(subData[which(subData$lowLim > subData$LR), "LR"], 1)
+      minY <- min(subData[subData$LR == randLR, "lowLim"])
+      maxY <- max(subData[subData$LR == randLR, "upLim"])
+      plot(lowLim ~ nboot, subData[subData$LR == randLR, ], type = "l", lwd = 2, ylim = c(minY, maxY))
+      lines(upLim ~ nboot, subData[subData$LR == randLR, ], type = "l", lwd = 2)
+      lines(LR ~ nboot, subData[subData$LR == randLR, ], type = "l", lwd = 2, col = "red")
+      
+    }
+    
+  }
+  
+  dev.off()
+  
+}
+
+# Separate with lines
+png(filename=paste("convergence_bootPI_separate_FULL.png", sep=""), type="cairo",
+    units="in", width=30, height=20, pointsize=12,
+    res=300)
+
+par(mfrow=c(4,6), cex = 1.1)
+
+for (i in c("Type1", "Type2", "Type3", "Type4")) {
+  
+  for (j in c("XY", "XZ", "YZ")) {
+    
+    plot(lowLim ~ nboot, checkFrame[checkFrame$neurType == i & checkFrame$Plane == j, ], type = "l", col = "lightgray")
+    lines(lowLim ~ nboot, lowMean[lowMean$neurType == i & lowMean$Plane == j, ], type = "l", lwd = 2)
+    
+    plot(upLim ~ nboot, checkFrame[checkFrame$neurType == i & checkFrame$Plane == j, ], type = "l", col = "lightgray")
+    lines(upLim ~ nboot, upMean[upMean$neurType == i & upMean$Plane == j, ], type = "l", lwd = 2)
+    
+  }
+  
+}
+
+dev.off()
+
+
+```
+
+
+####Estimate porcentual errors
+
+```{r, fig.width=8, fig.height=8}
+# # Load objects
+cvPI <- list()
+
+for (i in c("Type1", "Type2", "Type3", "Type4")) {
+  # Whithin each Type, element 1 are erros and element 2 are random samples, training and test
+  cvPI[[i]] <- readRDS(paste("ProyeccionAnalysis/proyec_boot10CVPI_", i, ".rds", sep = ""))
+}
+
+# Estimate and store in list
+peList <- list()
+peCount <- 1
+lrLength <- list()
+
+for (h in c("Type1", "Type2", "Type3", "Type4")) {
+  
+  for (i in 1:10) {
+    
+    for (j in c("XY", "XZ", "YZ")) {
+      
+      # Extract LR values
+      lr <- cvPI[[h]][[2]][[i]][[2]]
+      lr <- lr[lr$Plane == j, "LR"]
+      
+      for (k in 1:length(lr)) {
+        
+        # Divide absolute errors by LR to get porcentual errors
+        pe <- cvPI[[h]][[1]][[i]][[j]][[k]] / lr[k]
+        
+        # Store in list
+        peList[[peCount]] <- pe*100
+        peCount <-peCount + 1
+      }
+      
+    }
+    
+  }
+  
+  lrLength[h] <- length(lr)
+  
+}
+
+lrLength <- unlist(lrLength)
+peFrame <- data.frame(pe = unlist(peList),
+                      neurType = rep(c("Type1", "Type2", "Type3", "Type4"), times = lrLength*10*3*100),
+                      Plane = c(rep(rep(c("XY", "XZ", "YZ"), each = lrLength[1]*100), 10),
+                                rep(rep(c("XY", "XZ", "YZ"), each = lrLength[2]*100), 10),
+                                rep(rep(c("XY", "XZ", "YZ"), each = lrLength[3]*100), 10),
+                                rep(rep(c("XY", "XZ", "YZ"), each = lrLength[4]*100), 10)))
+
+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("ProyeccionAnalysis/")
+png(filename=paste("prederrorDistributions_ols_CV.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_ols_CV.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()
+```
+
+
+
+#####Plot errors vs cumul probability
+
+```{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
+  
+}
+
+# Plot with limited axis
+setwd("ProyeccionAnalysis/")
+png(filename=paste("errorProbabilityCum_bin", binProb, "_CV.png", sep=""), type="cairo",
+    units="in", width=20, height=5, pointsize=12,
+    res=300)
+
+par(mfrow = c(1,4))
+axCount <- 1
+
+for (i in seq(1,10,3)) {
+  plot(probSeq, probList[[i]], type = "l", col = "green", 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)
+  abline(v = 5, lty = "dashed", col = "gray")
+  abline(h = 0.95, lty = "dashed", col = "gray" )
+  legend("bottomright", c("XY", "XZ", "YZ"), lwd = 1.5, col = c("green", "red", "blue"))
+  axCount <- axCount + 1
+}
+
+dev.off()
+
+# Plot with free axis
+png(filename=paste("errorProbabilityCum_bin", binProb, "_freeAxis_CV.png", sep=""), type="cairo",
+    units="in", width=20, height=5, pointsize=12,
+    res=300)
+
+par(mfrow = c(1,4))
+axCount <- 1
+
+for (i in seq(1,10,3)) {
+  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)
+  abline(v = 5, lty = "dashed", col = "gray")
+  abline(h = 0.95, lty = "dashed", col = "gray" )
+  legend("bottomright", c("XY", "XZ", "YZ"), lwd = 1.5, col = c("green", "red", "blue"))
+  axCount <- axCount + 1
+}
+
+dev.off()
+```
+
+#####Random sample
+
+```{r}
+size <- 0.05
+peRand <- peFrame %>% group_by(interact)
+peSample <- sample_frac(peRand, size)
+```
+
+
+######Plot density and histogram
+
+```{r}
+setwd("ProyeccionAnalysis/")
+png(filename=paste("prederrorDistributions_ols_CV_rand", size, ".png", sep=""), type="cairo",
+    units="in", width=8, height=10, pointsize=12,
+    res=300)
+
+(p <- ggplot(peSample, 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_ols_CV_rand", size, ".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, peSample, nint = 300, type ="density", xlim = c(-50,50)))
+
+dev.off()
+```
+
+
+######Plot errors vs cumul probability
+
+```{r, fig.height=5, fig.width=20}
+# Inverse quantile
+quantInv <- function(distr, value) ecdf(distr)(value)
+
+peSample <- as.data.frame(peSample)
+probList <- list()
+binProb <- 0.5
+probSeq <- seq(binProb, 100, binProb)
+
+for (i in unique(peSample$interact)) {
+  
+  dataProb <- peSample[peSample$interact == i, "pe"]
+  probVec <- numeric()
+  
+  for (j in probSeq) {
+    
+    errProb <- quantInv(dataProb, j) - quantInv(dataProb, -j)
+    probVec <- c(probVec, errProb)
+    
+  }
+  
+  probList[[i]] <- probVec
+  
+}
+
+# Plot with limited axis
+setwd("ProyeccionAnalysis/")
+png(filename=paste("errorProbabilityCum_bin", binProb, "_CV_rand", size, ".png", sep=""), type="cairo",
+    units="in", width=20, height=5, pointsize=12,
+    res=300)
+
+par(mfrow = c(1,4))
+axCount <- 1
+
+for (i in seq(1,10,3)) {
+  plot(probSeq, probList[[i]], type = "l", col = "green", 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)
+  abline(v = 5, lty = "dashed", col = "gray")
+  abline(h = 0.95, lty = "dashed", col = "gray" )
+  legend("bottomright", c("XY", "XZ", "YZ"), lwd = 1.5, col = c("green", "red", "blue"))
+  axCount <- axCount + 1
+}
+
+dev.off()
+
+# Plot with free axis
+png(filename=paste("errorProbabilityCum_bin", binProb, "_freeAxis_CV_rand", size, ".png", sep=""), type="cairo",
+    units="in", width=20, height=5, pointsize=12,
+    res=300)
+
+par(mfrow = c(1,4))
+axCount <- 1
+
+for (i in seq(1,10,3)) {
+  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)
+  abline(v = 5, lty = "dashed", col = "gray")
+  abline(h = 0.95, lty = "dashed", col = "gray" )
+  legend("bottomright", c("XY", "XZ", "YZ"), lwd = 1.5, col = c("green", "red", "blue"))
+  axCount <- axCount + 1
+}
+
+dev.off()
+```
+
+####Estimate mean porcentual errors
+
+UNFINISHED
+
+```{r, fig.width=8, fig.height=8}
+# mm <- peFrame[peFrame$Plane == "XY", "pe"]
+# mm.mat <- matrix(mm, 100, 270)
+# meanMat <- apply(mm.mat, 1, mean)
+
+```

Code/Statistical Analysis/modelFree_stereology_MSE.Rmd

@@ -0,0 +1,772 @@
+---
+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(ggpubr)
+library(ggridges)
+library(mltools)
+library(data.table)
+library(caret)
+library(interactions)
+library(data.table)
+library(wesanderson)
+
+options(digits = 10)
+
+# https://rstudio-pubs-static.s3.amazonaws.com/297778_5fce298898d64c81a4127cf811a9d486.html
+
+
+###############################
+# Add continuous color legend #
+###############################
+
+legend.col <- function(col, lev){
+  
+  opar <- par
+  
+  n <- length(col)
+  
+  bx <- par("usr")
+  
+  box.cx <- c(bx[2] + (bx[2] - bx[1]) / 1000,
+              bx[2] + (bx[2] - bx[1]) / 1000 + (bx[2] - bx[1]) / 50)
+  box.cy <- c(bx[3], bx[3])
+  box.sy <- (bx[4] - bx[3]) / n
+  
+  xx <- rep(box.cx, each = 2)
+  
+  par(xpd = TRUE)
+  for(i in 1:n){
+    
+    yy <- c(box.cy[1] + (box.sy * (i - 1)),
+            box.cy[1] + (box.sy * (i)),
+            box.cy[1] + (box.sy * (i)),
+            box.cy[1] + (box.sy * (i - 1)))
+    polygon(xx, yy, col = col[i], border = col[i])
+    
+  }
+  par(new = TRUE)
+  plot(0, 0, type = "n",
+       ylim = c(min(lev), max(lev)),
+       yaxt = "n", ylab = "",
+       xaxt = "n", xlab = "",
+       frame.plot = FALSE)
+  axis(side = 4, las = 2, at = lev, tick = FALSE, line = .15)
+  par <- opar
+}
+```
+
+
+##Load data
+
+```{r}
+# Get file paths
+axonNames <- list.files(paste("EstereoDataPlanos/", sep=""), full.names=T)
+
+# Load matlab file
+axonFiles <- lapply(axonNames, function(x) readMat(x))
+names(axonFiles) <- c("Type1", "Type2", "Type3", "Type4")
+
+# Extract data
+# errorMatrix contains 4 arrays, one per neuron type
+# within each array, the dimensions corresponds to step:dist:neuron
+# this means that each array has as many matrices as neuron of a certain type
+dist         <- lapply(axonFiles, function(x) x$Distancias.Entre.Planos)[[1]]
+step         <- lapply(axonFiles, function(x) x$Step.Lengths)[[1]]
+errorMatrix  <- lapply(axonFiles, function(x) x$MATRIZ.ERROR.LENGTH)
+lpMatrix     <- lapply(axonFiles, function(x) x$MATRIZ.ESTIMATED.AXON.LENGTH)
+lrVals       <- lapply(axonFiles, function(x) x$AXON.REAL.LENGTH)
+
+# Get number of neurons per type
+numberTypes  <- unlist(lapply(errorMatrix, function(x) dim(x)[3]))
+
+# Vectorizing the arrays goes as follows: neuron, dist, step
+errorVector <- unlist(lapply(errorMatrix, function(x) as.vector(x)))
+lpVector <- unlist(lapply(lpMatrix, function(x) as.vector(x)))
+lrVector <- unlist(lapply(lrVals, function(x) as.vector(x)))
+
+# Store in data frames
+allData <- data.frame(id = rep(1:sum(numberTypes), each = 90),
+                      neurType = rep(c("Type1", "Type2", "Type3", "Type4"), times = 90*numberTypes),
+                      dist = rep(rep(dist, each = 9), sum(numberTypes)), 
+                      step = rep(rep(step, 10), sum(numberTypes)),
+                      error = abs(errorVector),
+                      error2 = errorVector,
+                      LRe = lpVector,
+                      LR = rep(lrVector, each = 90))
+allData$errorRaw <- allData$LR - allData$LRe
+allData$interact <- interaction(allData$step, allData$dist, lex.order = T)
+allData$interact2 <- interaction(allData$neurType, allData$step, allData$dist, lex.order = T)
+```
+
+
+## Estimate MSE
+
+```{r}
+mseDistep <- by(allData, allData$interact2, function(x) sqrt(mse(x$LRe, x$LR)))
+```
+
+
+```{r, fig.width=15}
+# mseDist <- by(allData, allData$dist, function(x) sqrt(mse(x$LRe, x$LR)))
+# mseStep <- by(allData, allData$step, function(x) sqrt(mse(x$LRe, x$LR)))
+# mseDistep <- by(allData, allData$interact2, function(x) sqrt(mse(x$LRe, x$LR)))
+# 
+# par(mfrow=c(1,3))
+# 
+# plot(seq(3,30,3), mseDist, 
+#      xaxt = "n", xlab = "Dist", ylab = "RMSE", 
+#      pch = 21, bg = "firebrick", col = NA, 
+#      main = "Test Error ~ Distancia",
+#      cex = 1.5, cex.axis = 1.5, cex.main = 2, cex.lab = 1.5)
+# lines(seq(3,30,3), mseDist)
+# axis(side=1, at=seq(3,30,3), cex.axis = 1.5)
+# grid()
+# 
+# plot(seq(70,150,10), mseStep, 
+#      xaxt = "n", xlab = "Step", ylab = "RMSE", 
+#      pch = 21, bg = "firebrick", col = NA, 
+#      main = "Test Error ~ Step",
+#      cex = 1.5, cex.axis = 1.5, cex.main = 2, cex.lab = 1.5)
+# lines(seq(70,150,10), mseStep)
+# axis(side=1, at=seq(70,150,10), cex.axis = 1.5)
+# grid()
+# 
+# plot(1:90, mseDistep, 
+#      xaxt = "n", xlab = "Step", ylab = "RMSE", 
+#      pch = 21, bg = rep(wes_palette("Zissou1", 10, type = "continuous"), 9), col = NA, 
+#      main = "Test Error ~ Step:Dist",
+#      cex = 1.5, cex.axis = 1.5, cex.main = 2, cex.lab = 1.5)
+# 
+# for (i in seq(0,90,10)) {
+#   lines((i+1):(i+10), mseDistep[(i+1):(i+10)])
+# }
+# 
+# axis(side=1, at=seq(5,90,10), las = 1, srt = 60, labels = seq(70, 150, 10), cex.axis = 1.5)
+# abline(v=seq(0,90,10), lty="dashed", col = "gray")
+# 
+# legend.col(wes_palette("Zissou1", 10, type = "continuous"), seq(3,30,3))
+```
+
+###RMSE ggplot (dashed lines)
+
+```{r, fig.height=5, fig.width=20}
+# Store MSE in data frame
+mseFrame <- data.frame(mse = as.numeric(mseDistep), 
+                       x = rep(1:90, 4), 
+                       step = rep(rep(seq(70,150,10), each = 10), 4), 
+                       dist = rep(rep(seq(3,30,3), 9), 4),
+                       axType = rep(c("Type1", "Type2", "Type3", "Type4"), each = 90))
+
+# Color palette
+pal <- wes_palette("Zissou1", 10, type = "continuous")
+rmseList <- list()
+
+# Plot and store in list
+for (i in unique(mseFrame$axType)) {
+  
+  mseType <- mseFrame[mseFrame$axType == i, ]
+  
+  rmseList[[i]] <- ggplot(data=mseType, aes(x=x, y=mse, group=step)) +
+    geom_line() +
+    geom_point(aes(colour = dist)) + 
+    scale_x_continuous(breaks=seq(5,90,10), labels = seq(70,150,10)) +
+    scale_y_continuous(breaks=seq(2000, 16000, 2000), limits = c(2000,16000)) +
+    scale_colour_gradientn(name = "Distance",
+                           colours = pal,
+                           breaks = seq(3,30,3)) +
+    # geom_vline(xintercept=seq(0,90,10), colour = "gray", linetype="dashed", ) +
+    # geom_smooth() +
+    theme_bw() + 
+    theme(panel.grid.major = element_blank(),
+          panel.grid.minor = element_blank(),
+          legend.key.size = unit(1.9, "cm"),
+          legend.key.width = unit(0.5,"cm"),
+          plot.title = element_text(hjust = 0.5, size = 22),
+          axis.text.x = element_text(size = 13, color = "black"),
+          axis.text.y = element_text(size = 13, color = "black"),
+          axis.title.x = element_text(size = 16),
+          axis.title.y = element_text(size = 16),
+          panel.background = element_rect(colour = "black",size = 1)) +
+    xlab("Step") +
+    ylab("RMSE") +
+    ggtitle(i)
+  
+}
+
+# Save as png
+png(filename=paste("rmseEstereo2.png", sep=""), type="cairo",
+    units="in", width=22, height=5, pointsize=12,
+    res=300)
+
+gridExtra::grid.arrange(grobs = rmseList, ncol = 4, nrow = 1)
+
+dev.off()
+```
+
+###RMSE ggplot (smooth bands)
+
+```{r, fig.height=5, fig.width=20}
+# Store MSE in data frame
+mseFrame <- data.frame(mse = as.numeric(mseDistep), 
+                       x = rep(1:90, 4), 
+                       step = rep(rep(seq(70,150,10), each = 10), 4), 
+                       dist = rep(rep(seq(3,30,3), 9), 4),
+                       axType = rep(c("Type1", "Type2", "Type3", "Type4"), each = 90))
+
+# Color palette
+pal <- wes_palette("Zissou1", 10, type = "continuous")
+rmseList <- list()
+
+# Plot and store in list
+for (i in unique(mseFrame$axType)) {
+  
+  mseType <- mseFrame[mseFrame$axType == i, ]
+  
+  rmseList[[i]] <- ggplot(data=mseType, aes(x=x, y=mse, group=step)) +
+    # geom_line() +
+    geom_smooth(colour = "black", method = "loess", span = 1.2, size = 0.5) +
+    geom_point(aes(colour = dist)) + 
+    geom_point(shape = 1, colour = "black") +
+    scale_x_continuous(breaks=seq(5,90,10), labels = seq(70,150,10)) +
+    scale_y_continuous(breaks=seq(2000, 16000, 2000), limits = c(2000,16000)) +
+    scale_colour_gradientn(name = "Distance",
+                           colours = pal,
+                           breaks = seq(3,30,3)) +
+    # geom_vline(xintercept=seq(0,90,10), colour = "gray", linetype="dashed", ) +
+    theme_bw() + 
+    theme(panel.grid.major = element_blank(),
+          panel.grid.minor = element_blank(),
+          legend.key.size = unit(1.9, "cm"),
+          legend.key.width = unit(0.5,"cm"),
+          plot.title = element_text(hjust = 0.5, size = 22),
+          axis.text.x = element_text(size = 13, color = "black"),
+          axis.text.y = element_text(size = 13, color = "black"),
+          axis.title.x = element_text(size = 16),
+          axis.title.y = element_text(size = 16),
+          panel.background = element_rect(colour = "black",size = 1)) +
+    xlab("Step") +
+    ylab("RMSE") +
+    ggtitle(i)
+  
+}
+
+# Save as png
+png(filename=paste("rmseEstereo_loess_border.png", sep=""), type="cairo",
+    units="in", width=22, height=5, pointsize=12,
+    res=300)
+
+gridExtra::grid.arrange(grobs = rmseList, ncol = 4, nrow = 1)
+
+dev.off()
+```
+
+###RMSE smooth + examples
+
+```{r, fig.height=5, fig.width=20}
+# Color palette
+pal <- wes_palette("Zissou1", 10, type = "continuous")
+list370 <- list()
+list30150 <- list()
+
+# Plot and store in list
+for (i in unique(allData$neurType)) {
+
+  #EXAMPLE 3.70
+  data2d <- allData[allData$neurType == i & allData$dist == 3 & allData$step == 70, ]
+  dataReal <- data.frame(LR = data2d$LR, LRe = data2d$LR)
+  
+  list370[[i]] <- ggplot(data=data2d, aes(x=LRe, y=LR)) +
+    geom_point(colour = "black", size = 0.6) + 
+    # geom_point(shape = 1, colour = "black") +
+    geom_point(data = dataReal, colour = "red", size = 0.6) + 
+    # geom_point(data = dataReal, shape = 1, colour = "black") +
+    theme_bw() + 
+    theme(panel.grid.major = element_blank(),
+          panel.grid.minor = element_blank(),
+          plot.title = element_text(hjust = 0.5, size = 12),
+          axis.title.x=element_text(size = 10),
+          axis.text.x=element_blank(),
+          axis.ticks.x=element_blank(),
+          axis.title.y=element_text(size = 10),
+          axis.text.y=element_blank(),
+          axis.ticks.y=element_blank(), 
+          panel.background = element_rect(colour = "black",size = 1)) +
+    ggtitle("70.3") + 
+    xlab("Predicted length") +
+    ylab("Real Length")
+  
+  #EXAMPLE 30.150
+  data2d <- allData[allData$neurType == i & allData$dist == 30 & allData$step == 150, ]
+  dataReal <- data.frame(LR = data2d$LR, LRe = data2d$LR)
+  
+  list30150[[i]] <- ggplot(data=data2d, aes(x=LRe, y=LR)) +
+    geom_point(colour = "black",  size = 0.6) + 
+    # geom_point(shape = 1, colour = "black") +
+    geom_point(data = dataReal, colour = "red", size = 0.6) + 
+    # geom_point(data = dataReal, shape = 1, colour = "black") +
+    theme_bw() + 
+    theme(panel.grid.major = element_blank(),
+          panel.grid.minor = element_blank(),
+          plot.title = element_text(hjust = 0.5, size = 12),
+          axis.title.x= element_text(size = 10),
+          axis.text.x=element_blank(),
+          axis.ticks.x=element_blank(),
+          axis.title.y= element_text(size = 10),
+          axis.text.y=element_blank(),
+          axis.ticks.y=element_blank(),
+          panel.background = element_rect(colour = "black",size = 1)) + 
+    ggtitle("150.30") + 
+    xlab("Real Model") +
+    ylab("Real Length")
+  
+}
+
+# Save as png
+for (i in unique(allData$neurType)) {
+  
+  png(filename=paste("rmseEstereo_loess_border_examples_", i, ".png", sep=""), type="cairo",
+      units="in", width=2.5, height=1.5, pointsize=12,
+      res=300)
+  
+  gridExtra::grid.arrange(grobs = c(list370[i], list30150[i]), ncol = 2, nrow = 1)
+  
+  
+  dev.off()
+  
+}
+
+```
+
+##Mean error
+
+```{r, fig.height=10, fig.width=20}
+# Estimate mean residual error
+meanFrame <- aggregate(error ~ dist + step + neurType, allData, mean)
+
+# Library
+library(latticeExtra)
+
+# Contour color palette
+# colfunc <- colorRampPalette(c("navy","royalblue","springgreen","gold","yellow"))
+# coul <- colfunc(1000)
+# library(viridisLite)
+coul <- pals::parula(10000)
+
+
+# Store heatmaps
+heatList <- list()
+smoothList <- list()
+
+for (i in c("Type1", "Type2", "Type3", "Type4")) {
+  
+  heatList[[i]] <- levelplot(error ~ dist * step, data = meanFrame[meanFrame$neurType == i, ], col.regions = coul,
+                             scales = list(y = list(at = unique(allData$step), labels = unique(allData$step), cex = 1.5),
+                                           x = list(at = unique(allData$dist), labels = unique(allData$dist)), cex = 1.5),
+                             ylim = c(155, 65),
+                             ylab = list(label = "Step", cex = 1.5),
+                             xlab = list(label = "Distance", cex = 1.5),
+                             main = list(label = paste(i, "Mean % Error"), cex = 2),
+                             colorkey = list(at = (seq(min(meanFrame[meanFrame$neurType == i, "error"]),
+                                                       max(meanFrame[meanFrame$neurType == i, "error"]),
+                                                       0.05)),
+                                             cex = 1.2)) 
+  
+  smoothList[[i]] <- levelplot(error ~ dist * step, data = meanFrame[meanFrame$neurType == i, ], col.regions = coul,
+                               scales = list(y = list(at = unique(allData$step), labels = unique(allData$step), cex = 1.3),
+                                             x = list(at = unique(allData$dist), labels = unique(allData$dist)), cex = 1.3),
+                               ylim = c(150,70),
+                               xlim = c(3,30),
+                               ylab = list(label = "Step", cex = 1.5),
+                               xlab = list(label = "Distance", cex = 1.5),
+                               main = list(label = paste(i, "Smooth"), cex = 1.5),
+                               cuts = 20, 
+                               colorkey = list(at = (seq(min(meanFrame[meanFrame$neurType == i, "error"]),
+                                                         max(meanFrame[meanFrame$neurType == i, "error"]),
+                                                         0.05)),
+                                               cex = 1.2),
+                               panel = panel.levelplot.points, cex = 0) + 
+    layer_(panel.2dsmoother(..., n = 200))
+  
+}
+
+
+# Plot
+
+png(filename=paste("meanErrorEstereo.png", sep=""), type="cairo",
+    units="in", width=22, height=10, pointsize=12,
+    res=300)
+
+gridExtra::grid.arrange(grobs = c(heatList, smoothList), ncol = 4, nrow = 2)
+
+dev.off()
+
+```
+
+##Mean error (Type123 vs 4)
+
+```{r, fig.height=10, fig.width=20}
+# Estimate mean residual error
+meanFrame <- aggregate(error ~ dist + step + neurType, allData, mean)
+
+# Library
+library(latticeExtra)
+
+# Contour color palette
+# colfunc <- colorRampPalette(c("navy","royalblue","springgreen","gold","yellow"))
+# coul <- colfunc(1000)
+# library(viridisLite)
+coul <- pals::parula(10000)
+
+
+# Store heatmaps
+heatList <- list()
+smoothList <- list()
+
+# Limits
+min123 <- min(meanFrame[meanFrame$neurType != "Type4", "error"])
+max123 <- max(meanFrame[meanFrame$neurType != "Type4", "error"])
+
+for (i in c("Type1", "Type2", "Type3", "Type4")) {
+  
+  if (i != "Type4") {
+    
+  heatList[[i]] <- levelplot(error ~ dist * step, data = meanFrame[meanFrame$neurType == i, ], col.regions = coul,
+                             scales = list(y = list(at = unique(allData$step), labels = unique(allData$step), cex = 1.5),
+                                           x = list(at = unique(allData$dist), labels = unique(allData$dist)), cex = 1.5),
+                             ylim = c(155, 65),
+                             ylab = list(label = "Step", cex = 1.5),
+                             xlab = list(label = "Distance", cex = 1.5),
+                             main = list(label = paste(i, "Mean % Error"), cex = 2),
+                             colorkey = list(at = (seq(min123,
+                                                       max123,
+                                                       0.05)),
+                                             cex = 1.2)) 
+  
+  smoothList[[i]] <- levelplot(error ~ dist * step, data = meanFrame[meanFrame$neurType == i, ], col.regions = coul,
+                               scales = list(y = list(at = unique(allData$step), labels = unique(allData$step), cex = 1.3),
+                                             x = list(at = unique(allData$dist), labels = unique(allData$dist)), cex = 1.3),
+                               ylim = c(150,70),
+                               xlim = c(3,30),
+                               ylab = list(label = "Step", cex = 1.5),
+                               xlab = list(label = "Distance", cex = 1.5),
+                               main = list(label = paste(i, "Smooth"), cex = 1.5),
+                               cuts = 20, 
+                               colorkey = list(at = (seq(min123,
+                                                         max123,
+                                                         0.05)),
+                                               cex = 1.2),
+                               panel = panel.levelplot.points, cex = 0) + 
+    layer_(panel.2dsmoother(..., n = 200))
+  
+  } else {
+    
+     heatList[[i]] <- levelplot(error ~ dist * step, data = meanFrame[meanFrame$neurType == i, ], col.regions = coul,
+                             scales = list(y = list(at = unique(allData$step), labels = unique(allData$step), cex = 1.5),
+                                           x = list(at = unique(allData$dist), labels = unique(allData$dist)), cex = 1.5),
+                             ylim = c(155, 65),
+                             ylab = list(label = "Step", cex = 1.5),
+                             xlab = list(label = "Distance", cex = 1.5),
+                             main = list(label = paste(i, "Mean % Error"), cex = 2),
+                             colorkey = list(at = (seq(min(meanFrame[meanFrame$neurType == i, "error"]),
+                                                       max(meanFrame[meanFrame$neurType == i, "error"]),
+                                                       0.05)),
+                                             cex = 1.2)) 
+  
+  smoothList[[i]] <- levelplot(error ~ dist * step, data = meanFrame[meanFrame$neurType == i, ], col.regions = coul,
+                               scales = list(y = list(at = unique(allData$step), labels = unique(allData$step), cex = 1.3),
+                                             x = list(at = unique(allData$dist), labels = unique(allData$dist)), cex = 1.3),
+                               ylim = c(150,70),
+                               xlim = c(3,30),
+                               ylab = list(label = "Step", cex = 1.5),
+                               xlab = list(label = "Distance", cex = 1.5),
+                               main = list(label = paste(i, "Smooth"), cex = 1.5),
+                               cuts = 20, 
+                               colorkey = list(at = (seq(min(meanFrame[meanFrame$neurType == i, "error"]),
+                                                         max(meanFrame[meanFrame$neurType == i, "error"]),
+                                                         0.05)),
+                                               cex = 1.2),
+                               panel = panel.levelplot.points, cex = 0) + 
+    layer_(panel.2dsmoother(..., n = 200))
+    
+  }
+
+  
+}
+
+
+# Plot
+
+png(filename=paste("meanErrorEstereo_DIFY.png", sep=""), type="cairo",
+    units="in", width=22, height=10, pointsize=12,
+    res=300)
+
+gridExtra::grid.arrange(grobs = c(heatList, smoothList), ncol = 4, nrow = 2)
+
+dev.off()
+
+```
+
+
+##Errors vs cumul probability
+
+```{r, fig.height=5, fig.width=20}
+# Inverse quantile
+quantInv <- function(distr, value) ecdf(distr)(value)
+
+# Function to apply to all axon types
+quantType <- function(peFrame, probSeq) {
+  
+  probList <- list()
+  
+  for (i in unique(peFrame$interact)) {
+    
+    dataProb <- peFrame[peFrame$interact == i, "error"]
+    probVec <- numeric()
+    
+    for (j in probSeq) {
+      
+      errProb <- quantInv(dataProb, j) - quantInv(dataProb, -j)
+      probVec <- c(probVec, errProb)
+      
+    }
+    
+    probList[[i]] <- probVec
+    
+  }
+  
+  return(probList)
+  
+}
+
+# Define errors for which to calculate probability
+binProb <- 0.5
+probSeq <- seq(binProb, 100, binProb)
+# Store axon types in lists
+frameList <- list(Type1 = allData[allData$neurType == "Type1", ],
+                        Type2 = allData[allData$neurType == "Type2", ],
+                        Type3 = allData[allData$neurType == "Type3", ],
+                        Type4 = allData[allData$neurType == "Type4", ])
+axProb <- lapply(frameList, function(x) quantType(x, probSeq))
+saveRDS(axProb, "errorProbs_RMSE.rds")
+```
+
+######Plot heatmap 5, 10
+
+```{r, fig.height=10, fig.width=15}
+library(latticeExtra)
+
+# Reformat as dataframe
+binProb <- 0.5
+probSeq <- seq(binProb, 100, binProb)
+axProb <- readRDS("errorProbs_RMSE.rds")
+probFrame <- data.frame(error = rep(probSeq, 90*4), 
+                        neurType = rep(c("Type1", "Type2", "Type3", "Type4"), each = length(probSeq)*90),
+                        dist = rep(rep(unique(allData$dist), each = length(probSeq)*9), 4),
+                        step = rep(rep(unique(allData$step), each = length(probSeq)), 10*4),
+                        prob = unlist(axProb))
+
+# Plot contour plot for 5% error
+coul <- viridis::magma(10000)
+
+# Store heatmaps
+heatProbList <- list()
+smoothProbList <- list()
+
+for (i in c("Type1", "Type2", "Type3", "Type4")) {
+  
+  levelList1 <- list()
+  levelList2 <- list()
+  
+  for (j in c(5, 10, 20)) {
+    
+    dataPlot <- probFrame[probFrame$neurType == i & probFrame$error == j, ]
+    
+    levelList1[[as.character(j)]] <- levelplot(prob ~ dist * step, data = dataPlot, 
+                                               col.regions = coul,
+                                               scales = list(y = list(at = unique(allData$step), labels = unique(allData$step), cex = 1.3),
+                                                             x = list(at = unique(allData$dist), labels = unique(allData$dist)), cex = 1.3),
+                                               ylim = c(155,65),
+                                               colorkey = list(at = (seq(min(dataPlot$prob),
+                                                                         max(dataPlot$prob),
+                                                                         0.0025))),
+                                               ylab = list(label = "Step", cex = 1.5),
+                                               xlab = list(label = "Distance", cex = 1.5),
+                                               main = list(label = paste(i, " P(Error <= ", j, " %)", sep = ""), cex = 2)) 
+    
+    levelList2[[as.character(j)]] <- levelplot(prob ~ dist * step, data = dataPlot, 
+                                               col.regions = coul,
+                                               scales = list(y = list(at = unique(allData$step), labels = unique(allData$step), cex = 1.3),
+                                                             x = list(at = unique(allData$dist), labels = unique(allData$dist)), cex = 1.3),
+                                               ylim = c(150,70),
+                                               xlim = c(3,30),
+                                               cuts = 20,
+                                               colorkey = list(at = (seq(min(dataPlot$prob),
+                                                                         max(dataPlot$prob),
+                                                                         0.0025))),
+                                               ylab = list(label = "Step", cex = 1.5),
+                                               xlab = list(label = "Distance", cex = 1.5),
+                                               main = list(label = paste(i, " P(Error <= ", j, " %)", sep = ""), cex = 2),
+                                               panel = panel.levelplot.points, cex = 0) + 
+      layer_(panel.2dsmoother(..., n = 200))
+  }
+  
+  heatProbList[[i]] <- levelList1
+  smoothProbList[[i]] <- levelList2
+}
+
+smoothProb5 <- lapply(smoothProbList, function(x) x[[1]])
+smoothProb10 <- lapply(smoothProbList, function(x) x[[2]])
+smoothProb20 <- lapply(smoothProbList, function(x) x[[3]])
+
+# Plot
+# setwd("EstereoAnalysis/")
+for (i in c("Type1", "Type2", "Type3", "Type4")) {
+
+  png(filename=paste("errorProbabilities_heatmap_", i, "_MSE.png", sep=""), type="cairo",
+      units="in", width=15, height=10, pointsize=12,
+      res=300)
+
+  gridExtra::grid.arrange(grobs = c(heatProbList[[i]], smoothProbList[[i]]), ncol = 3, nrow = 2)
+  
+  dev.off()
+
+}
+
+
+```
+
+
+##Common plot
+
+```{r, fig.height=40, fig.width=40}
+gridExtra::grid.arrange(grobs = c(rmseList, smoothList, smoothProb5, smoothProb10), ncol = 4, nrow = 4)
+```
+
+##Boxplot error~dist:step
+
+```{r, fig.width=20, fig.height=5}
+# Sort by step
+dataStep <- allData[order(allData$step), ]
+dataStep$x <- rep(1:10, 8559)
+
+# Color palette
+pal <- wes_palette("Zissou1", 10, type = "continuous")
+errorList <- list()
+
+# ggplot(dataStep, aes(x=factor(step), y=error, fill=factor(dist))) +
+#   geom_boxplot() +
+#   scale_fill_manual(values=pal) + 
+#   ylim(c(-5,100))
+
+# Plot and store in list
+for (i in unique(allData$neurType)) {
+  
+  errType <- dataStep[dataStep$neurType == i, ]
+  
+  meanErr <- aggregate(error ~ step + dist, errType, mean)
+  meanErr <- meanErr[order(meanErr$step), ]
+  
+  errorList[[i]] <- ggplot(errType, aes(x=factor(step), y=error)) +
+    geom_boxplot(aes(fill=interact), outlier.shape = NA) +
+    scale_fill_manual(values=pal) + 
+    # geom_line() +
+    # geom_boxplot() + 
+    # scale_x_continuous(breaks=seq(5,90,10), labels = seq(70,150,10)) +
+    # scale_y_continuous(breaks=seq(2000, 16000, 2000), limits = c(2000,16000)) +
+    # scale_colour_gradientn(name = "Distance",
+    #                        colours = pal,
+    #                        breaks = seq(3,30,3)) +
+    # geom_vline(xintercept=seq(0,90,10), colour = "gray", linetype="dashed", ) +
+    geom_smooth(method="loess", se=TRUE) +
+    theme_bw() + 
+    theme(panel.grid.major = element_blank(),
+          panel.grid.minor = element_blank(),
+          legend.key.size = unit(1, "cm"),
+          legend.key.width = unit(0.5,"cm"),
+          plot.title = element_text(hjust = 0.5, size = 22),
+          axis.text.x = element_text(size = 13, color = "black"),
+          axis.text.y = element_text(size = 13, color = "black"),
+          axis.title.x = element_text(size = 16),
+          axis.title.y = element_text(size = 16)) +
+    xlab("Step") +
+    ylab("RMSE") +
+    ggtitle(i) +
+    scale_y_continuous(limits = c(0,200))
+  
+}
+
+# Save as png
+# png(filename=paste("rmseEstereo2.png", sep=""), type="cairo",
+#     units="in", width=22, height=5, pointsize=12,
+#     res=300)
+
+gridExtra::grid.arrange(grobs = errorList, ncol = 4, nrow = 1)
+
+# dev.off()
+```
+
+
+```{r}
+data=data.frame(date=as.Date(c("2011-02-10","2011-02-10","2011-02-10","2011-02-10","2011-02-10",
+                         "2011-02-10","2011-02-10","2011-02-10","2011-02-10","2011-02-10",
+                         "2011-02-20","2011-02-20","2011-02-20","2011-02-20","2011-02-20",
+                         "2011-02-20","2011-02-20","2011-02-20","2011-02-20","2011-02-20",
+                         "2011-02-28","2011-02-28","2011-02-28","2011-02-28","2011-02-28",
+                         "2011-02-28","2011-02-28","2011-02-28","2011-02-28","2011-02-28",
+                         "2011-03-10","2011-03-10","2011-03-10","2011-03-10","2011-03-10",
+                         "2011-03-10","2011-03-10","2011-03-10","2011-03-10","2011-03-10"),format="%Y-%m-%d"),
+            spore=c(0,1,0,1,0,
+                    1,2,0,1,1,
+                    8,5,6,12,7,
+                    7,5,4,7,6,
+                    18,24,25,32,14,
+                    24,16,18,23,30,
+                    27,26,36,31,22,
+                    28,29,37,30,24),
+            house = rep(c("Int", "Ext"), each = 5, times = 4)
+            )
+data$interact <- interaction(data$date, data$house)
+
+
+ggplot(data,aes(x=date,y=spore)) + 
+  geom_boxplot(aes(fill=interact), outlier.shape = NA) + 
+  # scale_fill_manual(values=rep(pal,9)) +
+  geom_smooth() + 
+  # scale_y_continuous(limits = c(0,200)) +
+  theme(legend.position = "none")
+
+
+ggplot(allData,aes(x=step,y=error)) + 
+  geom_boxplot(aes(fill=interact), outlier.shape = NA) + 
+  scale_fill_manual(values=rep(pal,9)) +
+  geom_smooth(method="loess") + 
+  scale_y_continuous(limits = c(0,200)) +
+  theme(legend.position = "none")
+```
+
+

Code/Statistical Analysis/modelFree_stereology_errors_CV.Rmd

@@ -0,0 +1,788 @@
+---
+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(ggpubr)
+library(ggridges)
+library(mltools)
+library(data.table)
+library(caret)
+library(interactions)
+
+options(digits = 10)
+
+# https://rstudio-pubs-static.s3.amazonaws.com/297778_5fce298898d64c81a4127cf811a9d486.html
+
+
+#############################
+# Bootstrap coefficients LM #
+#############################
+
+# For models with interactions
+
+bootCoefintLM <- function(mydata, modForm, nreps) {
+  
+  # Fit model
+  modLM <- lm(formula(modForm), data = mydata)
+
+  # Ensure reproducibility
+  # set.seed(12345) 
+  
+  # Set up a matrix to store the results (1 intercept + 1 predictor = 2)
+  coefmat <- matrix(NA, nreps, length(coef(modLM)))
+  
+  # We need the fitted values and the residuals from the model
+  resids <- residuals(modLM)
+  preds <- fitted(modLM)
+  
+  # Repeatedly generate bootstrapped responses
+  for(i in 1:nreps) {
+    booty <- preds + sample(resids, rep=TRUE)
+    modLM2 <- update(modLM, booty ~ .)
+    coefmat[i,] <- coef(modLM2)
+    print(i)
+  }
+  
+  # Create a dataframe to store predictors values
+  colnames(coefmat) <- names(coef(modLM2))
+  coefmat <- data.frame(coefmat)
+  
+  return(coefmat)
+}
+
+
+################
+# 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, do.coef = FALSE)$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 # STEREOLOGY VERSION
+##############################
+
+the.replication.type <- function(reg, s, x_Np1 = 0, x_Np2 = 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
+  dist <- model.frame(reg)[,2]
+  step <- model.frame(reg)[,3]
+  bs.reg <- lm(y.star ~ dist:step - 1)
+
+  # Create bootstrapped adjusted residuals
+  bs.lev <- influence(bs.reg, do.coef = FALSE)$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*x_Np2
+  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("EstereoDataPlanos/", sep=""), full.names=T)
+
+# Load matlab file
+axonFiles <- lapply(axonNames, function(x) readMat(x))
+names(axonFiles) <- c("Type1", "Type2", "Type3", "Type4")
+
+# Extract data
+# errorMatrix contains 4 arrays, one per neuron type
+# within each array, the dimensions corresponds to step:dist:neuron
+# this means that each array has as many matrices as neuron of a certain type
+dist         <- lapply(axonFiles, function(x) x$Distancias.Entre.Planos)[[1]]
+step         <- lapply(axonFiles, function(x) x$Step.Lengths)[[1]]
+errorMatrix  <- lapply(axonFiles, function(x) x$MATRIZ.ERROR.LENGTH)
+lpMatrix     <- lapply(axonFiles, function(x) x$MATRIZ.ESTIMATED.AXON.LENGTH)
+lrVals       <- lapply(axonFiles, function(x) x$AXON.REAL.LENGTH)
+
+# Get number of neurons per type
+numberTypes  <- unlist(lapply(errorMatrix, function(x) dim(x)[3]))
+
+# Vectorizing the arrays goes as follows: neuron, dist, step
+errorVector <- unlist(lapply(errorMatrix, function(x) as.vector(x)))
+lpVector <- unlist(lapply(lpMatrix, function(x) as.vector(x)))
+lrVector <- unlist(lapply(lrVals, function(x) as.vector(x)))
+
+# Store in data frames
+allData <- data.frame(id = rep(1:sum(numberTypes), each = 90),
+                      neurType = rep(c("Type1", "Type2", "Type3", "Type4"), times = 90*numberTypes),
+                      dist = rep(rep(dist, each = 9), sum(numberTypes)), 
+                      step = rep(rep(step, 10), sum(numberTypes)),
+                      error = abs(errorVector),
+                      LP = lpVector,
+                      LR = rep(lrVector, each = 90))
+```
+
+
+##Bootstrap PI - leave dist out
+
+```{r}
+# Read data
+frameList <- list()
+
+for (h in c("Type1", "Type2", "Type3", "Type4")) {
+  
+  bootout <- readRDS(paste("EstereoAnalysis/stereo_bootPI_leaveOut_", h, ".rds", sep=""))
+  y.p <- bootout[[1]]
+  epList <- bootout[[2]]
+  
+  # Add errors to predictions
+  epList.sum <- list()
+  
+  for (i in 1:90) {
+    epList.sum[[i]] <- epList[[i]] + y.p[[i]]
+  }
+  
+  names(epList.sum) <- paste(rep(unique(allData$dist), each = 9), rep(unique(allData$step), 10), sep=".")
+  
+  # Reformat as dataframe
+  peFrame <- data.frame(neurType = rep(h, 90*100),
+                        dist = rep(unique(allData$dist), each = 9*100),
+                        step = rep(rep(unique(allData$step), each = 100), 10),
+                        pe = unlist(epList.sum))
+  peFrame$interact <- interaction(peFrame$dist, peFrame$step, lex.order = TRUE)
+  peFrame$abspe <- abs(peFrame$pe)
+  
+  # Store in list
+  frameList[[h]] <- peFrame
+  
+}
+
+# Join into a common data frame
+peFrame <- do.call(rbind, frameList)
+```
+
+
+###Mean error
+
+```{r, fig.height=10, fig.width=20}
+meanFrame <- aggregate(abspe ~ dist + step + neurType, peFrame, mean)
+
+# Library
+library(latticeExtra)
+
+# Contour color palette
+colfunc <- colorRampPalette(c("navy","royalblue","springgreen","gold","yellow"))
+colfin <- colfunc(1000)
+library(viridisLite)
+coul <- viridis(1000)
+
+
+# Store heatmaps
+heatList <- list()
+smoothList <- list()
+
+for (i in c("Type1", "Type2", "Type3", "Type4")) {
+  
+  heatList[[i]] <- levelplot(abspe ~ dist * step, data = meanFrame[meanFrame$neurType == i, ], col.regions = coul,
+                             scales = list(y = list(at = unique(allData$step), labels = unique(allData$step)),
+                                           x = list(at = unique(allData$dist), labels = unique(allData$dist))),
+                             ylim = c(155, 65),
+                             colorkey = list(at = (seq(min(meanFrame[meanFrame$neurType == i, "abspe"]),
+                                                       max(meanFrame[meanFrame$neurType == i, "abspe"]),
+                                                       0.05))),
+                             main = paste("Mean Abs Error,", i)) 
+  
+  smoothList[[i]] <- levelplot(abspe ~ dist * step, data = meanFrame[meanFrame$neurType == i, ], col.regions = coul,
+                               scales = list(y = list(at = unique(allData$step), labels = unique(allData$step)),
+                                             x = list(at = unique(allData$dist), labels = unique(allData$dist))),
+                               ylim = c(150,70),
+                               xlim = c(3,30),
+                               cuts = 20, 
+                               colorkey = list(at = (seq(min(meanFrame[meanFrame$neurType == i, "abspe"]),
+                                                         max(meanFrame[meanFrame$neurType == i, "abspe"]),
+                                                         0.05))),
+                               main = paste("Mean Abs Error,", i),
+                               panel = panel.levelplot.points, cex = 0) + 
+    layer_(panel.2dsmoother(..., n = 200))
+  
+}
+
+
+# Plot
+setwd("EstereoAnalysis/")
+
+png(filename="meanProbabilities_heatmap_CVleaveout.png", type="cairo",
+    units="in", width=20, height=10, pointsize=12,
+    res=300)
+
+  gridExtra::grid.arrange(grobs = c(heatList, smoothList), ncol = 4, nrow = 2)
+
+dev.off()
+```
+
+###Plot density
+
+```{r, fig.height=10, fig.width=20}
+# Plot
+# peGroup <- peFrame %>% group_by(interact)
+# peSample <- sample_n(peGroup, 100)
+
+plotList <- list()
+
+for (m in c("Type1", "Type2", "Type3", "Type4")) {
+  
+  plotList[[m]] <- ggplot(peFrame[peFrame$neurType == m, ], 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)) +
+    ggtitle(m)
+  
+}
+
+setwd("EstereoAnalysis/")
+png(filename=paste("prederrorDistributions_ols_CVleaveout.png", sep=""), type="cairo",
+    units="in", width=20, height=10, pointsize=12,
+    res=300)
+
+ggarrange(plotlist = plotList, ncol = 4, nrow = 1)
+
+dev.off()
+```
+
+###Errors vs cumul probability
+
+```{r, fig.height=5, fig.width=20}
+# Inverse quantile
+quantInv <- function(distr, value) ecdf(distr)(value)
+
+# Function to apply to all axon types
+quantType <- function(peFrame, probSeq) {
+  
+  probList <- list()
+  
+  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
+    
+  }
+  
+  return(probList)
+  
+}
+
+# Define errors for which to calculate probability
+binProb <- 0.5
+probSeq <- seq(binProb, 100, binProb)
+# Store axon types in lists
+axProb <- lapply(frameList, function(x) quantType(x, probSeq))
+# Save 
+saveRDS(axProb, "errorProbs_CVleaveout.rds")
+```
+
+######Plot heatmap
+
+```{r, fig.height=5, fig.width=15}
+library(latticeExtra)
+
+# Reformat as dataframe
+binProb <- 0.5
+probSeq <- seq(binProb, 100, binProb)
+axProb <- readRDS("errorProbs_CVleaveout.rds")
+probFrame <- data.frame(error = rep(probSeq, 90*4), 
+                        neurType = rep(c("Type1", "Type2", "Type3", "Type4"), each = length(probSeq)*90),
+                        dist = rep(rep(unique(allData$dist), each = length(probSeq)*9), 4),
+                        step = rep(rep(unique(allData$step), each = length(probSeq)), 10*4),
+                        prob = unlist(axProb))
+
+# Plot conttour plot for 5% error
+library(viridisLite)
+coul <- viridis(1000)
+
+# Store heatmaps
+heatList <- list()
+smoothList <- list()
+
+for (i in c("Type1", "Type2", "Type3", "Type4")) {
+  
+  levelList1 <- list()
+  levelList2 <- list()
+  
+  for (j in c(5, 10, 20)) {
+    
+    dataPlot <- probFrame[probFrame$neurType == i & probFrame$error == j, ]
+    
+    levelList1[[as.character(j)]] <- levelplot(prob ~ dist * step, data = dataPlot, 
+                                               col.regions = coul,
+                                               scales = list(y = list(at = unique(allData$step), labels = unique(allData$step)),
+                                                             x = list(at = unique(allData$dist), labels = unique(allData$dist))),
+                                               ylim = c(155,65),
+                                               colorkey = list(at = (seq(min(dataPlot$prob),
+                                                                         max(dataPlot$prob),
+                                                                         0.001))),
+                                               main = paste("P(%Error <=", j, ")", sep = "")) 
+    
+    levelList2[[as.character(j)]] <- levelplot(prob ~ dist * step, data = dataPlot, 
+                                               col.regions = coul,
+                                               scales = list(y = list(at = unique(allData$step), labels = unique(allData$step)),
+                                                             x = list(at = unique(allData$dist), labels = unique(allData$dist))),
+                                               ylim = c(150,70),
+                                               xlim = c(3,30),
+                                               cuts = 20,
+                                               colorkey = list(at = (seq(min(dataPlot$prob),
+                                                                         max(dataPlot$prob),
+                                                                         0.001))),
+                                               main = paste("P(%Error <=", j, ")", sep = ""),
+                                               panel = panel.levelplot.points, cex = 0) + 
+      layer_(panel.2dsmoother(..., n = 200))
+  }
+  
+  heatList[[i]] <- levelList1
+  smoothList[[i]] <- levelList2
+}
+
+# Plot
+setwd("EstereoAnalysis/")
+for (i in c("Type1", "Type2", "Type3", "Type4")) {
+
+  png(filename=paste("errorProbabilities_heatmap_", i, "_CVleaveout.png", sep=""), type="cairo",
+      units="in", width=15, height=10, pointsize=12,
+      res=300)
+
+  gridExtra::grid.arrange(grobs = c(heatList[[i]], smoothList[[i]]), ncol = 3, nrow = 2)
+  
+  dev.off()
+
+}
+
+
+```
+
+
+##Bootstrap PI - leave dist and step out
+
+```{r}
+# Read data
+frameList <- list()
+
+for (h in c("Type1", "Type2", "Type3", "Type4")) {
+  
+  boot2out <- readRDS(paste("EstereoAnalysis/stereo_bootPI_leave2Out_", h, ".rds", sep=""))
+  y.p <- boot2out[[1]]
+  epList <- boot2out[[2]]
+  
+  # Add errors to predictions
+  epList.sum <- list()
+  
+  for (i in 1:90) {
+    epList.sum[[i]] <- epList[[i]] + y.p[[i]]
+  }
+  
+  names(epList.sum) <- paste(rep(unique(allData$dist), each = 9), rep(unique(allData$step), 10), sep=".")
+  
+  # Reformat as dataframe
+  peFrame <- data.frame(neurType = rep(h, 90*100),
+                        dist = rep(unique(allData$dist), each = 9*100),
+                        step = rep(rep(unique(allData$step), each = 100), 10),
+                        pe = unlist(epList.sum))
+  peFrame$interact <- interaction(peFrame$dist, peFrame$step, lex.order = TRUE)
+  peFrame$abspe <- abs(peFrame$pe)
+  
+  # Store in list
+  frameList[[h]] <- peFrame
+  
+}
+
+# Join into a common data frame
+peFrame <- do.call(rbind, frameList)
+```
+
+###Mean error
+
+```{r, fig.height=10, fig.width=20}
+meanFrame <- aggregate(abspe ~ dist + step + neurType, peFrame, mean)
+
+# Library
+library(latticeExtra)
+
+# Contour color palette
+colfunc <- colorRampPalette(c("navy","royalblue","springgreen","gold","yellow"))
+colfin <- colfunc(1000)
+library(viridisLite)
+coul <- viridis(1000)
+
+
+# Store heatmaps
+heatList <- list()
+smoothList <- list()
+
+for (i in c("Type1", "Type2", "Type3", "Type4")) {
+  
+  heatList[[i]] <- levelplot(abspe ~ dist * step, data = meanFrame[meanFrame$neurType == i, ], col.regions = coul,
+                             scales = list(y = list(at = unique(allData$step), labels = unique(allData$step)),
+                                           x = list(at = unique(allData$dist), labels = unique(allData$dist))),
+                             ylim = c(155, 65),
+                             colorkey = list(at = (seq(min(meanFrame[meanFrame$neurType == i, "abspe"]),
+                                                       max(meanFrame[meanFrame$neurType == i, "abspe"]),
+                                                       0.05))),
+                             main = paste("Mean Abs Error,", i)) 
+  
+  smoothList[[i]] <- levelplot(abspe ~ dist * step, data = meanFrame[meanFrame$neurType == i, ], col.regions = coul,
+                               scales = list(y = list(at = unique(allData$step), labels = unique(allData$step)),
+                                             x = list(at = unique(allData$dist), labels = unique(allData$dist))),
+                               ylim = c(150,70),
+                               xlim = c(3,30),
+                               cuts = 20, 
+                               colorkey = list(at = (seq(min(meanFrame[meanFrame$neurType == i, "abspe"]),
+                                                         max(meanFrame[meanFrame$neurType == i, "abspe"]),
+                                                         0.05))),
+                               main = paste("Mean Abs Error,", i),
+                               panel = panel.levelplot.points, cex = 0) + 
+    layer_(panel.2dsmoother(..., n = 200))
+  
+}
+
+
+# Plot
+setwd("EstereoAnalysis/")
+
+png(filename="meanProbabilities_heatmap_CVleave2out.png", type="cairo",
+    units="in", width=20, height=10, pointsize=12,
+    res=300)
+
+  gridExtra::grid.arrange(grobs = c(heatList, smoothList), ncol = 4, nrow = 2)
+
+dev.off()
+```
+
+###Plot density
+
+```{r, fig.height=10, fig.width=20}
+# Plot
+# peGroup <- peFrame %>% group_by(interact)
+# peSample <- sample_n(peGroup, 100)
+
+plotList <- list()
+
+for (m in c("Type1", "Type2", "Type3", "Type4")) {
+  
+  plotList[[m]] <- ggplot(peFrame[peFrame$neurType == m, ], 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)) +
+    ggtitle(m)
+  
+}
+
+setwd("EstereoAnalysis/")
+png(filename=paste("prederrorDistributions_ols_CVleave2out.png", sep=""), type="cairo",
+    units="in", width=20, height=10, pointsize=12,
+    res=300)
+
+ggarrange(plotlist = plotList, ncol = 4, nrow = 1)
+
+dev.off()
+```
+
+
+##Bootstrap PI - CV
+
+
+```{r}
+# Read data
+frameList <- list()
+ycvList <- list()
+cvTotal <- list()
+
+for (g in c("Type1", "Type2", "Type3", "Type4")) {
+  
+  bootCV <- readRDS(paste("EstereoAnalysis/stereo_boot10CVPI_", g, ".rds", sep=""))
+  y.cv <- bootCV[[1]]
+  cvList <- bootCV[[2]]
+  
+  epCV <- list()
+  
+  for (h in 1:10) {
+    
+    # Add errors to predictions
+    epList.sum <- list()
+    
+    for (i in 1:90) {
+      
+      epList.sum[[i]] <- cvList[[h]][[i]] + y.cv[[h]][[i]]
+      
+    }
+    
+    names(epList.sum) <- paste(rep(unique(allData$dist), each = 9), rep(unique(allData$step), 10), sep=".")
+    epCV[[h]] <- epList.sum
+    
+  }
+  
+  # Reformat as dataframe
+  peFrame <- data.frame(neurType = rep(g, 90*100*10),
+                        dist = rep(rep(unique(allData$dist), each = 9*100), 10),
+                        step = rep(rep(rep(unique(allData$step), each = 100), 10), 10),
+                        pe = unlist(unlist(epCV)))
+  peFrame$interact <- interaction(peFrame$dist, peFrame$step, lex.order = TRUE)
+  peFrame$abspe <- abs(peFrame$pe)
+  
+  # Store in list
+  frameList[[g]] <- peFrame
+  ycvList[[g]] <- y.cv
+  cvTotal[[g]] <- cvList
+  
+}
+
+# Join into a common data frame
+peFrame <- do.call(rbind, frameList)
+```
+
+###Mean error
+
+```{r, fig.height=10, fig.width=20}
+meanFrame <- aggregate(abspe ~ dist + step + neurType, peFrame, mean)
+
+# Library
+library(latticeExtra)
+
+# Contour color palette
+colfunc <- colorRampPalette(c("navy","royalblue","springgreen","gold","yellow"))
+colfin <- colfunc(1000)
+library(viridisLite)
+coul <- viridis(1000)
+
+
+# Store heatmaps
+heatList <- list()
+smoothList <- list()
+
+for (i in c("Type1", "Type2", "Type3", "Type4")) {
+  
+  heatList[[i]] <- levelplot(abspe ~ dist * step, data = meanFrame[meanFrame$neurType == i, ], col.regions = coul,
+                             scales = list(y = list(at = unique(allData$step), labels = unique(allData$step)),
+                                           x = list(at = unique(allData$dist), labels = unique(allData$dist))),
+                             ylim = c(155, 65),
+                             colorkey = list(at = (seq(min(meanFrame[meanFrame$neurType == i, "abspe"]),
+                                                       max(meanFrame[meanFrame$neurType == i, "abspe"]),
+                                                       0.05))),
+                             main = paste("Mean Abs Error,", i)) 
+  
+  smoothList[[i]] <- levelplot(abspe ~ dist * step, data = meanFrame[meanFrame$neurType == i, ], col.regions = coul,
+                               scales = list(y = list(at = unique(allData$step), labels = unique(allData$step)),
+                                             x = list(at = unique(allData$dist), labels = unique(allData$dist))),
+                               ylim = c(150,70),
+                               xlim = c(3,30),
+                               cuts = 20, 
+                               colorkey = list(at = (seq(min(meanFrame[meanFrame$neurType == i, "abspe"]),
+                                                         max(meanFrame[meanFrame$neurType == i, "abspe"]),
+                                                         0.05))),
+                               main = paste("Mean Abs Error,", i),
+                               panel = panel.levelplot.points, cex = 0) + 
+    layer_(panel.2dsmoother(..., n = 200))
+  
+}
+
+
+# Plot
+setwd("EstereoAnalysis/")
+
+png(filename="meanProbabilities_heatmap_CV.png", type="cairo",
+    units="in", width=20, height=10, pointsize=12,
+    res=300)
+
+  gridExtra::grid.arrange(grobs = c(heatList, smoothList), ncol = 4, nrow = 2)
+
+dev.off()
+```
+
+###Plot density
+
+```{r, fig.height=10, fig.width=20}
+# Plot
+# peGroup <- peFrame %>% group_by(interact)
+# peSample <- sample_n(peGroup, 100)
+
+plotList <- list()
+
+for (m in c("Type1", "Type2", "Type3", "Type4")) {
+  
+  plotList[[m]] <- ggplot(peFrame[peFrame$neurType == m, ], 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)) +
+    ggtitle(m)
+  
+}
+
+setwd("EstereoAnalysis/")
+png(filename=paste("prederrorDistributions_ols_CV.png", sep=""), type="cairo",
+    units="in", width=20, height=10, pointsize=12,
+    res=300)
+
+ggarrange(plotlist = plotList, ncol = 4, nrow = 1)
+
+dev.off()
+```