Rubio_etal_2021/Code/Statistical Analysis/modelFree_stereology_CV_LRLP.Rmd

---
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)

options(digits = 10)

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

##############################
# Bootstrap PI per axon type # STEREOLOGY VERSION, LR vs LP
##############################

the.replication.type2 <- function(reg, s, disti, stepj, trainset, 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, do.coef = FALSE)$hat
  bs.s   <- residuals(bs.reg)
  
  bs.levAdd <- bs.lev[which(trainset$dist == disti & trainset$step == stepj)]
  bs.sAdd <- bs.s[which(trainset$dist == disti & trainset$step == stepj)]
  
  bs.Stud <- bs.sAdd/sqrt(1-bs.levAdd)
  bs.add  <- bs.Stud - mean(bs.Stud)

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

##############################
# Bootstrap PI per axon type # STEREOLOGY VERSION, LR vs LP
##############################

the.replication.type3 <- 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)

  # Estimation
  return(unname(coef(bs.reg)*x_Np1))
  
}

##############################
# 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),
                      error2 = errorVector,
                      LP = lpVector,
                      LR = rep(lrVector, each = 90))
allData$errorRaw <- allData$LR - allData$LP
```


##Bootstrap PI - Prediction CV - LR vs LP

This code is meant to be run independently for every neurType and in parallel

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

# Subset axon type
dataCV <- allData[allData$neurType == "Type1", ]

# List to store CV repetitions
cvList <- list ()
cvRand <- list()
y.cv <- list()

for (h in 1:10) {

  # Stratified random training/test sets
  stratDat <- stratified(dataCV, c("dist", "step"), 0.7, bothSets = TRUE, replace = FALSE)
  trainSet <- stratDat$SET1
  testSet <- stratDat$SET2
 
  # Store train and test set
  cvRand[[h]] <- list(trainSet = trainSet, testSet = 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)
  my.s.resid <- residuals(my.reg)
  
  # List to store errors
  epList <- list()
  # Vector to store predicted LR
  y.p <- numeric()
  
  
  # Bootstrapping on test set
  for (i in unique(allData$dist)) {
    
    for (j in unique(allData$step)) {
      
      # Subset data
      dataSub <- testSet[testSet$dist == i & testSet$step == j, ]
      
      kList <- list()
      kCount <- 1
      
      # Predict on unseen LP
      for (k in dataSub[ , "LP"]) {

        # Do bootstrap with 100 replications
        kList[[kCount]] <- replicate(n=100, the.replication.type3(reg = my.reg, s = my.s.resid, x_Np1 = k))
        y.p <- c(y.p, kList[[kCount]])
        kCount <- kCount + 1

        
      }
      
       # Store in list
      epList[[paste(as.character(i), as.character(j), sep = ".")]] <- kList
      print(paste("Step", j, "finished"))
      
    }
    
    print(paste("Dist", i, "finished"))
    
  }
  
  cvList[[h]] <- epList
  y.cv[[h]] <- y.p
  print(paste("CV", h, "finished"))
  
}


saveRDS(list(y.cv, cvList, cvRand), paste("stereo_boot10CVPI_LRLP_pred_Type1.rds", sep = ""))
```


##Bootstrap PI - CV - LR vs LP

This code is meant to be run independently for every neurType and in parallel

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

# Subset axon type
dataCV <- allData[allData$neurType == "Type1", ]

# List to store CV repetitions
cvList <- list ()
cvRand <- list()
y.cv <- list()

