Rubio_etal_2021/Code/Statistical Analysis/modelFree_stereology_MSE.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)
library(wesanderson)

options(digits = 10)

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


###############################
# Add continuous color legend #
###############################

legend.col <- function(col, lev){
  
  opar <- par
  
  n <- length(col)
  
  bx <- par("usr")
  
  box.cx <- c(bx[2] + (bx[2] - bx[1]) / 1000,
              bx[2] + (bx[2] - bx[1]) / 1000 + (bx[2] - bx[1]) / 50)
  box.cy <- c(bx[3], bx[3])
  box.sy <- (bx[4] - bx[3]) / n
  
  xx <- rep(box.cx, each = 2)
  
  par(xpd = TRUE)
  for(i in 1:n){
    
    yy <- c(box.cy[1] + (box.sy * (i - 1)),
            box.cy[1] + (box.sy * (i)),
            box.cy[1] + (box.sy * (i)),
            box.cy[1] + (box.sy * (i - 1)))
    polygon(xx, yy, col = col[i], border = col[i])
    
  }
  par(new = TRUE)
  plot(0, 0, type = "n",
       ylim = c(min(lev), max(lev)),
       yaxt = "n", ylab = "",
       xaxt = "n", xlab = "",
       frame.plot = FALSE)
  axis(side = 4, las = 2, at = lev, tick = FALSE, line = .15)
  par <- opar
}
```


##Load data

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

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

# Extract data
# errorMatrix contains 4 arrays, one per neuron type
# within each array, the dimensions corresponds to step:dist:neuron
# this means that each array has as many matrices as neuron of a certain type
dist         <- lapply(axonFiles, function(x) x$Distancias.Entre.Planos)[[1]]
step         <- lapply(axonFiles, function(x) x$Step.Lengths)[[1]]
errorMatrix  <- lapply(axonFiles, function(x) x$MATRIZ.ERROR.LENGTH)
lpMatrix     <- lapply(axonFiles, function(x) x$MATRIZ.ESTIMATED.AXON.LENGTH)
lrVals       <- lapply(axonFiles, function(x) x$AXON.REAL.LENGTH)

# Get number of neurons per type
numberTypes  <- unlist(lapply(errorMatrix, function(x) dim(x)[3]))

# Vectorizing the arrays goes as follows: neuron, dist, step
errorVector <- unlist(lapply(errorMatrix, function(x) as.vector(x)))
lpVector <- unlist(lapply(lpMatrix, function(x) as.vector(x)))
lrVector <- unlist(lapply(lrVals, function(x) as.vector(x)))

# Store in data frames
allData <- data.frame(id = rep(1:sum(numberTypes), each = 90),
                      neurType = rep(c("Type1", "Type2", "Type3", "Type4"), times = 90*numberTypes),
                      dist = rep(rep(dist, each = 9), sum(numberTypes)), 
                      step = rep(rep(step, 10), sum(numberTypes)),
                      error = abs(errorVector),
                      error2 = errorVector,
                      LRe = lpVector,
                      LR = rep(lrVector, each = 90))
allData$errorRaw <- allData$LR - allData$LRe
allData$interact <- interaction(allData$step, allData$dist, lex.order = T)
allData$interact2 <- interaction(allData$neurType, allData$step, allData$dist, lex.order = T)
```


## Estimate MSE

```{r}
mseDistep <- by(allData, allData$interact2, function(x) sqrt(mse(x$LRe, x$LR)))
```


