Rubio_etal_2021/Code/Statistical Analysis/modelFree_spheres_MSE_figures.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(grid)
library(gridExtra)
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("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)
```


## Estimate  MSE

```{r}
msediamsep <- by(allData, allData$interact2, function(x) sqrt(mse(x$LRe, x$LR)))
# msediamsep <- by(allData, allData$interact2, function(x) cor(x$LRe, x$LR)^2)
```

###RMSE ggplot (smooth bands)

```{r}
#######################
# Organize data frame #
#######################

# Store MSE in data frame
mseFrame <- data.frame(mse = as.numeric(msediamsep), 
                       x = rep(1:81, 4), 
                       step = rep(rep(seq(70,150,10), each = 9), 4), 
                       diams = rep(rep(seq(10,50,5), 9), 4),
                       axType = rep(c("Type1", "Type2", "Type3", "Type4"), each = 81))

# Find max and min for Types1-3
min123 <- min(mseFrame[mseFrame$axType != "Type4", "mse"])
max123 <- max(mseFrame[mseFrame$axType != "Type4", "mse"])


# Find max and min for Type4
min4 <- min(mseFrame[mseFrame$axType == "Type4", "mse"])
max4 <- max(mseFrame[mseFrame$axType == "Type4", "mse"])

# Normalize
mseFrame$maxMSE <- c(rep(max123, 81*3), rep(max4, 81))
mseFrame$mseNorm <- mseFrame$mse/mseFrame$maxMSE

# Color palette
pal <- wes_palette("Zissou1", 9, type = "continuous")
rmseList <- list()

# Find max and min for plotting
min123plot <- min(mseFrame[mseFrame$axType != "Type4", "mseNorm"])
max123plot <- max(mseFrame[mseFrame$axType != "Type4", "mseNorm"])
min4plot <- min(mseFrame[mseFrame$axType == "Type4", "mseNorm"])
max4plot <- max(mseFrame[mseFrame$axType == "Type4", "mseNorm"])

# Color palette
pal <- wes_palette("Zissou1", 9, type = "continuous")
rmseList <- list()

####################
# RMSE Types 1,2,3 #
####################

mseType <- mseFrame[mseFrame$axType != "Type4", ]

mse123 <- ggplot(data=mseType, aes(x=x, y=mseNorm, group=step)) +
  # geom_line() +
  facet_grid(. ~ axType) +
  geom_smooth(colour = "black", method = "loess", span = 1.2, size = 0.5) +
  geom_point(aes(colour = diams), size = 3) +
  # geom_point(shape = 1, colour = "black") +
  scale_x_continuous(breaks=seq(4.5,79.5,9), labels = seq(70,150,10)) +
  scale_colour_gradientn(name = "Diameter",
                         colours = pal,
                         breaks = seq(10,50,5)) +
  # 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"),
        strip.text = element_text(size = 20),
        plot.title = element_text(hjust = 0.5, size = 22),
        axis.text.x = element_text(size = 15, color = "black"),
        axis.text.y = element_text(size = 15, color = "black"),
        axis.title.x = element_text(size = 18),
        axis.title.y = element_text(size = 18),
        panel.background = element_rect(colour = "black",size = 1)) +
  xlab("Step") +
  ylab("Model error")

###############
# RMSE Type 4 #
###############

mseType <- mseFrame[mseFrame$axType == "Type4", ]
# mseType$axType <- factor(mseFrame$axType, levels = c("Type4", "Type1", "Type2", "Type3"))

mse4 <- ggplot(data=mseType, aes(x=x, y=mseNorm, group=step)) +
  # geom_line() +
  facet_grid(. ~ axType) +
  geom_smooth(colour = "black", method = "loess", span = 1.2, size = 0.5) +
  geom_point(aes(colour = diams), size = 3) +
  # geom_point(shape = 1, colour = "black") +
  scale_x_continuous(breaks=seq(4.5,79.5,9), labels = seq(70,150,10)) +
  scale_colour_gradientn(name = "Diameter",
                         colours = pal,
                         breaks = seq(10,50,5)) +
  # 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"),
        strip.text = element_text(size = 20),
        plot.title = element_text(hjust = 0.5, size = 20),
        axis.text.x = element_text(size = 15, color = "black"),
        axis.text.y = element_text(size = 15, color = "black"),
        axis.title.x = element_text(size = 18),
        axis.title.y = element_text(size = 18),
        panel.background = element_rect(colour = "black",size = 1)) +
  xlab("Step") +
  ylab("Model error")

