123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757758759760761762763764765766767768769770771772773774775776777778779780781782783784785786787788789790791792793794795796797798799800801802803804805806807808809810811812813814815816817818819820821822823824825826827828829830831832833834835836837838839840841842843 |
- ---
- 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(ggridges)
- library(hrbrthemes)
- library(viridis)
- library(ggpubr)
- library(mltools)
- library(data.table)
- library(caret)
- library(interactions)
- # library(performance)
- library(MASS)
- library(geepack)
- library(pstools)
- library(MXM)
- library(rmcorr)
- library(multcomp)
- library(Hmisc)
- options(digits = 10)
- # https://rstudio-pubs-static.s3.amazonaws.com/297778_5fce298898d64c81a4127cf811a9d486.html
- ################
- # Bootstrap PI #
- ################
- the.replication <- function(reg, s, x_Np1 = 0){
- # Make bootstrap residuals
- ep.star <- sample(s, size=length(reg$residuals),replace=TRUE)
- # Make bootstrap Y
- y.star <- fitted(reg) + ep.star
- # Do bootstrap regression
- lp <- model.frame(reg)[,2]
- nt <- model.frame(reg)[,3]
- bs.reg <- lm(y.star ~ lp:nt - 1)
- # Create bootstrapped adjusted residuals
- bs.lev <- influence(bs.reg)$hat
- bs.s <- residuals(bs.reg)/sqrt(1-bs.lev)
- bs.s <- bs.s - mean(bs.s)
- # Calculate draw on prediction error
- # xb.xb <- coef(my.reg)["(Intercept)"] - coef(bs.reg)["(Intercept)"]
- xb.xb <- (coef(my.reg)[1] - coef(bs.reg)[1])*x_Np1
- return(unname(xb.xb + sample(bs.s, size=1)))
-
- }
- ############################
- # Bootstrap PI generalized #
- ############################
- the.replication.gen <- function(reg, s, axType, x_Np1 = 0){
- # Make bootstrap residuals
- ep.star <- sample(s, size=length(reg$residuals), replace=TRUE)
- # Make bootstrap Y
- y.star <- fitted(reg) + ep.star
- # Do bootstrap regression
- lp <- model.frame(reg)[,2]
- nt <- model.frame(reg)[,3]
- bs.reg <- lm(y.star ~ lp:nt - 1)
- # Create bootstrapped adjusted residuals
- bs.lev <- influence(bs.reg)$hat
- bs.s <- residuals(bs.reg)/sqrt(1-bs.lev)
- bs.s <- bs.s - mean(bs.s)
- # Calculate draw on prediction error
- # xb.xb <- coef(my.reg)["(Intercept)"] - coef(bs.reg)["(Intercept)"]
- xb.xb <- (coef(my.reg)[axType] - coef(bs.reg)[axType])*x_Np1
- return(unname(xb.xb + sample(bs.s, size=1)))
-
- }
- ##############################
- # Bootstrap PI per axon type #
- ##############################
- the.replication.type <- 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)
- # Create bootstrapped adjusted residuals
- bs.lev <- influence(bs.reg)$hat
- bs.s <- residuals(bs.reg)/sqrt(1-bs.lev)
- bs.s <- bs.s - mean(bs.s)
- # Calculate draw on prediction error
- # xb.xb <- coef(my.reg)["(Intercept)"] - coef(bs.reg)["(Intercept)"]
- xb.xb <- (coef(my.reg) - coef(bs.reg))*x_Np1
- return(unname(xb.xb + sample(bs.s, size=1)))
-
- }
- ##############################
- # Stratified random sampling #
- ##############################
- stratified <- function(df, group, size, select = NULL,
- replace = FALSE, bothSets = FALSE) {
- if (is.null(select)) {
- df <- df
- } else {
- if (is.null(names(select))) stop("'select' must be a named list")
- if (!all(names(select) %in% names(df)))
- stop("Please verify your 'select' argument")
- temp <- sapply(names(select),
- function(x) df[[x]] %in% select[[x]])
- df <- df[rowSums(temp) == length(select), ]
- }
- df.interaction <- interaction(df[group], drop = TRUE)
- df.table <- table(df.interaction)
- df.split <- split(df, df.interaction)
- if (length(size) > 1) {
- if (length(size) != length(df.split))
- stop("Number of groups is ", length(df.split),
- " but number of sizes supplied is ", length(size))
- if (is.null(names(size))) {
- n <- setNames(size, names(df.split))
- message(sQuote("size"), " vector entered as:\n\nsize = structure(c(",
- paste(n, collapse = ", "), "),\n.Names = c(",
- paste(shQuote(names(n)), collapse = ", "), ")) \n\n")
- } else {
- ifelse(all(names(size) %in% names(df.split)),
- n <- size[names(df.split)],
- stop("Named vector supplied with names ",
- paste(names(size), collapse = ", "),
- "\n but the names for the group levels are ",
- paste(names(df.split), collapse = ", ")))
- }
- } else if (size < 1) {
- n <- round(df.table * size, digits = 0)
- } else if (size >= 1) {
- if (all(df.table >= size) || isTRUE(replace)) {
- n <- setNames(rep(size, length.out = length(df.split)),
- names(df.split))
- } else {
- message(
- "Some groups\n---",
- paste(names(df.table[df.table < size]), collapse = ", "),
- "---\ncontain fewer observations",
- " than desired number of samples.\n",
- "All observations have been returned from those groups.")
- n <- c(sapply(df.table[df.table >= size], function(x) x = size),
- df.table[df.table < size])
- }
- }
- temp <- lapply(
- names(df.split),
- function(x) df.split[[x]][sample(df.table[x],
- n[x], replace = replace), ])
- set1 <- do.call("rbind", temp)
-
- if (isTRUE(bothSets)) {
- set2 <- df[!rownames(df) %in% rownames(set1), ]
- list(SET1 = set1, SET2 = set2)
- } else {
- set1
- }
- }
- ```
- ##Load data
- ```{r}
- # Get file paths
- axonNames <- list.files(paste("ProyeccionData/", sep=""), full.names=T)
- # Load matlab file
- axonFiles <- lapply(axonNames, function(x) readMat(x))
- names(axonFiles) <- c("Type1", "Type2", "Type3", "Type4")
- # VAR NAMES
- # REAL_LENGTH, PROYECTION_LENGTH_XY, PROYECTION_LENGTH_XZ, PROYECTION_LENGTH_YZ
- # Extract data
- realLength <- lapply(axonFiles, function(x) x$REAL.LENGTH)
- planeXY <- lapply(axonFiles, function(x) x$PROYECTION.LENGTH.XY)
- planeXZ <- lapply(axonFiles, function(x) x$PROYECTION.LENGTH.XZ)
- planeYZ <- lapply(axonFiles, function(x) x$PROYECTION.LENGTH.YZ)
- # Get number of neurons per type
- numberTypes <- unlist(lapply(realLength, function(x) length(x)))
- # Store in data frames by plane
- xyData <- data.frame(id = 1:length(unlist(realLength)),
- LR = unlist(realLength), LP = unlist(planeXY), Plane = rep("XY", sum(numberTypes)),
- neurType = rep(c("Type1", "Type2", "Type3", "Type4"), numberTypes))
- xzData <- data.frame(id = 1:length(unlist(realLength)),
- LR = unlist(realLength), LP = unlist(planeXZ), Plane = rep("XZ", sum(numberTypes)),
- neurType = rep(c("Type1", "Type2", "Type3", "Type4"), numberTypes))
- yzData <- data.frame(id = 1:length(unlist(realLength)),
- LR = unlist(realLength), LP = unlist(planeYZ), Plane = rep("YZ", sum(numberTypes)),
- neurType = rep(c("Type1", "Type2", "Type3", "Type4"), numberTypes))
- # Merge into a single data frame and sort by id
- allData <- do.call("rbind", list(xyData, xzData, yzData))
- allData <- allData[order(allData$id),]
- # Add different number for every neuron-plane combination
- allData$NeurPlane <- c(rep(c(1,2,3), numberTypes[1]), rep(c(4,5,6), numberTypes[2]),
- rep(c(7,8,9), numberTypes[3]), rep(c(10,11,12), numberTypes[4]))
- allData$interact <-interaction(allData$neurType, allData$Plane, lex.order=T)
- allData$idR <- 1:dim(allData)[1]
- # Create data frame w/o Type4
- dataNo4 <- allData[allData$neurType != "Type4", ]; dataNo4$neurType <- factor(dataNo4$neurType)
- # Data in short format
- dataShort <- data.frame(id = 1:length(unlist(realLength)),
- LR = unlist(realLength), XY = unlist(planeXY), XZ = unlist(planeXZ), YZ = unlist(planeYZ),
- neurType = rep(c("Type1", "Type2", "Type3", "Type4"), numberTypes))
- ```
- ##CV-Bootstrap Pred Interval
- RUN ONLY ONCE, AFTERWARDS LOAD THE FILE IN MEAN ERRORS CHUNK
- ```{r}
- # Ensure reproducibility
- set.seed(123456)
- axList <- list()
- for (h in unique(allData$neurType)) {
-
- # Subset axon type
- dataCV <- allData[allData$neurType == h, ]
-
- # List to store CV repetitions
- cvList <- list ()
- cvRand <- list()
-
- # Stratified random 5-fold sets
- setList <- list()
- sampSize <- 0.2*length(unique(dataCV$id))
-
- groupDat1 <- dataCV %>% group_by(Plane)
- set1<- sample_n(groupDat1, sampSize)
-
- groupDat2 <- groupDat1[-set1$idR, ]
- set2 <- sample_n(groupDat2, sampSize)
-
- groupDat3 <- groupDat2[-set2$idR, ]
- set3 <- sample_n(groupDat3, sampSize)
-
- groupDat4 <- groupDat3[-set3$idR, ]
- set4 <- sample_n(groupDat4, sampSize)
-
- groupDat5 <- groupDat4[-set4$idR, ]
- set5 <- sample_n(groupDat5, sampSize)
-
- # Store sets (ungrouped)
- setList <- lapply(list(set1, set2, set3, set4, set5), ungroup)
- axList[[h]]$sets <- setList
-
- # List to store CV
- cvList <- list()
-
- for (i in 1:5) {
-
- # Stratified random training/test sets
- trainSet <- do.call(rbind, setList[-i])
- testSet <- setList[[i]]
-
- # 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)
-
- # List to store planes
- planeList <- list()
-
- # Bootstrapping on test set
- for (j in c("XY", "XZ", "YZ")) {
-
- # Subset plane
- dataSub <- testSet[testSet$Plane == j, ]
-
- # List to store errors
- epList <- list()
- epCount <- 1
-
- for (k in dataSub[ , "LP"]) {
-
- # Do bootstrap with 100 replications
- ep.draws <- replicate(n=100, the.replication.type(reg = my.reg, s = my.s.resid, x_Np1 = k))
-
- # Store in list
- epList[[epCount]] <- ep.draws
- epCount <- epCount + 1
-
- }
-
- # Store in list
- planeList[[j]] <- epList
-
- print(paste("CV", i, "Plane", j, "finished"))
-
- }
-
- # Store in list
- cvList[[i]] <- planeList
-
- }
-
- # Store in list
- axList[[h]]$cv <- cvList
-
- print(paste("Axon", h, "finished"))
-
- }
- saveRDS(axList, paste("projection_5CV_allAxons.rds", sep = ""))
- ```
- ##Estimate porcentual errors
- ```{r, fig.width=8, fig.height=8}
- # Load 5-fold CV
- axList <- readRDS("ProyeccionAnalysis/projection_5CV_allAxons.rds")
- # Estimate and store in list
- peList <- list()
- peCount <- 1
- lrLength <- list()
- for (h in c("Type1", "Type2", "Type3", "Type4")) {
-
- for (i in 1:5) {
-
- for (j in c("XY", "XZ", "YZ")) {
-
- # Extract LR values
- lr <- axList[[h]]$sets[[i]]
- lr <- lr[lr$Plane == j, "LR"]
-
- # Extract LP values
- lp <- axList[[h]]$cv[[i]][[j]][[1]]
-
- for (k in 1:dim(lr)[1]) {
-
- # Divide errors by LR to get porcentual errors
- pe <- lp[k, ] / lr$LR[k]
-
- # Store in list
- peList[[peCount]] <- pe*100
- peCount <-peCount + 1
- }
-
- }
-
- }
-
- lrLength[h] <- dim(lr)[1]
-
- }
- # Store in data frame
- lrLength <- unlist(lrLength)
- peFrame <- data.frame(pe = unlist(peList),
- neurType = rep(c("Type1", "Type2", "Type3", "Type4"), times = lrLength*5*3*100),
- Plane = c(rep(rep(c("XY", "XZ", "YZ"), each = lrLength[1]*100), 5),
- rep(rep(c("XY", "XZ", "YZ"), each = lrLength[2]*100), 5),
- rep(rep(c("XY", "XZ", "YZ"), each = lrLength[3]*100), 5),
- rep(rep(c("XY", "XZ", "YZ"), each = lrLength[4]*100), 5)))
- peFrame$interact <- interaction(peFrame$neurType, peFrame$Plane, lex.order=T)
- ```
- ####Plot density and histogram
- ```{r}
- # Plot
- # peGroup <- peFrame %>% group_by(interact)
- # peSample <- sample_n(peGroup, 100)
- setwd("ProyeccionFigures/")
- png(filename=paste("prederrorDistributions_5CV.png", sep=""), type="cairo",
- units="in", width=8, height=10, pointsize=12,
- res=300)
- (p <- ggplot(peFrame, 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))
- dev.off()
- png(filename=paste("prederrorHistogram_5CV.png", sep=""), type="cairo",
- units="in", width=8, height=10, pointsize=12,
- res=300)
- # Histogram
- # nint 300 bootPI
- # nint 500 bootPI_full
- # nint 400 bootPI_full_replace
- (h <- histogram( ~ pe | Plane + neurType, peFrame, nint = 300, type ="density", xlim = c(-50,50)))
- dev.off()
- ```
- #####Insets
- ```{r}
- # Find mean bandwidth
- bwFrame <- aggregate(pe ~ neurType + Plane, peFrame, function(x) density(x)$bw)
- # Named color vector
- colVec <- c("limegreen", "blue", "red")
- names(colVec) <- c("XY", "XZ", "YZ")
- setwd("ProyeccionFigures/")
- for (i in unique(peFrame$neurType)) {
-
- png(filename=paste("insetProb_", i, ".png", sep=""), type="cairo",
- units="in", width=5, height=5, pointsize=12,
- res=300)
- par(mfrow = c(3,2),
- oma = c(5,4,0,0) + 0.1,
- mar = c(0,0,1,1) + 0.1)
-
- for (j in unique(peFrame$Plane)) {
-
- # Estimate density
- dens <- density(peFrame[peFrame$neurType == i & peFrame$Plane == j, "pe"], bw = mean(bwFrame$pe))
- # 5%
- x1 <- min(which(dens$x >= -5))
- x2 <- max(which(dens$x < 5))
-
- plot(dens, ylim = c(0, 0.15), xlim = c(-20,20), zero.line = FALSE, main = "", ylab = "", xlab = "", axes = FALSE)
- with(dens, polygon(x=c(x[c(x1,x1:x2,x2)]), y= c(0, y[x1:x2], 0), col=colVec[[j]], border = NA))
- lines(dens)
-
- # Add y axis
- axis(side = 2, at=seq(-0.05,0.15,0.05), labels = seq(-0.05,0.15,0.05), cex.axis = 1.2)
- # Different x axis for bottom plot
- if (j == "YZ") {
- axis(side = 1, at=seq(-30,30,10), labels = seq(-30,30,10), cex.axis = 1.2)
- } else {
- axis(side = 1, at=seq(-30,30,10), labels = FALSE)
- }
-
-
- # 10%
- x1 <- min(which(dens$x >= -10))
- x2 <- max(which(dens$x < 10))
-
- plot(dens, ylim = c(0, 0.15), xlim = c(-20,20), zero.line = FALSE, main = "", ylab = "", xlab = "", axes = FALSE)
- with(dens, polygon(x=c(x[c(x1,x1:x2,x2)]), y= c(0, y[x1:x2], 0), col=colVec[[j]], border = NA))
- lines(dens)
-
- # Add y axis
- axis(side = 2, at=seq(-0.05,0.15,0.05), labels = FALSE)
- # Different x axis for bottom plot
- if (j == "YZ") {
- axis(side = 1, at=seq(-30,30,10), labels = seq(-30,30,10), cex.axis = 1.2)
- } else {
- axis(side = 1, at=seq(-30,30,10), labels = FALSE)
- }
-
- }
-
-
- title(xlab = "% Error",
- ylab = "Density",
- outer = TRUE, line = 3,
- cex.lab = 1.5)
-
- dev.off()
-
- }
- ```
- ####Errors vs cumul probability
- RUN ONLY ONCE, AFTERWARDS LOAD THE FILE IN PLOT CHUNKS
- ```{r, fig.height=5, fig.width=20}
- # Inverse quantile
- quantInv <- function(distr, value) ecdf(distr)(value)
- #ecdf plot example
- # plot(ecdf(peFrame[peFrame$interact == "Type1.XY", "pe"]))
- # lines(ecdf(peFrame[peFrame$interact == "Type1.XZ", "pe"]), col = "red")
- # lines(ecdf(peFrame[peFrame$interact == "Type1.YZ", "pe"]), col = "blue")
- probList <- list()
- binProb <- 0.5
- probSeq <- seq(binProb, 100, binProb)
- 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
-
- }
- setwd("ProyeccionAnalysis/")
- saveRDS(probList, "cumProb_5CV.rds")
- ```
- #####Plot
- ```{r, fig.height=5, fig.width=20}
- # Plot with limited axis
- setwd("ProyeccionFigures/")
- probList <- readRDS("cumProb_5CV.rds")
- png(filename=paste("errorProbabilityCum_bin", binProb, "_5CV.png", sep=""), type="cairo",
- units="in", width=20, height=5, pointsize=12,
- res=300)
- par(mfrow = c(1,4))
- hlist <- list()
- axCount <- 1
- for (i in seq(1,10,3)) {
- # Extract height for 5% error
- h5 <- sapply(probList[i:(i+2)], function(x) x[5/binProb])
- max5 <- max(h5)
-
- # Extract height for 10% error
- h10 <- sapply(probList[i:(i+2)], function(x) x[10/binProb])
- max10 <- max(h10)
-
- # Plot cumulative curves
- plot(probSeq, probList[[i]], type = "l", col = "limegreen", lwd = 1.5, ylim = c(0,1), xlim = c(0.5, 40),
- ylab = "Cumulative probability", xlab = "% Error", main = paste("Axon Type", axCount))
- lines(probSeq, probList[[i + 1]], col = "red", lwd = 1.5)
- lines(probSeq, probList[[i + 2]], col = "blue", lwd = 1.5)
-
- # Plot vertical segments at 5% and 10%
- segments(5, -0.1, 5, max5, lty = "dashed", col = "gray")
- segments(10, -0.1, 10, max10, lty = "dashed", col = "gray")
-
- # Plot horizontal segments at cumul prob for 5% and 10%
- transp <- 0.8
- segments(-2, h5[1], 5, h5[1], lty = "dashed", col = alpha("limegreen", transp))
- segments(-2, h5[2], 5, h5[2], lty = "dashed", col = alpha("red", transp))
- segments(-2, h5[3], 5, h5[3], lty = "dashed", col = alpha("blue", transp))
-
- segments(-2, h10[1], 10, h10[1], lty = "dashed", col = alpha("limegreen", transp))
- segments(-2, h10[2], 10, h10[2], lty = "dashed", col = alpha("red", transp))
- segments(-2, h10[3], 10, h10[3], lty = "dashed", col = alpha("blue", transp))
-
- # Add legend
- # legend("bottomright", c("XY", "XZ", "YZ"), lwd = 1.5, col = c("limegreen", "red", "blue"))
-
- # Store to check later
- hlist[[paste("Type", axCount, sep="")]] <- data.frame(h5 = h5, h10 = h10)
- axCount <- axCount + 1
- }
- dev.off()
- # Plot with free axis
- png(filename=paste("errorProbabilityCum_bin", binProb, "_freeAxis_5CV.png", sep=""), type="cairo",
- units="in", width=20, height=5, pointsize=12,
- res=300)
- par(mfrow = c(1,4))
- hlist <- list()
- axCount <- 1
- for (i in seq(1,10,3)) {
- # Extract height for 5% error
- h5 <- sapply(probList[i:(i+2)], function(x) x[5/binProb])
- max5 <- max(h5)
-
- # Extract height for 10% error
- h10 <- sapply(probList[i:(i+2)], function(x) x[10/binProb])
- max10 <- max(h10)
-
- # Plot cumulative curves
- plot(probSeq, probList[[i]], type = "l", col = "green", lwd = 1.5,
- ylab = "Cumulative probability", xlab = "% Error", main = paste("Axon Type", axCount))
- lines(probSeq, probList[[i + 1]], col = "red", lwd = 1.5)
- lines(probSeq, probList[[i + 2]], col = "blue", lwd = 1.5)
-
- # Plot vertical segments at 5% and 10%
- segments(5, -0.1, 5, max5, lty = "dashed", col = "gray")
- segments(10, -0.1, 10, max10, lty = "dashed", col = "gray")
-
- # Plot horizontal segments at cumul prob for 5% and 10%
- segments(-2, h5[1], 5, h5[1], lty = "dashed", col = "yellowgreen")
- segments(-2, h5[2], 5, h5[2], lty = "dashed", col = "tomato")
- segments(-2, h5[3], 5, h5[3], lty = "dashed", col = "lightskyblue")
-
- segments(-2, h10[1], 10, h10[1], lty = "dashed", col = "yellowgreen")
- segments(-2, h10[2], 10, h10[2], lty = "dashed", col = "tomato")
- segments(-2, h10[3], 10, h10[3], lty = "dashed", col = "lightskyblue")
-
- # Add legend
- legend("bottomright", c("XY", "XZ", "YZ"), lwd = 1.5, col = c("green", "red", "blue"))
-
- # Store to check later
- hlist[[paste("Type", axCount, sep="")]] <- data.frame(h5 = h5, h10 = h10)
- axCount <- axCount + 1
- }
- dev.off()
- ```
- ##Full figure
- ```{r, fig.width=20}
- ## show the regions that have been allocated to each plot
- # layout.show(8)
- #############
- # Load data #
- #############
- probList <- readRDS("ProyeccionAnalysis/cumProb_5CV.rds")
- ##############
- # MAIN PLOTS #
- ##############
- setwd("ProyeccionFigures/")
- png(filename=paste("errorProbabilityCum_bin", binProb, "_5CV_FULL.png", sep=""), type="cairo",
- units="in", width=20, height=5, pointsize=12,
- res=300)
- # Graphical parameters
- par(oma = c(1,4,1,1),
- mar = c(4,0,3,0) + 0.5,
- cex.axis = 1.5,
- cex.lab = 1.5,
- cex.main = 2)
- layout(matrix(c(1,1,2,2,3,3,0,4,4,
- 1,5,2,7,3,9,0,4,11,
- 1,6,2,8,3,10,0,4,12,
- 1,0,2,0,3,0,0,4,0), 4, 9, byrow = TRUE),
- widths = c(2.5,2.5,2.5,2.5,2.5,2.5,1,3,2.5),
- heights = c(1.6,1.6,1.6,1.2))
- # layout.show(12)
- # Plot loop
- hlist <- list()
- axCount <- 1
- for (i in seq(1,10,3)) {
-
- # Extract height for 5% error
- h5 <- sapply(probList[i:(i+2)], function(x) x[5/binProb])
-
- # Extract height for 10% error
- h10 <- sapply(probList[i:(i+2)], function(x) x[10/binProb])
-
- # Plot cumulative error and Axis labels
-
- if (axCount == 1) {
-
- plot(probSeq, probList[[i]], type = "l", col = "limegreen", lwd = 1.5, ylim = c(0,1), xlim = c(0.5, 40),
- main = paste("Axon Type", axCount), yaxt = "n", xaxt = "n", xlab = "", ylab = "")
- lines(probSeq, probList[[i + 1]], col = "red", lwd = 1.5)
- lines(probSeq, probList[[i + 2]], col = "blue", lwd = 1.5)
- title(ylab = "Cumulative probability",
- outer = TRUE, line=3)
-
-
- } else if (axCount == 2) {
-
- plot(probSeq, probList[[i]], type = "l", col = "limegreen", lwd = 1.5, ylim = c(0,1), xlim = c(0.5, 40),
- main = paste("Axon Type", axCount), yaxt = "n", xaxt = "n", xlab = "", ylab = "")
- lines(probSeq, probList[[i + 1]], col = "red", lwd = 1.5)
- lines(probSeq, probList[[i + 2]], col = "blue", lwd = 1.5)
- title(xlab="% Error", line=3)
-
-
- } else if (axCount == 3) {
-
- plot(probSeq, probList[[i]], type = "l", col = "limegreen", lwd = 1.5, ylim = c(0,1), xlim = c(0.5, 40),
- main = paste("Axon Type", axCount), yaxt = "n", xaxt = "n", xlab = "", ylab = "")
- lines(probSeq, probList[[i + 1]], col = "red", lwd = 1.5)
- lines(probSeq, probList[[i + 2]], col = "blue", lwd = 1.5)
-
- } else if (axCount == 4) {
-
- par(mar = c(4,5,3,0) + 0.5)
- plot(probSeq, probList[[i]], type = "l", col = "limegreen", lwd = 1.5, ylim = c(0,1), xlim = c(0.5, 40),
- main = paste("Axon Type", axCount), yaxt = "n", xaxt = "n", xlab = "", ylab = "")
- lines(probSeq, probList[[i + 1]], col = "red", lwd = 1.5)
- lines(probSeq, probList[[i + 2]], col = "blue", lwd = 1.5)
- title(ylab = "Cumulative probability", xlab="% Error", line=3)
-
- }
-
- # Add x axis
- axis(side = 1, at=seq(-10,40,10), labels = seq(-10,40,10))
-
- # Different y axis for Type1 and Type4
- if (axCount == 1) {
-
- axis(side = 2, at=seq(-0.2,1,0.2), labels = seq(-0.2,1,0.2))
-
- } else if (axCount == 4) {
-
- axis(side = 1, at=seq(-10,40,10), labels = seq(-10,40,10))
- axis(side = 2, at=seq(-0.2,1,0.2), labels = seq(-0.2,1,0.2))
-
- }
-
- # Plot vertical segments at 5% and 10% for XY plane
- segments(5, -0.1, 5, h5[1], lty = "dashed", col = "darkgray")
- segments(10, -0.1, 10, h10[1], lty = "dashed", col = "darkgray")
-
- # Plot horizontal segments at cumul prob for 5% and 10% for XY plane
- transp <- 0.8
- segments(-2, h5[1], 5, h5[1], lty = "dashed", col = alpha("darkgray", transp))
- segments(-2, h10[1], 10, h10[1], lty = "dashed", col = alpha("darkgray", transp))
-
- # Store to check later
- hlist[[paste("Type", axCount, sep="")]] <- data.frame(h5 = h5, h10 = h10)
- axCount <- axCount + 1
- }
- ###############
- # PLOT INSETS #
- ###############
- # Find mean bandwidth
- bwFrame <- aggregate(pe ~ neurType + Plane, peFrame, function(x) density(x)$bw)
- for (i in unique(peFrame$neurType)) {
-
- # Estimate density
- dens <- density(peFrame[peFrame$neurType == i & peFrame$Plane == "XY", "pe"], bw = mean(bwFrame$pe))
- # 5%
- x1 <- min(which(dens$x >= -5))
- x2 <- max(which(dens$x < 5))
-
- par(mar = c(0.2,5,0.2,2))
- plot(dens, ylim = c(0, 0.15), xlim = c(-20,20), zero.line = FALSE, main = "", ylab = "", xlab = "", axes = FALSE)
- with(dens, polygon(x=c(x[c(x1,x1:x2,x2)]), y= c(0, y[x1:x2], 0), col="limegreen", border = NA))
- lines(dens)
-
- # Add y axis
- axis(side = 2, at=seq(-0.05,0.15,0.05), labels = seq(-0.05,0.15,0.05), cex.axis = 1.2)
-
-
- # 10%
- x1 <- min(which(dens$x >= -10))
- x2 <- max(which(dens$x < 10))
-
- plot(dens, ylim = c(0, 0.15), xlim = c(-20,20), zero.line = FALSE, main = "", ylab = "", xlab = "", axes = FALSE)
- with(dens, polygon(x=c(x[c(x1,x1:x2,x2)]), y= c(0, y[x1:x2], 0), col="limegreen", border = NA))
- lines(dens)
-
- # Add y axis
- axis(side = 2, at=seq(-0.05,0.15,0.05), labels = seq(-0.05,0.15,0.05), cex.axis = 1.2)
- # Add x axis
- axis(side = 1, at=seq(-30,30,10), labels = seq(-30,30,10), cex.axis = 1.2)
-
-
- }
- dev.off()
- ```
|