```{r, fig.width=15}
# mseDist <- by(allData, allData$dist, function(x) sqrt(mse(x$LRe, x$LR)))
# mseStep <- by(allData, allData$step, function(x) sqrt(mse(x$LRe, x$LR)))
# mseDistep <- by(allData, allData$interact2, function(x) sqrt(mse(x$LRe, x$LR)))
# 
# par(mfrow=c(1,3))
# 
# plot(seq(3,30,3), mseDist, 
#      xaxt = "n", xlab = "Dist", ylab = "RMSE", 
#      pch = 21, bg = "firebrick", col = NA, 
#      main = "Test Error ~ Distancia",
#      cex = 1.5, cex.axis = 1.5, cex.main = 2, cex.lab = 1.5)
# lines(seq(3,30,3), mseDist)
# axis(side=1, at=seq(3,30,3), cex.axis = 1.5)
# grid()
# 
# plot(seq(70,150,10), mseStep, 
#      xaxt = "n", xlab = "Step", ylab = "RMSE", 
#      pch = 21, bg = "firebrick", col = NA, 
#      main = "Test Error ~ Step",
#      cex = 1.5, cex.axis = 1.5, cex.main = 2, cex.lab = 1.5)
# lines(seq(70,150,10), mseStep)
# axis(side=1, at=seq(70,150,10), cex.axis = 1.5)
# grid()
# 
# plot(1:90, mseDistep, 
#      xaxt = "n", xlab = "Step", ylab = "RMSE", 
#      pch = 21, bg = rep(wes_palette("Zissou1", 10, type = "continuous"), 9), col = NA, 
#      main = "Test Error ~ Step:Dist",
#      cex = 1.5, cex.axis = 1.5, cex.main = 2, cex.lab = 1.5)
# 
# for (i in seq(0,90,10)) {
#   lines((i+1):(i+10), mseDistep[(i+1):(i+10)])
# }
# 
# axis(side=1, at=seq(5,90,10), las = 1, srt = 60, labels = seq(70, 150, 10), cex.axis = 1.5)
# abline(v=seq(0,90,10), lty="dashed", col = "gray")
# 
# legend.col(wes_palette("Zissou1", 10, type = "continuous"), seq(3,30,3))
```

###RMSE ggplot (dashed lines)

```{r, fig.height=5, fig.width=20}
# Store MSE in data frame
mseFrame <- data.frame(mse = as.numeric(mseDistep), 
                       x = rep(1:90, 4), 
                       step = rep(rep(seq(70,150,10), each = 10), 4), 
                       dist = rep(rep(seq(3,30,3), 9), 4),
                       axType = rep(c("Type1", "Type2", "Type3", "Type4"), each = 90))

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

# Plot and store in list
for (i in unique(mseFrame$axType)) {
  
  mseType <- mseFrame[mseFrame$axType == i, ]
  
  rmseList[[i]] <- ggplot(data=mseType, aes(x=x, y=mse, group=step)) +
    geom_line() +
    geom_point(aes(colour = dist)) + 
    scale_x_continuous(breaks=seq(5,90,10), labels = seq(70,150,10)) +
    scale_y_continuous(breaks=seq(2000, 16000, 2000), limits = c(2000,16000)) +
    scale_colour_gradientn(name = "Distance",
                           colours = pal,
                           breaks = seq(3,30,3)) +
    # geom_vline(xintercept=seq(0,90,10), colour = "gray", linetype="dashed", ) +
    # geom_smooth() +
    theme_bw() + 
    theme(panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          legend.key.size = unit(1.9, "cm"),
          legend.key.width = unit(0.5,"cm"),
          plot.title = element_text(hjust = 0.5, size = 22),
          axis.text.x = element_text(size = 13, color = "black"),
          axis.text.y = element_text(size = 13, color = "black"),
          axis.title.x = element_text(size = 16),
          axis.title.y = element_text(size = 16),
          panel.background = element_rect(colour = "black",size = 1)) +
    xlab("Step") +
    ylab("RMSE") +
    ggtitle(i)
  
}

# Save as png
png(filename=paste("rmseEstereo2.png", sep=""), type="cairo",
    units="in", width=22, height=5, pointsize=12,
    res=300)

gridExtra::grid.arrange(grobs = rmseList, ncol = 4, nrow = 1)

dev.off()
```

###RMSE ggplot (smooth bands)