###########################
# EXAMPLES Types 1,2,3, 4 #
###########################

#List to store
exList <- list()

# Type 1
for (i in unique(allData$neurType)) {
  
  data2d <- allData[allData$neurType == i & (allData$interact == "70.10" | allData$interact == "150.50"), ]
  dataReal <- data.frame(LR = data2d$LR, LRe = data2d$LR)
  
  if (i == "Type1" | i == "Type4") {
    
     exList[[i]] <- ggplot(data=data2d, aes(x=LRe, y=LR)) +
    facet_grid(. ~ interact) +
    geom_point(colour = "black", size = 1) + 
    # geom_point(shape = 1, colour = "black") +
    geom_point(data = dataReal, colour = "red", size = 1) + 
    # geom_point(data = dataReal, shape = 1, colour = "black") +
    theme_bw() + 
    theme(panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          strip.text = element_text(size = 20),
          plot.title = element_text(hjust = 0.5, size = 12),
          axis.title.x=element_text(size = 14),
          axis.text.x=element_blank(),
          axis.ticks.x=element_blank(),
          axis.title.y=element_text(size = 14),
          axis.text.y=element_blank(),
          axis.ticks.y=element_blank(), 
          panel.background = element_rect(colour = "black",size = 1),
          plot.margin=unit(c(0.2,-0.2,0.2,-0.2), "cm")) +
    xlab("Predicted length") +
    ylab("Real Length")
     
  } else {
    
     exList[[i]] <- ggplot(data=data2d, aes(x=LRe, y=LR)) +
    facet_grid(. ~ interact) +
    geom_point(colour = "black", size = 1) + 
    # geom_point(shape = 1, colour = "black") +
    geom_point(data = dataReal, colour = "red", size = 1) + 
    # geom_point(data = dataReal, shape = 1, colour = "black") +
    theme_bw() + 
    theme(panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          strip.text = element_text(size = 20),
          plot.title = element_text(hjust = 0.5, size = 12),
          axis.title.x=element_text(size = 14),
          axis.text.x=element_blank(),
          axis.ticks.x=element_blank(),
          axis.title.y=element_blank(),
          axis.text.y=element_blank(),
          axis.ticks.y=element_blank(), 
          panel.background = element_rect(colour = "black",size = 1),
          plot.margin=unit(c(0.2,-0.2,0.2,-0.2), "cm")) +
    xlab("Predicted length")
  
     
  }
 
}


```

###Arrange

```{r, fig.width=15, fig.height=8}
# https://cran.r-project.org/web/packages/gridExtra/vignettes/arrangeGrob.html

setwd("EsferasFigures/")
png(filename=paste("rmse123_arranged.png", sep=""), type="cairo",
      units="in", width=15, height=8, pointsize=12,
      res=300)

gridExtra::grid.arrange(mse123, exList[[1]], exList[[2]], exList[[3]], ncol = 7, nrow = 2,
             layout_matrix = rbind(rep(1, 7), c(NA,2,NA,3,NA,4,NA)),
             widths = c(0.8, 5.8, 0.28, 5.5, 0.28, 5.5, 1.375),
             heights = c(6, 3))

dev.off()
```

```{r, fig.width=20, fig.height=8}
# https://cran.r-project.org/web/packages/gridExtra/vignettes/arrangeGrob.html

setwd("EsferasFigures/")
png(filename=paste("rmseFull_arranged.png", sep=""), type="cairo",
      units="in", width=20, height=8, pointsize=12,
      res=300)

gridExtra::grid.arrange(mse123, exList[[1]], exList[[2]], exList[[3]], mse4, exList[[4]], ncol = 11, nrow = 2,
             layout_matrix = rbind(c(rep(1, 7), NA, 5, 5, 5), 
                                   c(NA, 2, NA, 3, NA, 4, NA, NA, NA, 6, NA)),
             widths = c(0.95, 5.8, 0.28, 5.5, 0.28, 5.5, 1.55, 1, 0.9, 5.85, 1.5),
             heights = c(6, 3))

