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