```{r, fig.height=5, fig.width=20}
# Store MSE in data frame
mseFrame <- data.frame(mse = as.numeric(mseDistep), 
                       x = rep(1:90, 4), 
                       step = rep(rep(seq(70,150,10), each = 10), 4), 
                       dist = rep(rep(seq(3,30,3), 9), 4),
                       axType = rep(c("Type1", "Type2", "Type3", "Type4"), each = 90))

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

# Plot and store in list
for (i in unique(mseFrame$axType)) {
  
  mseType <- mseFrame[mseFrame$axType == i, ]
  
  rmseList[[i]] <- ggplot(data=mseType, aes(x=x, y=mse, group=step)) +
    # geom_line() +
    geom_smooth(colour = "black", method = "loess", span = 1.2, size = 0.5) +
    geom_point(aes(colour = dist)) + 
    geom_point(shape = 1, colour = "black") +
    scale_x_continuous(breaks=seq(5,90,10), labels = seq(70,150,10)) +
    scale_y_continuous(breaks=seq(2000, 16000, 2000), limits = c(2000,16000)) +
    scale_colour_gradientn(name = "Distance",
                           colours = pal,
                           breaks = seq(3,30,3)) +
    # geom_vline(xintercept=seq(0,90,10), colour = "gray", linetype="dashed", ) +
    theme_bw() + 
    theme(panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          legend.key.size = unit(1.9, "cm"),
          legend.key.width = unit(0.5,"cm"),
          plot.title = element_text(hjust = 0.5, size = 22),
          axis.text.x = element_text(size = 13, color = "black"),
          axis.text.y = element_text(size = 13, color = "black"),
          axis.title.x = element_text(size = 16),
          axis.title.y = element_text(size = 16),
          panel.background = element_rect(colour = "black",size = 1)) +
    xlab("Step") +
    ylab("RMSE") +
    ggtitle(i)
  
}

# Save as png
png(filename=paste("rmseEstereo_loess_border.png", sep=""), type="cairo",
    units="in", width=22, height=5, pointsize=12,
    res=300)

gridExtra::grid.arrange(grobs = rmseList, ncol = 4, nrow = 1)

dev.off()
```

###RMSE smooth + examples

```{r, fig.height=5, fig.width=20}
# Color palette
pal <- wes_palette("Zissou1", 10, type = "continuous")
list370 <- list()
list30150 <- list()

# Plot and store in list
for (i in unique(allData$neurType)) {

  #EXAMPLE 3.70
  data2d <- allData[allData$neurType == i & allData$dist == 3 & allData$step == 70, ]
  dataReal <- data.frame(LR = data2d$LR, LRe = data2d$LR)
  
  list370[[i]] <- ggplot(data=data2d, aes(x=LRe, y=LR)) +
    geom_point(colour = "black", size = 0.6) + 
    # geom_point(shape = 1, colour = "black") +
    geom_point(data = dataReal, colour = "red", size = 0.6) + 
    # geom_point(data = dataReal, shape = 1, colour = "black") +
    theme_bw() + 
    theme(panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          plot.title = element_text(hjust = 0.5, size = 12),
          axis.title.x=element_text(size = 10),
          axis.text.x=element_blank(),
          axis.ticks.x=element_blank(),
          axis.title.y=element_text(size = 10),
          axis.text.y=element_blank(),
          axis.ticks.y=element_blank(), 
          panel.background = element_rect(colour = "black",size = 1)) +
    ggtitle("70.3") + 
    xlab("Predicted length") +
    ylab("Real Length")
  
  #EXAMPLE 30.150
  data2d <- allData[allData$neurType == i & allData$dist == 30 & allData$step == 150, ]
  dataReal <- data.frame(LR = data2d$LR, LRe = data2d$LR)
  
  list30150[[i]] <- ggplot(data=data2d, aes(x=LRe, y=LR)) +
    geom_point(colour = "black",  size = 0.6) + 
    # geom_point(shape = 1, colour = "black") +
    geom_point(data = dataReal, colour = "red", size = 0.6) + 
    # geom_point(data = dataReal, shape = 1, colour = "black") +
    theme_bw() + 
    theme(panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          plot.title = element_text(hjust = 0.5, size = 12),
          axis.title.x= element_text(size = 10),
          axis.text.x=element_blank(),
          axis.ticks.x=element_blank(),
          axis.title.y= element_text(size = 10),
          axis.text.y=element_blank(),
          axis.ticks.y=element_blank(),
          panel.background = element_rect(colour = "black",size = 1)) + 
    ggtitle("150.30") + 
    xlab("Real Model") +
    ylab("Real Length")
  
}

# Save as png
for (i in unique(allData$neurType)) {
  
  png(filename=paste("rmseEstereo_loess_border_examples_", i, ".png", sep=""), type="cairo",
      units="in", width=2.5, height=1.5, pointsize=12,
      res=300)
  
  gridExtra::grid.arrange(grobs = c(list370[i], list30150[i]), ncol = 2, nrow = 1)
  
  
  dev.off()
  
}

```