dev.off()

```

###Save mse

UNCOMMENT IF ARRANGE DOES NOT WORK

```{r, fig.width=15, fig.height=5}
# # Types 1,2,3
# setwd("EstereoFigures/")
# png(filename=paste("rmse123.png", sep=""), type="cairo",
#       units="in", width=15, height=5, pointsize=12,
#       res=300)
# mse123
# dev.off()
# 
# # Type 4
# png(filename=paste("rmse4.png", sep=""), type="cairo",
#       units="in", width=6, height=5, pointsize=12,
#       res=300)
# mse4
# dev.off()
# 
# # Examples
# png(filename=paste("examplesRMSE_Type1.png", sep=""), type="cairo",
#     units="in", width=3.5, height=2.25, pointsize=12,
#     res=300)
# exList[["Type1"]]
# dev.off()
# 
# png(filename=paste("examplesRMSE_Type2.png", sep=""), type="cairo",
#     units="in", width=3.5, height=2.25, pointsize=12,
#     res=300)
# exList[["Type2"]]
# dev.off()
# 
# png(filename=paste("examplesRMSE_Type3.png", sep=""), type="cairo",
#     units="in", width=3.5, height=2.25, pointsize=12,
#     res=300)
# exList[["Type3"]]
# dev.off()
# 
# png(filename=paste("examplesRMSE_Type4.png", sep=""), type="cairo",
#     units="in", width=3.5, height=2.25, pointsize=12,
#     res=300)
# exList[["Type4"]]
# dev.off()

```

```{r, fig.width=3.5, fig.height=2.25}
setwd("EstereoFigures/")



```


## Estimate  R2

```{r}
msediamsep <- by(allData, allData$interact2, function(x) cor(x$LRe, x$LR)^2)
```

###R2 ggplot (smooth bands)

```{r}
#######################
# Organize data frame #
#######################

# Store MSE in data frame
mseFrame <- data.frame(mse = as.numeric(msediamsep), 
                       x = rep(1:81, 4), 
                       step = rep(rep(seq(70,150,10), each = 9), 4), 
                       diams = rep(rep(seq(10,50,5), 9), 4),
                       axType = rep(c("Type1", "Type2", "Type3", "Type4"), each = 81))

# Find max and min for Types1-3
min123 <- min(mseFrame[mseFrame$axType != "Type4", "mse"])
max123 <- max(mseFrame[mseFrame$axType != "Type4", "mse"])


# Find max and min for Type4
min4 <- min(mseFrame[mseFrame$axType == "Type4", "mse"])
max4 <- max(mseFrame[mseFrame$axType == "Type4", "mse"])

# Normalize
mseFrame$maxMSE <- c(rep(max123, 81*3), rep(max4, 81))
mseFrame$mseNorm <- mseFrame$mse/mseFrame$maxMSE

# Color palette
pal <- wes_palette("Zissou1", 9, type = "continuous")
rmseList <- list()

# Find max and min for plotting
min123plot <- min(mseFrame[mseFrame$axType != "Type4", "mseNorm"])
max123plot <- max(mseFrame[mseFrame$axType != "Type4", "mseNorm"])
min4plot <- min(mseFrame[mseFrame$axType == "Type4", "mseNorm"])
max4plot <- max(mseFrame[mseFrame$axType == "Type4", "mseNorm"])

# Color palette
pal <- wes_palette("Zissou1", 9, type = "continuous")
rmseList <- list()

####################
# RMSE Types 1,2,3 #
####################

mseType <- mseFrame[mseFrame$axType != "Type4", ]