for (h in 1) {

  # Stratified random training/test sets
  stratDat <- stratified(dataCV, c("dist", "step"), 0.7, bothSets = TRUE, replace = FALSE)
  trainSet <- stratDat$SET1
  testSet <- stratDat$SET2
 
  # Store train and test set
  cvRand[[h]] <- list(trainSet = trainSet, testSet = 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)
  my.s.resid <- residuals(my.reg)
  
  # List to store errors
  epList <- list()
  # Vector to store predicted LR
  y.p <- numeric()
  
  
  # Bootstrapping on test set
  for (i in unique(allData$dist)) {
    
    for (j in unique(allData$step)) {
      
      # Subset data
      dataSub <- testSet[testSet$dist == i & testSet$step == j, ]
      
      kList <- list()
      kCount <- 1
      
      # Predict on unseen LP
      for (k in dataSub[ , "LP"]) {
        
        # Store predicted error
        y.p <- c(y.p, coef(my.reg)*k)
        
        # Create adjusted residuals
        leverage <- influence(my.reg, do.coef = FALSE)$hat
        my.s.resid <- residuals(my.reg)/sqrt(1-leverage)
        my.s.resid <- my.s.resid - mean(my.s.resid)
        
        # Do bootstrap with 100 replications
        kList[[kCount]] <- replicate(n=100, the.replication.type2(reg = my.reg, s = my.s.resid, disti = i, stepj = j, trainset = trainSet, 
                                                                  x_Np1 = k))
        kCount <- kCount + 1
        
      }
      
       # Store in list
      epList[[paste(as.character(i), as.character(j), sep = ".")]] <- kList
      print(paste("Step", j, "finished"))
      
    }
    
    print(paste("Dist", i, "finished"))
    
  }
  
  cvList[[h]] <- epList
  names(y.p) <- paste(rep(unique(allData$dist), each = 9), rep(unique(allData$step), 10), sep=".")
  y.cv[[h]] <- y.p
  print(paste("CV", h, "finished"))
  
}


saveRDS(list(y.cv, cvList, cvRand), paste("stereo_boot10CVPI_LRLP_notStud_Type1.rds", sep = ""))
```

##Plotting for bsAdd

```{r}
bsAdd <- readRDS("stereo_boot10CVPI_LRLP_bsAdd_Type1.rds")
y.cv <- bsAdd[[1]]
cvList <- bsAdd[[2]]
cvRand <- bsAdd[[3]]

estimList <- list()

for (h in 1) {
  lrSet <- cvRand[[h]][["testSet"]]
  lrOrder <- lrSet[with(lrSet, order(dist, step)), ]
  
  errorSet <- unlist(cvList[[h]])
  
  errorFrame <- data.frame(dist = rep(lrOrder$dist, each = 100), 
                           step = rep(lrOrder$step, each = 100),
                           LR = rep(lrOrder$LR, each = 100),
                           error = errorSet)
  
  estimList[[h]] <- errorFrame
}

estimTotal <- do.call(rbind, estimList)
estimTotal$pe <- (estimTotal$error/estimTotal$LR)*100
estimTotal$abspe <- abs(estimTotal$pe)
estimTotal$interact <- interaction(estimTotal$dist, estimTotal$step, lex.order = TRUE)

estimMean <- aggregate(abspe ~ dist + step, estimTotal, mean)
```

###Mean

```{r, fig.height=10, fig.width=20}
meanFrame <- estimMean
meanFrame$neurType <- rep("Type1", dim(meanFrame)[1])

# 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")) {
  
  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()
```

###Density

```{r,  fig.width=5, fig.height=20}
ggplot(estimTotal, 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))
  
```

###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))
axProb <- quantType(estimTotal, probSeq)
# Save 
# saveRDS(axProb, "errorProbs_CVleaveout.rds")
```

######Plot heatmap

```{r, width=15, height=10}
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))

probFrame <- data.frame(error = rep(probSeq, 90), 
                        neurType = rep("Type1", length(probSeq)*90),
                        dist = rep(unique(allData$dist), each = length(probSeq)*9),
                        step = rep(rep(unique(allData$step), each = length(probSeq)), 10),
                        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")) {
  
  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")) {

  png(filename=paste("prueba_", 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()

}


```


###PI 1

```{r}
jj <- allData[allData$neurType == "Type1", ]
boxplot(LP ~ dist + step, jj)
boxplot(LP ~ dist + step, cvRand[[1]]$trainSet)
boxplot(LP ~ dist + step, cvRand[[1]]$testSet)

boxplot(LR ~ dist + step, jj)
boxplot(LR ~ dist + step, cvRand[[1]]$trainSet)
boxplot(LR ~ dist + step, cvRand[[1]]$testSet)
```


