Rubio_etal_2021/Code/Statistical Analysis/modelFree_stereology_errors_CV.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)

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