mse123 <- ggplot(data=mseType, aes(x=x, y=mseNorm, group=step)) +
  # geom_line() +
  facet_grid(. ~ axType) +
  geom_smooth(colour = "black", method = "loess", span = 1.2, size = 0.5) +
  geom_point(aes(colour = diams), size = 3) +
  # geom_point(shape = 1, colour = "black") +
  scale_x_continuous(breaks=seq(4.5,79.5,9), labels = seq(70,150,10)) +
  scale_colour_gradientn(name = "Diameter",
                         colours = pal,
                         breaks = seq(10,50,5)) +
  # 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"),
        strip.text = element_text(size = 20),
        plot.title = element_text(hjust = 0.5, size = 22),
        axis.text.x = element_text(size = 15, color = "black"),
        axis.text.y = element_text(size = 15, color = "black"),
        axis.title.x = element_text(size = 18),
        axis.title.y = element_text(size = 18),
        panel.background = element_rect(colour = "black",size = 1)) +
  xlab("Step") +
  ylab("R2")

###############
# RMSE Type 4 #
###############

mseType <- mseFrame[mseFrame$axType == "Type4", ]
# mseType$axType <- factor(mseFrame$axType, levels = c("Type4", "Type1", "Type2", "Type3"))

mse4 <- ggplot(data=mseType, aes(x=x, y=mseNorm, group=step)) +
  # geom_line() +
  facet_grid(. ~ axType) +
  geom_smooth(colour = "black", method = "loess", span = 1.2, size = 0.5) +
  geom_point(aes(colour = diams), size = 3) +
  # geom_point(shape = 1, colour = "black") +
  scale_x_continuous(breaks=seq(4.5,79.5,9), labels = seq(70,150,10)) +
  scale_colour_gradientn(name = "Diameter",
                         colours = pal,
                         breaks = seq(10,50,5)) +
  # 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"),
        strip.text = element_text(size = 20),
        plot.title = element_text(hjust = 0.5, size = 20),
        axis.text.x = element_text(size = 15, color = "black"),
        axis.text.y = element_text(size = 15, color = "black"),
        axis.title.x = element_text(size = 18),
        axis.title.y = element_text(size = 18),
        panel.background = element_rect(colour = "black",size = 1)) +
  xlab("Step") +
  ylab("R2")

###########################
# EXAMPLES Types 1,2,3, 4 #
###########################

#List to store
exList <- list()

# Type 1
for (i in unique(allData$neurType)) {
  
  data2d <- allData[allData$neurType == i & (allData$interact == "70.10" | allData$interact == "150.50"), ]
  dataReal <- data.frame(LR = data2d$LR, LRe = data2d$LR)
  
  if (i == "Type1" | i == "Type4") {
    
     exList[[i]] <- ggplot(data=data2d, aes(x=LRe, y=LR)) +
    facet_grid(. ~ interact) +
    geom_point(colour = "black", size = 1) + 
    # geom_point(shape = 1, colour = "black") +
    geom_point(data = dataReal, colour = "red", size = 1) + 
    # geom_point(data = dataReal, shape = 1, colour = "black") +
    theme_bw() + 
    theme(panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          strip.text = element_text(size = 20),
          plot.title = element_text(hjust = 0.5, size = 12),
          axis.title.x=element_text(size = 14),
          axis.text.x=element_blank(),
          axis.ticks.x=element_blank(),
          axis.title.y=element_text(size = 14),
          axis.text.y=element_blank(),
          axis.ticks.y=element_blank(), 
          panel.background = element_rect(colour = "black",size = 1),
          plot.margin=unit(c(0.2,-0.2,0.2,-0.2), "cm")) +
    xlab("Predicted length") +
    ylab("Real Length")
     
  } else {
    
     exList[[i]] <- ggplot(data=data2d, aes(x=LRe, y=LR)) +
    facet_grid(. ~ interact) +
    geom_point(colour = "black", size = 1) + 
    # geom_point(shape = 1, colour = "black") +
    geom_point(data = dataReal, colour = "red", size = 1) + 
    # geom_point(data = dataReal, shape = 1, colour = "black") +
    theme_bw() + 
    theme(panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          strip.text = element_text(size = 20),
          plot.title = element_text(hjust = 0.5, size = 12),
          axis.title.x=element_text(size = 14),
          axis.text.x=element_blank(),
          axis.ticks.x=element_blank(),
          axis.title.y=element_blank(),
          axis.text.y=element_blank(),
          axis.ticks.y=element_blank(), 
          panel.background = element_rect(colour = "black",size = 1),
          plot.margin=unit(c(0.2,-0.2,0.2,-0.2), "cm")) +
    xlab("Predicted length")
  
     
  }
 
}