##Mean error

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

# Library
library(latticeExtra)

# Contour color palette
# colfunc <- colorRampPalette(c("navy","royalblue","springgreen","gold","yellow"))
# coul <- colfunc(1000)
# library(viridisLite)
coul <- pals::parula(10000)


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

for (i in c("Type1", "Type2", "Type3", "Type4")) {
  
  heatList[[i]] <- levelplot(error ~ dist * step, data = meanFrame[meanFrame$neurType == i, ], col.regions = coul,
                             scales = list(y = list(at = unique(allData$step), labels = unique(allData$step), cex = 1.5),
                                           x = list(at = unique(allData$dist), labels = unique(allData$dist)), cex = 1.5),
                             ylim = c(155, 65),
                             ylab = list(label = "Step", cex = 1.5),
                             xlab = list(label = "Distance", cex = 1.5),
                             main = list(label = paste(i, "Mean % Error"), cex = 2),
                             colorkey = list(at = (seq(min(meanFrame[meanFrame$neurType == i, "error"]),
                                                       max(meanFrame[meanFrame$neurType == i, "error"]),
                                                       0.05)),
                                             cex = 1.2)) 
  
  smoothList[[i]] <- levelplot(error ~ dist * step, data = meanFrame[meanFrame$neurType == i, ], col.regions = coul,
                               scales = list(y = list(at = unique(allData$step), labels = unique(allData$step), cex = 1.3),
                                             x = list(at = unique(allData$dist), labels = unique(allData$dist)), cex = 1.3),
                               ylim = c(150,70),
                               xlim = c(3,30),
                               ylab = list(label = "Step", cex = 1.5),
                               xlab = list(label = "Distance", cex = 1.5),
                               main = list(label = paste(i, "Smooth"), cex = 1.5),
                               cuts = 20, 
                               colorkey = list(at = (seq(min(meanFrame[meanFrame$neurType == i, "error"]),
                                                         max(meanFrame[meanFrame$neurType == i, "error"]),
                                                         0.05)),
                                               cex = 1.2),
                               panel = panel.levelplot.points, cex = 0) + 
    layer_(panel.2dsmoother(..., n = 200))
  
}


# Plot

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

gridExtra::grid.arrange(grobs = c(heatList, smoothList), ncol = 4, nrow = 2)

dev.off()