```{r, fig.width=45, fig.height=50}
piLimits <- list()
cv1 <- cvList[[1]]

for (i in 1:90) {
  piLimits[[i]] <- lapply(cv1[[i]], function(x) quantile(x, c(0.025, 0.975)))
}

names(piLimits) <- names(cv1)
limVec <- unlist(unlist(piLimits))
oddPos <- seq(1, length(limVec), 2)

lowLims <- limVec[oddPos]
upLims <- limVec[-oddPos]
yp <- unlist(y.cv)

testDat <- cvRand[[1]]$testSet
testOrd <- testDat[order(testDat$dist, testDat$step), ]
testOrd$interact <- interaction(testOrd$dist, testOrd$step, lex.order = T)
lp <-testOrd$LP
names(lp) <- testOrd$interact
lr <- testOrd$LR
names(lr) <- testOrd$interact

  
# ySort <- sort(yp)
# lpSort <- lp[order(yp)]
# lrSort <- lr[order(yp)]
# lowSort <- lowLims[order(yp)]
# upSort <- upLims[order(yp)]

lpSort <- sort(lp)
ySort <- yp[order(lp)]
lrSort <- lr[order(lp)]
lowSort <- lowLims[order(lp)]
upSort <- upLims[order(lp)]

# plot(lpSort, ySort, type = "l")
# polygon(c(lpSort, rev(lpSort)), c(ySort + lowSort, rev(ySort + upSort)), col = gray(0.5), border = NULL, lwd = 0.01)

distep <- as.numeric(unlist(strsplit(jj, "[.]")))

framePol <- data.frame(dist = distep[seq(1, length(distep), 2)],
                       step = distep[-seq(1, length(distep), 2)],
                       yp = ySort,
                       lp = lpSort,
                       lr = lrSort,
                       lowLim = ySort + lowSort,
                       upLim = ySort + upSort)

par(mfcol=c(9,10))

for (i in sort(unique(framePol$dist))) {
  
  for (j in sort(unique(framePol$step))) {
    
    ySort <- framePol[framePol$dist == i & framePol$step == j, "yp"]
    lpSort <- framePol[framePol$dist == i & framePol$step == j, "lp"]
    lrSort <- framePol[framePol$dist == i & framePol$step == j, "lr"]
    lowSort <- framePol[framePol$dist == i & framePol$step == j, "lowLim"]
    upSort <- framePol[framePol$dist == i & framePol$step == j, "upLim"]
    
    # ySort <- frollmean(framePol[framePol$dist == i & framePol$step == j, "yp"], 5)
    # lpSort <- frollmean(framePol[framePol$dist == i & framePol$step == j, "lp"], 5)
    # lrSort <- frollmean(framePol[framePol$dist == i & framePol$step == j, "lr"], 5)
    # lowSort <- frollmean(framePol[framePol$dist == i & framePol$step == j, "lowLim"], 5)
    # upSort <- frollmean(framePol[framePol$dist == i & framePol$step == j, "upLim"], 5)
    
    # lowesLP <- lowess(lpSort, ySort)
    # lowesLR <- lowess(lrSort, ySort)
    lowesPIlow <- lowess(lpSort, lowSort)
    lowesPIup <- lowess(lpSort, upSort)
    
    plot(lpSort, ySort, ylim = c(0,350000), pch = "",
         main = paste(i, j, sep="."), ylab = "LR", xlab = "LP", cex.main = 2.5, cex.axis = 2.5)
    polygon(c(lowesPIlow[[1]], rev(lowesPIup[[1]])), c(lowesPIlow[[2]], rev(lowesPIup[[2]])), col = gray(0.6), border = NULL, lwd = 0.01)
    # polygon(c(lpSort, rev(lpSort)), c(lowSort, rev(upSort)), col = gray(0.6), border = NULL, lwd = 0.01)

    lines(lpSort, ySort, type = "l")
    lines(lpSort, lrSort, col = "red")
    
  }
  
}
```