```

###Arrange

```{r, fig.width=20, fig.height=8}
# https://cran.r-project.org/web/packages/gridExtra/vignettes/arrangeGrob.html

setwd("EsferasFigures/")
png(filename=paste("r2_Full_arranged.png", sep=""), type="cairo",
      units="in", width=20, height=8, pointsize=12,
      res=300)

gridExtra::grid.arrange(mse123, exList[[1]], exList[[2]], exList[[3]], mse4, exList[[4]], ncol = 11, nrow = 2,
             layout_matrix = rbind(c(rep(1, 7), NA, 5, 5, 5), 
                                   c(NA, 2, NA, 3, NA, 4, NA, NA, NA, 6, NA)),
             widths = c(0.95, 5.8, 0.28, 5.5, 0.28, 5.5, 1.55, 1, 0.9, 5.85, 1.5),
             heights = c(6, 3))

dev.off()

```

###Save R2

UNCOMMENT IF ARRANGE DOES NOT WORK

```{r, fig.width=15, fig.height=5}
# # Types 1,2,3
# setwd("EstereoFigures/")
# png(filename=paste("r2_123.png", sep=""), type="cairo",
#       units="in", width=15, height=5, pointsize=12,
#       res=300)
# mse123
# dev.off()
# 
# # Type 4
# png(filename=paste("r2_4.png", sep=""), type="cairo",
#       units="in", width=6, height=5, pointsize=12,
#       res=300)
# mse4
# dev.off()

```


##Mean error (Type123 vs 4)

```{r, fig.height=10, fig.width=20}
# Estimate mean residual error
meanFrame <- aggregate(error ~ diams + step + neurType, allData, mean)

# Library
library(latticeExtra)
# coul <- pals::parula(10000)
coul <- pals::jet(100000)

# Store heatmaps
heatList <- list()
smoothList <- list()

# Limits
min123 <- min(meanFrame[meanFrame$neurType != "Type4", "error"])
max123 <- max(meanFrame[meanFrame$neurType != "Type4", "error"])

# Type 1
smoothList[["Type1"]] <- levelplot(error ~ diams * step, data = meanFrame[meanFrame$neurType == "Type1", ], col.regions = coul,
                               scales = list(y = list(at = unique(allData$step), labels = unique(allData$step), cex = 1.5),
                                             x = list(at = unique(allData$diams), labels = unique(allData$diams)), cex = 1.5),
                               ylim = c(150,70),
                               xlim = c(10,50),
                               ylab = list(label = "Step", cex = 1.5),
                               xlab = list(label = "Diameter", cex = 1.5),
                               main = list(label = "Type1", cex = 1.5),
                               at = seq(min123,
                                       max123,
                                       0.05),
                               cuts = 100, 
                               colorkey = FALSE,
                               panel = panel.levelplot.points, cex = 0) + 
    layer_(panel.2dsmoother(..., n = 200))

# Type 2
smoothList[["Type2"]] <- levelplot(error ~ diams * step, data = meanFrame[meanFrame$neurType == "Type2", ], col.regions = coul,
                             scales = list(y = list(at = unique(allData$step), labels = NULL, cex = 1.5),
                                           x = list(at = unique(allData$diams), labels = unique(allData$diams)), cex = 1.5),
                             ylim = c(150,70),
                             xlim = c(10,50),
                             ylab="",
                             xlab = list(label = "Diameter", cex = 1.5),
                             main = list(label = "Type2", cex = 1.5),
                             at = seq(min123,
                                      max123,
                                      0.025),
                             cuts = 100, 
                             colorkey = FALSE,
                             panel = panel.levelplot.points, cex = 0) + 
  layer_(panel.2dsmoother(..., n = 200))