```

##Mean error (Type123 vs 4)

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

# Library
library(latticeExtra)

# Contour color palette
# colfunc <- colorRampPalette(c("navy","royalblue","springgreen","gold","yellow"))
# coul <- colfunc(1000)
# library(viridisLite)
coul <- pals::parula(10000)


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

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

for (i in c("Type1", "Type2", "Type3", "Type4")) {
  
  if (i != "Type4") {
    
  heatList[[i]] <- levelplot(error ~ dist * step, data = meanFrame[meanFrame$neurType == i, ], col.regions = coul,
                             scales = list(y = list(at = unique(allData$step), labels = unique(allData$step), cex = 1.5),
                                           x = list(at = unique(allData$dist), labels = unique(allData$dist)), cex = 1.5),
                             ylim = c(155, 65),
                             ylab = list(label = "Step", cex = 1.5),
                             xlab = list(label = "Distance", cex = 1.5),
                             main = list(label = paste(i, "Mean % Error"), cex = 2),
                             colorkey = list(at = (seq(min123,
                                                       max123,
                                                       0.05)),
                                             cex = 1.2)) 
  
  smoothList[[i]] <- levelplot(error ~ dist * step, data = meanFrame[meanFrame$neurType == i, ], col.regions = coul,
                               scales = list(y = list(at = unique(allData$step), labels = unique(allData$step), cex = 1.3),
                                             x = list(at = unique(allData$dist), labels = unique(allData$dist)), cex = 1.3),
                               ylim = c(150,70),
                               xlim = c(3,30),
                               ylab = list(label = "Step", cex = 1.5),
                               xlab = list(label = "Distance", cex = 1.5),
                               main = list(label = paste(i, "Smooth"), cex = 1.5),
                               cuts = 20, 
                               colorkey = list(at = (seq(min123,
                                                         max123,
                                                         0.05)),
                                               cex = 1.2),
                               panel = panel.levelplot.points, cex = 0) + 
    layer_(panel.2dsmoother(..., n = 200))
  
  } else {
    
     heatList[[i]] <- levelplot(error ~ dist * step, data = meanFrame[meanFrame$neurType == i, ], col.regions = coul,
                             scales = list(y = list(at = unique(allData$step), labels = unique(allData$step), cex = 1.5),
                                           x = list(at = unique(allData$dist), labels = unique(allData$dist)), cex = 1.5),
                             ylim = c(155, 65),
                             ylab = list(label = "Step", cex = 1.5),
                             xlab = list(label = "Distance", cex = 1.5),
                             main = list(label = paste(i, "Mean % Error"), cex = 2),
                             colorkey = list(at = (seq(min(meanFrame[meanFrame$neurType == i, "error"]),
                                                       max(meanFrame[meanFrame$neurType == i, "error"]),
                                                       0.05)),
                                             cex = 1.2)) 
  
  smoothList[[i]] <- levelplot(error ~ dist * step, data = meanFrame[meanFrame$neurType == i, ], col.regions = coul,
                               scales = list(y = list(at = unique(allData$step), labels = unique(allData$step), cex = 1.3),
                                             x = list(at = unique(allData$dist), labels = unique(allData$dist)), cex = 1.3),
                               ylim = c(150,70),
                               xlim = c(3,30),
                               ylab = list(label = "Step", cex = 1.5),
                               xlab = list(label = "Distance", cex = 1.5),
                               main = list(label = paste(i, "Smooth"), cex = 1.5),
                               cuts = 20, 
                               colorkey = list(at = (seq(min(meanFrame[meanFrame$neurType == i, "error"]),
                                                         max(meanFrame[meanFrame$neurType == i, "error"]),
                                                         0.05)),
                                               cex = 1.2),
                               panel = panel.levelplot.points, cex = 0) + 
    layer_(panel.2dsmoother(..., n = 200))
    
  }

  
}


# Plot

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

gridExtra::grid.arrange(grobs = c(heatList, smoothList), ncol = 4, nrow = 2)

dev.off()

```


##Errors vs cumul probability

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

# Function to apply to all axon types
quantType <- function(peFrame, probSeq) {
  
  probList <- list()
  
  for (i in unique(peFrame$interact)) {
    
    dataProb <- peFrame[peFrame$interact == i, "error"]
    probVec <- numeric()
    
    for (j in probSeq) {
      
      errProb <- quantInv(dataProb, j) - quantInv(dataProb, -j)
      probVec <- c(probVec, errProb)
      
    }
    
    probList[[i]] <- probVec
    
  }
  
  return(probList)
  
}

# Define errors for which to calculate probability
binProb <- 0.5
probSeq <- seq(binProb, 100, binProb)
# Store axon types in lists
frameList <- list(Type1 = allData[allData$neurType == "Type1", ],
                        Type2 = allData[allData$neurType == "Type2", ],
                        Type3 = allData[allData$neurType == "Type3", ],
                        Type4 = allData[allData$neurType == "Type4", ])