###PI 2

```{r, fig.width=15, fig.height=15}
piLimits <- list()
cv1 <- cvList[[1]]

for (i in 1:90) {
  piLimits[[i]] <- lapply(cv1[[i]], function(x) quantile(x, c(0.025, 0.975)))
}

names(piLimits) <- names(cv1)
limVec <- unlist(unlist(piLimits))
oddPos <- seq(1, length(limVec), 2)

lowLims <- limVec[oddPos]
upLims <- limVec[-oddPos]
yp <- unlist(y.cv)

testDat <- cvRand[[1]]$testSet
testOrd <- testDat[order(testDat$dist, testDat$step), ]
lp <-testOrd$LP
lr <- testOrd$LR

  
ySort <- sort(yp)
lpSort <- lp[order(yp)]
lrSort <- lr[order(yp)]
lowSort <- lowLims[order(yp)]
upSort <- upLims[order(yp)]

# plot(lpSort, ySort, type = "l")
# polygon(c(lpSort, rev(lpSort)), c(ySort + lowSort, rev(ySort + upSort)), col = gray(0.5), border = NULL, lwd = 0.01)

distep <- as.numeric(unlist(strsplit(jj, "[.]")))

framePol <- data.frame(dist = distep[seq(1, length(distep), 2)],
                       step = distep[-seq(1, length(distep), 2)],
                       yp = ySort,
                       lp = lpSort,
                       lr = lrSort,
                       lowLim = ySort + lowSort,
                       upLim = ySort + upSort)

coul <- viridisLite::viridis(10)

par(mfrow=c(3,3))

for (j in sort(unique(framePol$step))) {
  
  colCount <- 1
  xLim <- c(min(framePol[framePol$step == j, "lp"]), max(framePol[framePol$step == j, "lp"]))
  
  for (i in sort(unique(framePol$dist))) {
    
    ySort <- framePol[framePol$dist == i & framePol$step == j, "yp"]
    lpSort <- framePol[framePol$dist == i & framePol$step == j, "lp"]
    lrSort <- framePol[framePol$dist == i & framePol$step == j, "lr"]
    lowSort <- framePol[framePol$dist == i & framePol$step == j, "lowLim"]
    upSort <- framePol[framePol$dist == i & framePol$step == j, "upLim"]
    
        
    lowesPIlow <- lowess(lpSort, lowSort)
    lowesPIup <- lowess(lpSort, upSort)
    
    if (colCount == 1) {
      
      plot(lpSort, ySort, ylim = c(0,350000), xlim = xLim, pch = "",
           main = j, ylab = "LR", xlab = "LP", cex.main = 2.5, cex.axis = 1.5, cex.lab = 1.5)
      
      polygon(c(lowesPIlow[[1]], rev(lowesPIup[[1]])), c(lowesPIlow[[2]], rev(lowesPIup[[2]])), col = alpha(coul[colCount], 0.5), border = NULL)

      lines(lpSort, ySort, type = "l")
      lines(lpSort, lrSort, col = "red")
      
      
    } else {
      
      polygon(c(lowesPIlow[[1]], rev(lowesPIup[[1]])), c(lowesPIlow[[2]], rev(lowesPIup[[2]])), col = alpha(coul[colCount], 0.5), border = NULL)

      lines(lpSort, ySort, type = "l")
      lines(lpSort, lrSort, col = "red")
      
      
    }
    
   abline(h = 23697.93, lty = "dashed", col ="gray")
   colCount <- colCount + 1

    
  }
  
}
```


##Plotting for prueba