# Type 3
smoothList[["Type3"]] <- levelplot(error ~ diams * step, data = meanFrame[meanFrame$neurType == "Type3", ], col.regions = coul,
                             scales = list(y = list(at = unique(allData$step), labels = NULL, cex = 1.5),
                                           x = list(at = unique(allData$diams), labels = unique(allData$diams)), cex = 1.5),
                             ylim = c(150,70),
                             xlim = c(10,50),
                             ylab="",
                             xlab = list(label = "Diameter", cex = 1.5),
                             main = list(label = "Type3", cex = 1.5),
                             at = seq(min123,
                                      max123,
                                      0.025),
                             cuts = 100, 
                             colorkey = list(at = seq(min123,
                                                      max123,
                                                      0.025),
                                             cex = 1.2),
                             panel = panel.levelplot.points, cex = 0) + 
  layer_(panel.2dsmoother(..., n = 200))
  
# Type 4
smoothList[["Type4"]] <- levelplot(error ~ diams * step, data = meanFrame[meanFrame$neurType == "Type4", ], col.regions = coul,
                               scales = list(y = list(at = unique(allData$step), labels = unique(allData$step), cex = 1.5),
                                             x = list(at = unique(allData$diams), labels = unique(allData$diams)), cex = 1.5),
                               ylim = c(150,70),
                               xlim = c(10,50),
                               ylab = list(label = "Step", cex = 1.5),
                               xlab = list(label = "Diameter", cex = 1.5),
                               main = list(label = "Type4", cex = 1.5),
                               cuts = 100, 
                               colorkey = list(at = (seq(min(meanFrame[meanFrame$neurType == "Type4", "error"]),
                                                         max(meanFrame[meanFrame$neurType == "Type4", "error"]),
                                                         0.05)),
                                               cex = 1.2),
                               panel = panel.levelplot.points, cex = 0) + 
    layer_(panel.2dsmoother(..., n = 200))
```

####Arrange

```{r, fig.width=20, fig.height=5}
# Plot
setwd("EsferasFigures/")
png(filename=paste("meanError_Full_jet.png", sep=""), type="cairo",
    units="in", width=20, height=5, pointsize=12,
    res=300)

gridExtra::grid.arrange(grobs = smoothList, ncol = 5, nrow = 1, 
                        layout_matrix = rbind(c(1,2,3,NA,4)),
                        widths = c(5,4.5,5,0.5,5.5))

dev.off()
```

####Joint arrange

```{r}
coul <- pals::jet(100000)

# Types 1,2,3
m123 <- levelplot(error ~ diams * step | neurType, data = meanFrame[meanFrame$neurType != "Type4", ], col.regions = coul,
          scales = list(y = list(at = unique(allData$step), labels = unique(allData$step), cex = 1.5),
                        x = list(at = unique(allData$diams), labels = unique(allData$diams)), cex = 1.5),
          ylim = c(150,70),
          xlim = c(10,50),
          ylab = list(label = "Step", cex = 1.5),
          xlab = list(label = "Diameter", cex = 1.5),
          at = seq(min123,
                   max123,
                   0.05),
          cuts = 100, 
          colorkey = list(labels=list(cex=1.5)),
          layout = c(3, 1),
          par.settings = list(strip.background=list(col="white")),
          par.strip.text=list(cex=1.5),
          panel = panel.levelplot.points, cex = 0) + 
  layer_(panel.2dsmoother(..., n = 200))

# Type4
m4 <- levelplot(error ~ diams * step | neurType, data = meanFrame[meanFrame$neurType == "Type4", ], col.regions = coul,
          scales = list(y = list(at = unique(allData$step), labels = unique(allData$step), cex = 1.5),
                        x = list(at = unique(allData$diams), labels = unique(allData$diams), cex = 1.5, alternating = 3)),
          ylim = c(150,70),
          xlim = c(10,50),
          ylab = list(label = "Step", cex = 1.5),
          xlab = list(label = "Diameter", cex = 1.5),
          cuts = 100, 
          colorkey = list(labels=list(cex=1.5)),
          par.settings = list(strip.background=list(col="white")),
          par.strip.text=list(cex=1.5),
          panel = panel.levelplot.points, cex = 0) + 
  layer_(panel.2dsmoother(..., n = 200))

# Plot
setwd("EsferasFigures/")
png(filename=paste("meanError_Full_joint_jet.png", sep=""), type="cairo",
    units="in", width=24, height=6.1, pointsize=12,
    res=300)