axProb <- lapply(frameList, function(x) quantType(x, probSeq))
saveRDS(axProb, "errorProbs_RMSE.rds")
```

######Plot heatmap 5, 10

```{r, fig.height=10, fig.width=15}
library(latticeExtra)

# Reformat as dataframe
binProb <- 0.5
probSeq <- seq(binProb, 100, binProb)
axProb <- readRDS("errorProbs_RMSE.rds")
probFrame <- data.frame(error = rep(probSeq, 90*4), 
                        neurType = rep(c("Type1", "Type2", "Type3", "Type4"), each = length(probSeq)*90),
                        dist = rep(rep(unique(allData$dist), each = length(probSeq)*9), 4),
                        step = rep(rep(unique(allData$step), each = length(probSeq)), 10*4),
                        prob = unlist(axProb))

# Plot contour plot for 5% error
coul <- viridis::magma(10000)

# Store heatmaps
heatProbList <- list()
smoothProbList <- list()

for (i in c("Type1", "Type2", "Type3", "Type4")) {
  
  levelList1 <- list()
  levelList2 <- list()
  
  for (j in c(5, 10, 20)) {
    
    dataPlot <- probFrame[probFrame$neurType == i & probFrame$error == j, ]
    
    levelList1[[as.character(j)]] <- levelplot(prob ~ dist * step, data = dataPlot, 
                                               col.regions = coul,
                                               scales = list(y = list(at = unique(allData$step), labels = unique(allData$step), cex = 1.3),
                                                             x = list(at = unique(allData$dist), labels = unique(allData$dist)), cex = 1.3),
                                               ylim = c(155,65),
                                               colorkey = list(at = (seq(min(dataPlot$prob),
                                                                         max(dataPlot$prob),
                                                                         0.0025))),
                                               ylab = list(label = "Step", cex = 1.5),
                                               xlab = list(label = "Distance", cex = 1.5),
                                               main = list(label = paste(i, " P(Error <= ", j, " %)", sep = ""), cex = 2)) 
    
    levelList2[[as.character(j)]] <- levelplot(prob ~ dist * step, data = dataPlot, 
                                               col.regions = coul,
                                               scales = list(y = list(at = unique(allData$step), labels = unique(allData$step), cex = 1.3),
                                                             x = list(at = unique(allData$dist), labels = unique(allData$dist)), cex = 1.3),
                                               ylim = c(150,70),
                                               xlim = c(3,30),
                                               cuts = 20,
                                               colorkey = list(at = (seq(min(dataPlot$prob),
                                                                         max(dataPlot$prob),
                                                                         0.0025))),
                                               ylab = list(label = "Step", cex = 1.5),
                                               xlab = list(label = "Distance", cex = 1.5),
                                               main = list(label = paste(i, " P(Error <= ", j, " %)", sep = ""), cex = 2),
                                               panel = panel.levelplot.points, cex = 0) + 
      layer_(panel.2dsmoother(..., n = 200))
  }
  
  heatProbList[[i]] <- levelList1
  smoothProbList[[i]] <- levelList2
}

smoothProb5 <- lapply(smoothProbList, function(x) x[[1]])
smoothProb10 <- lapply(smoothProbList, function(x) x[[2]])
smoothProb20 <- lapply(smoothProbList, function(x) x[[3]])

# Plot
# setwd("EstereoAnalysis/")
for (i in c("Type1", "Type2", "Type3", "Type4")) {

  png(filename=paste("errorProbabilities_heatmap_", i, "_MSE.png", sep=""), type="cairo",
      units="in", width=15, height=10, pointsize=12,
      res=300)

  gridExtra::grid.arrange(grobs = c(heatProbList[[i]], smoothProbList[[i]]), ncol = 3, nrow = 2)
  
  dev.off()

}


```


##Common plot

```{r, fig.height=40, fig.width=40}
gridExtra::grid.arrange(grobs = c(rmseList, smoothList, smoothProb5, smoothProb10), ncol = 4, nrow = 4)
```

##Boxplot error~dist:step

```{r, fig.width=20, fig.height=5}
# Sort by step
dataStep <- allData[order(allData$step), ]
dataStep$x <- rep(1:10, 8559)

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