```{r}
estimList <- list()

for (h in 1:10) {
  lrSet <- cvRand[[h]][["testSet"]]
  lrOrder <- lrSet[with(lrSet, order(dist, step)), ]
  ySet <- y.cv[[h]]
  
  estimFrame <- data.frame(dist = lrOrder$dist, step = lrOrder$step, LR = lrOrder$LR, pred = ySet)
  
  estimList[[h]] <- estimFrame
}

estimTotal <- do.call(rbind, estimList)
estimTotal$pe <- ((estimTotal$LR - estimTotal$pred)/estimTotal$LR)*100
estimTotal$abspe <- abs(estimTotal$pe)
estimTotal$interact <- interaction(estimFrame$dist, estimFrame$step, lex.order = TRUE)

estimMean <- aggregate(abspe ~ dist + step, estimTotal, mean)
```

###Mean

```{r, fig.height=10, fig.width=20}
meanFrame <- estimMean
meanFrame$neurType <- rep("Type1", dim(meanFrame)[1])

# 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")) {
  
  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()
```

###Density

```{r,  fig.width=5, fig.height=20}
ggplot(estimTotal, 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))
  
```

###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))
axProb <- quantType(estimTotal, probSeq)
# Save 
# saveRDS(axProb, "errorProbs_CVleaveout.rds")
```

######Plot heatmap

```{r, width=15, height=10}
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))

probFrame <- data.frame(error = rep(probSeq, 90), 
                        neurType = rep("Type1", length(probSeq)*90),
                        dist = rep(unique(allData$dist), each = length(probSeq)*9),
                        step = rep(rep(unique(allData$step), each = length(probSeq)), 10),
                        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")) {
  
  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")) {

  png(filename=paste("prueba_", 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()

}


```

##Plotting for preds

```{r}
# predsData <- readRDS("D:/Sergio/Cremoto Drive/stereo_boot10CVPI_LRLP_predReal100_Type1.rds")
# jjMse <-  jj %>%
#      group_by(dist, step) %>% 
#      summarise(val=mltools::mse(pred, LR))

testSets <- lapply(predsData, function(x) do.call(rbind, x$testSets))
testJoin <- do.call(rbind, testSets)
testJoin$pe <- ((testJoin$LR - testJoin$pred)/testJoin$LR)*100
testJoin$abspe <- abs(testJoin$pe)
testJoin$interact <- interaction(testJoin$dist, testJoin$step, lex.order = TRUE)
testMean <- aggregate(abspe ~ dist + step, testJoin, mean)
```

```{r}
# predsData <- readRDS("D:/Sergio/Cremoto Drive/stereo_boot10CVPI_LRLP_predReal100_Type1.rds")
# jjMse <-  jj %>%
#      group_by(dist, step) %>% 
#      summarise(val=mltools::mse(pred, LR))

testSets <- lapply(predsData, function(x) do.call(rbind, x$testSets))
testPE <- lapply(testSets, function(x) x$pe <- ((x$LR - x$pred)/x$LR)*100)
testJoin <- do.call(rbind, testSets)
testJoin$pe <- ((testJoin$LR - testJoin$pred)/testJoin$LR)*100
testJoin$abspe <- abs(testJoin$pe)
testJoin$interact <- interaction(testJoin$dist, testJoin$step, lex.order = TRUE)
testMean <- aggregate(abspe ~ dist + step, testJoin, mean)
```

###Mean

```{r, fig.height=10, fig.width=20}
meanFrame <- testMean
meanFrame$neurType <- rep("Type1", dim(meanFrame)[1])

# 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")) {
  
  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()
```

###Density

```{r,  fig.width=5, fig.height=20}
ggplot(testJoin, 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))
  
```

###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$pe, j) - quantInv(dataProb$pe, -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))
axProb <- quantType(testJoin, probSeq)
# Save 
# saveRDS(axProb, "errorProbs_CVleaveout.rds")
```

######Plot heatmap

```{r, width=15, height=10}
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))

probFrame <- data.frame(error = rep(probSeq, 90), 
                        neurType = rep("Type1", length(probSeq)*90),
                        dist = rep(unique(allData$dist), each = length(probSeq)*9),
                        step = rep(rep(unique(allData$step), each = length(probSeq)), 10),
                        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")) {
  
  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")) {

  png(filename=paste("prueba_", 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()

}


```


##Plotting for predMean

```{r}
predsData <- readRDS("D:/Sergio/Cremoto Drive/stereo_boot10CVPI_LRLP_predRealMean_GEE_Type1.rds")
# jjMse <-  jj %>%
#      group_by(dist, step) %>% 
#      summarise(val=mltools::mse(pred, LR))

testJoin <- do.call(rbind, predsData[[1]])
testJoin$interact <- interaction(testJoin$dist, testJoin$step, lex.order = TRUE)
testMean <- aggregate(abspe ~ dist + step, testJoin, mean)

peJoin <- do.call(rbind, predsData[[2]])
peJoin$interact <- interaction(peJoin$dist, peJoin$step, lex.order = TRUE)
```

###Mean

```{r, fig.height=10, fig.width=20}
meanFrame <- testMean
meanFrame$neurType <- rep("Type1", dim(meanFrame)[1])

# 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")) {
  
  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()
```

###Density

```{r,  fig.width=5, fig.height=20}
# https://cran.r-project.org/web/packages/ggridges/vignettes/introduction.html

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

```

```{r, fig.width=50, fig.height=45}
histogram(~ pe | factor(dist) + factor(step), data = peJoin)
```


###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$pe, j) - quantInv(dataProb$pe, -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))
axProb <- quantType(peJoin, probSeq)
# Save 
# saveRDS(axProb, "errorProbs_CVleaveout.rds")
```

######Plot heatmap

```{r, width=15, height=10}
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))