gridExtra::grid.arrange(m123, m4, ncol = 3, nrow = 1, 
                        layout_matrix = rbind(c(1,NA,2)),
                        widths = c(17,1.2,7.4))

dev.off()
```

##Errors vs cumul probability

RUN ONLY ONCE, THEN LOAD RDS FILE IN PLOT CHUNK

```{r, fig.height=5, fig.width=20}
# Inverse quantile
quantInv <- function(diamsr, value) ecdf(diamsr)(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))
setwd("EsferasAnalysis/")
saveRDS(axProb, "errorProbs_RMSE.rds")
```

###Load probs

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

####Plot conditional heatmap

```{r, fig.height=15, fig.width=15}
# Color scale
coul <- viridis::magma(10000)

###############
# Types 1,2,3 #
###############

# Create data frame to plot
dataPlot <- probFrame[probFrame$neurType != "Type4" & (probFrame$error == 5 | probFrame$error == 10 | probFrame$error == 15), ]
dataPlot$condError <- factor(dataPlot$error, levels = c("15", "10", "5"))
levels(dataPlot$condError) <- c("15%", "10%", "5%")

lp123 <- levelplot(prob ~ diams * step | neurType + condError, data = dataPlot, 
                   col.regions = coul,
                   scales = list(y = list(at = unique(allData$step), labels = unique(allData$step), cex = 1.5),
                                 x = list(at = unique(allData$diams), labels = unique(allData$diams), cex = 1.5)),
                   ylim = c(150,70),
                   xlim = c(10,50),
                   cuts = 100,
                   colorkey = list(at = (seq(min(dataPlot$prob),
                                             max(dataPlot$prob),
                                             0.0025)),
                                   labels = list(cex = 1.5),
                                   title = "P(%error <= x)"),
                   ylab = list(label = "Step", cex = 2),
                   xlab = list(label = "Diameter", cex = 2),
                   par.settings = list(strip.background=list(col=c("slategray1","white"))),
                   par.strip.text=list(cex=1.5),
                   panel = panel.levelplot.points, cex = 0) + 
  layer_(panel.2dsmoother(..., n = 200))

#########
# Type4 #
#########

# Create data frame to plot
dataPlot <- probFrame[probFrame$neurType == "Type4" & (probFrame$error == 5 | probFrame$error == 10 | probFrame$error == 15), ]
dataPlot$condError <- factor(dataPlot$error, levels = c("15", "10", "5"))
levels(dataPlot$condError) <- c("15%", "10%", "5%")

lp4 <- levelplot(prob ~ diams * step | neurType + condError, data = dataPlot, 
                   col.regions = coul,
                   scales = list(y = list(at = unique(allData$step), labels = unique(allData$step), cex = 1.5),
                                 x = list(at = unique(allData$diams), labels = unique(allData$diams), cex = 1.5, alternating = 3)),
                   ylim = c(150,70),
                   xlim = c(10,50),
                   cuts = 100,
                   colorkey = list(at = (seq(min(dataPlot$prob),
                                             max(dataPlot$prob),
                                             0.0025)),
                                   labels = list(cex = 1.5),
                                   title = "P(%error <= x)"),
                   ylab = list(label = "Step", cex = 2),
                   xlab = list(label = "Diameter", cex = 2),
                   # main = list(label = paste("P(%Error <= x)", sep = ""), cex = 2),
                   par.settings = list(strip.background=list(col=c("slategray1","white"))),
                   par.strip.text=list(cex=1.5),
                   panel = panel.levelplot.points, cex = 0) + 
  layer_(panel.2dsmoother(..., n = 200))
```

####Arrange

```{r}
# Plot
setwd("EsferasFigures/")
png(filename=paste("probError_full_jet.png", sep=""), type="cairo",
    units="in", width=24, height=15, pointsize=12,
    res=300)

gridExtra::grid.arrange(lp123, lp4, ncol = 3, nrow = 1, 
                        layout_matrix = rbind(c(1,NA,2)),
                        widths = c(17,1.2,7.5))

dev.off()
```


##Boxplot error~diams: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(diams))) +
#   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 + diams, 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 = "Diameter",
    #                        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")
```