# ggplot(dataStep, aes(x=factor(step), y=error, fill=factor(dist))) +
#   geom_boxplot() +
#   scale_fill_manual(values=pal) + 
#   ylim(c(-5,100))

# Plot and store in list
for (i in unique(allData$neurType)) {
  
  errType <- dataStep[dataStep$neurType == i, ]
  
  meanErr <- aggregate(error ~ step + dist, errType, mean)
  meanErr <- meanErr[order(meanErr$step), ]
  
  errorList[[i]] <- ggplot(errType, aes(x=factor(step), y=error)) +
    geom_boxplot(aes(fill=interact), outlier.shape = NA) +
    scale_fill_manual(values=pal) + 
    # geom_line() +
    # geom_boxplot() + 
    # scale_x_continuous(breaks=seq(5,90,10), labels = seq(70,150,10)) +
    # scale_y_continuous(breaks=seq(2000, 16000, 2000), limits = c(2000,16000)) +
    # scale_colour_gradientn(name = "Distance",
    #                        colours = pal,
    #                        breaks = seq(3,30,3)) +
    # geom_vline(xintercept=seq(0,90,10), colour = "gray", linetype="dashed", ) +
    geom_smooth(method="loess", se=TRUE) +
    theme_bw() + 
    theme(panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          legend.key.size = unit(1, "cm"),
          legend.key.width = unit(0.5,"cm"),
          plot.title = element_text(hjust = 0.5, size = 22),
          axis.text.x = element_text(size = 13, color = "black"),
          axis.text.y = element_text(size = 13, color = "black"),
          axis.title.x = element_text(size = 16),
          axis.title.y = element_text(size = 16)) +
    xlab("Step") +
    ylab("RMSE") +
    ggtitle(i) +
    scale_y_continuous(limits = c(0,200))
  
}

# Save as png
# png(filename=paste("rmseEstereo2.png", sep=""), type="cairo",
#     units="in", width=22, height=5, pointsize=12,
#     res=300)

gridExtra::grid.arrange(grobs = errorList, ncol = 4, nrow = 1)

# dev.off()
```


```{r}
data=data.frame(date=as.Date(c("2011-02-10","2011-02-10","2011-02-10","2011-02-10","2011-02-10",
                         "2011-02-10","2011-02-10","2011-02-10","2011-02-10","2011-02-10",
                         "2011-02-20","2011-02-20","2011-02-20","2011-02-20","2011-02-20",
                         "2011-02-20","2011-02-20","2011-02-20","2011-02-20","2011-02-20",
                         "2011-02-28","2011-02-28","2011-02-28","2011-02-28","2011-02-28",
                         "2011-02-28","2011-02-28","2011-02-28","2011-02-28","2011-02-28",
                         "2011-03-10","2011-03-10","2011-03-10","2011-03-10","2011-03-10",
                         "2011-03-10","2011-03-10","2011-03-10","2011-03-10","2011-03-10"),format="%Y-%m-%d"),
            spore=c(0,1,0,1,0,
                    1,2,0,1,1,
                    8,5,6,12,7,
                    7,5,4,7,6,
                    18,24,25,32,14,
                    24,16,18,23,30,
                    27,26,36,31,22,
                    28,29,37,30,24),
            house = rep(c("Int", "Ext"), each = 5, times = 4)
            )
data$interact <- interaction(data$date, data$house)


ggplot(data,aes(x=date,y=spore)) + 
  geom_boxplot(aes(fill=interact), outlier.shape = NA) + 
  # scale_fill_manual(values=rep(pal,9)) +
  geom_smooth() + 
  # scale_y_continuous(limits = c(0,200)) +
  theme(legend.position = "none")


ggplot(allData,aes(x=step,y=error)) + 
  geom_boxplot(aes(fill=interact), outlier.shape = NA) + 
  scale_fill_manual(values=rep(pal,9)) +
  geom_smooth(method="loess") + 
  scale_y_continuous(limits = c(0,200)) +
  theme(legend.position = "none")
```