probFrame <- data.frame(error = rep(probSeq, 90), 
                        neurType = rep("Type1", length(probSeq)*90),
                        dist = rep(unique(allData$dist), each = length(probSeq)*9),
                        step = rep(rep(unique(allData$step), each = length(probSeq)), 10),
                        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")) {
  
  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")) {

  png(filename=paste("prueba_", 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()

}


```


##Plotting for CVloo

```{r}
predsData <- readRDS("D:/Sergio/Cremoto Drive/stereo_bootCVlooPI_LRLP_predReal_Type1.rds")
# jjMse <-  jj %>%
#      group_by(dist, step) %>% 
#      summarise(val=mltools::mse(pred, LR))

testJoin <- predsData
testJoin$interact <- interaction(testJoin$dist, testJoin$step, lex.order = TRUE)
testMean <- aggregate(abspe ~ dist + step, testJoin, mean)

peJoin <- predsData
peJoin$interact <- interaction(peJoin$dist, peJoin$step, lex.order = TRUE)
```

###Mean

```{r, fig.height=10, fig.width=20}
meanFrame <- testMean
meanFrame$neurType <- rep("Type1", dim(meanFrame)[1])

# 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")) {
  
  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()
```

###Density

```{r,  fig.width=5, fig.height=20}
# https://cran.r-project.org/web/packages/ggridges/vignettes/introduction.html

print(ggplot(peJoin, aes(x = `pe`, y = `interact`, fill = ifelse(..x.. >=-5 & ..x.. <=5, ">= -5", "< 5"))) +
  scale_fill_manual(values = c("gray70", "coral2"), name = NULL) +
  stat_density_ridges(geom = "density_ridges_gradient") +
  theme_ridges() + 
  theme(legend.position = "none") + 
  xlab("% Error") + 
  xlim(c(-40,40)))

```

```{r, fig.width=50, fig.height=45}
histogram(~ pe | factor(dist) + factor(step), data = peJoin)
```


###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$pe, j) - quantInv(dataProb$pe, -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))
axProb <- quantType(peJoin, probSeq)
# Save 
# saveRDS(axProb, "errorProbs_CVleaveout.rds")
```

######Plot heatmap

```{r, width=15, height=10}
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))

probFrame <- data.frame(error = rep(probSeq, 90), 
                        neurType = rep("Type1", length(probSeq)*90),
                        dist = rep(unique(allData$dist), each = length(probSeq)*9),
                        step = rep(rep(unique(allData$step), each = length(probSeq)), 10),
                        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")) {
  
  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")) {

  png(filename=paste("prueba_", 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()

}


```