modelBased_projections_errors_CV.Rmd 25 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757758759760761762763764765766767768769770771772773774775776777778779780781782783784785786787788789790791792793794795796797798799800801802803804805806807808809810811812813814815816817818819820821822823824825826827828829830831832833834835836837838839840841842843844845846
  1. ---
  2. title: "Axonal Length"
  3. author: "Sergio D?ez Hermano"
  4. date: '`r format(Sys.Date(),"%e de %B, %Y")`'
  5. output:
  6. html_document:
  7. highlight: tango
  8. toc: yes
  9. toc_depth: 4
  10. toc_float:
  11. collapsed: no
  12. ---
  13. ```{r setup, include=FALSE}
  14. require(knitr)
  15. # include this code chunk as-is to set options
  16. 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)
  17. Sys.setlocale("LC_TIME", "C")
  18. ```
  19. ##Load libraries and functions
  20. ```{r}
  21. library(R.matlab)
  22. library(lattice)
  23. library(dplyr)
  24. library(ggplot2)
  25. library(ggridges)
  26. library(hrbrthemes)
  27. library(viridis)
  28. library(ggpubr)
  29. library(mltools)
  30. library(data.table)
  31. library(caret)
  32. library(interactions)
  33. library(performance)
  34. library(MASS)
  35. library(geepack)
  36. library(pstools)
  37. library(MXM)
  38. library(rmcorr)
  39. library(multcomp)
  40. library(Hmisc)
  41. options(digits = 10)
  42. # https://rstudio-pubs-static.s3.amazonaws.com/297778_5fce298898d64c81a4127cf811a9d486.html
  43. ################
  44. # Bootstrap PI #
  45. ################
  46. the.replication <- function(reg, s, x_Np1 = 0){
  47. # Make bootstrap residuals
  48. ep.star <- sample(s, size=length(reg$residuals),replace=TRUE)
  49. # Make bootstrap Y
  50. y.star <- fitted(reg) + ep.star
  51. # Do bootstrap regression
  52. lp <- model.frame(reg)[,2]
  53. nt <- model.frame(reg)[,3]
  54. bs.reg <- lm(y.star ~ lp:nt - 1)
  55. # Create bootstrapped adjusted residuals
  56. bs.lev <- influence(bs.reg)$hat
  57. bs.s <- residuals(bs.reg)/sqrt(1-bs.lev)
  58. bs.s <- bs.s - mean(bs.s)
  59. # Calculate draw on prediction error
  60. # xb.xb <- coef(my.reg)["(Intercept)"] - coef(bs.reg)["(Intercept)"]
  61. xb.xb <- (coef(my.reg)[1] - coef(bs.reg)[1])*x_Np1
  62. return(unname(xb.xb + sample(bs.s, size=1)))
  63. }
  64. ############################
  65. # Bootstrap PI generalized #
  66. ############################
  67. the.replication.gen <- function(reg, s, axType, x_Np1 = 0){
  68. # Make bootstrap residuals
  69. ep.star <- sample(s, size=length(reg$residuals), replace=TRUE)
  70. # Make bootstrap Y
  71. y.star <- fitted(reg) + ep.star
  72. # Do bootstrap regression
  73. lp <- model.frame(reg)[,2]
  74. nt <- model.frame(reg)[,3]
  75. bs.reg <- lm(y.star ~ lp:nt - 1)
  76. # Create bootstrapped adjusted residuals
  77. bs.lev <- influence(bs.reg)$hat
  78. bs.s <- residuals(bs.reg)/sqrt(1-bs.lev)
  79. bs.s <- bs.s - mean(bs.s)
  80. # Calculate draw on prediction error
  81. # xb.xb <- coef(my.reg)["(Intercept)"] - coef(bs.reg)["(Intercept)"]
  82. xb.xb <- (coef(my.reg)[axType] - coef(bs.reg)[axType])*x_Np1
  83. return(unname(xb.xb + sample(bs.s, size=1)))
  84. }
  85. ##############################
  86. # Bootstrap PI per axon type #
  87. ##############################
  88. the.replication.type <- function(reg, s, x_Np1 = 0){
  89. # Make bootstrap residuals
  90. ep.star <- sample(s, size=length(reg$residuals),replace=TRUE)
  91. # Make bootstrap Y
  92. y.star <- fitted(reg) + ep.star
  93. # Do bootstrap regression
  94. lp <- model.frame(reg)[,2]
  95. bs.reg <- lm(y.star ~ lp - 1)
  96. # Create bootstrapped adjusted residuals
  97. bs.lev <- influence(bs.reg)$hat
  98. bs.s <- residuals(bs.reg)/sqrt(1-bs.lev)
  99. bs.s <- bs.s - mean(bs.s)
  100. # Calculate draw on prediction error
  101. # xb.xb <- coef(my.reg)["(Intercept)"] - coef(bs.reg)["(Intercept)"]
  102. xb.xb <- (coef(my.reg) - coef(bs.reg))*x_Np1
  103. return(unname(xb.xb + sample(bs.s, size=1)))
  104. }
  105. ##############################
  106. # Stratified random sampling #
  107. ##############################
  108. stratified <- function(df, group, size, select = NULL,
  109. replace = FALSE, bothSets = FALSE) {
  110. if (is.null(select)) {
  111. df <- df
  112. } else {
  113. if (is.null(names(select))) stop("'select' must be a named list")
  114. if (!all(names(select) %in% names(df)))
  115. stop("Please verify your 'select' argument")
  116. temp <- sapply(names(select),
  117. function(x) df[[x]] %in% select[[x]])
  118. df <- df[rowSums(temp) == length(select), ]
  119. }
  120. df.interaction <- interaction(df[group], drop = TRUE)
  121. df.table <- table(df.interaction)
  122. df.split <- split(df, df.interaction)
  123. if (length(size) > 1) {
  124. if (length(size) != length(df.split))
  125. stop("Number of groups is ", length(df.split),
  126. " but number of sizes supplied is ", length(size))
  127. if (is.null(names(size))) {
  128. n <- setNames(size, names(df.split))
  129. message(sQuote("size"), " vector entered as:\n\nsize = structure(c(",
  130. paste(n, collapse = ", "), "),\n.Names = c(",
  131. paste(shQuote(names(n)), collapse = ", "), ")) \n\n")
  132. } else {
  133. ifelse(all(names(size) %in% names(df.split)),
  134. n <- size[names(df.split)],
  135. stop("Named vector supplied with names ",
  136. paste(names(size), collapse = ", "),
  137. "\n but the names for the group levels are ",
  138. paste(names(df.split), collapse = ", ")))
  139. }
  140. } else if (size < 1) {
  141. n <- round(df.table * size, digits = 0)
  142. } else if (size >= 1) {
  143. if (all(df.table >= size) || isTRUE(replace)) {
  144. n <- setNames(rep(size, length.out = length(df.split)),
  145. names(df.split))
  146. } else {
  147. message(
  148. "Some groups\n---",
  149. paste(names(df.table[df.table < size]), collapse = ", "),
  150. "---\ncontain fewer observations",
  151. " than desired number of samples.\n",
  152. "All observations have been returned from those groups.")
  153. n <- c(sapply(df.table[df.table >= size], function(x) x = size),
  154. df.table[df.table < size])
  155. }
  156. }
  157. temp <- lapply(
  158. names(df.split),
  159. function(x) df.split[[x]][sample(df.table[x],
  160. n[x], replace = replace), ])
  161. set1 <- do.call("rbind", temp)
  162. if (isTRUE(bothSets)) {
  163. set2 <- df[!rownames(df) %in% rownames(set1), ]
  164. list(SET1 = set1, SET2 = set2)
  165. } else {
  166. set1
  167. }
  168. }
  169. ```
  170. ##Load data
  171. ```{r}
  172. # Get file paths
  173. axonNames <- list.files(paste("ProyeccionData/", sep=""), full.names=T)
  174. # Load matlab file
  175. axonFiles <- lapply(axonNames, function(x) readMat(x))
  176. names(axonFiles) <- c("Type1", "Type2", "Type3", "Type4")
  177. # VAR NAMES
  178. # REAL_LENGTH, PROYECTION_LENGTH_XY, PROYECTION_LENGTH_XZ, PROYECTION_LENGTH_YZ
  179. # Extract data
  180. realLength <- lapply(axonFiles, function(x) x$REAL.LENGTH)
  181. planeXY <- lapply(axonFiles, function(x) x$PROYECTION.LENGTH.XY)
  182. planeXZ <- lapply(axonFiles, function(x) x$PROYECTION.LENGTH.XZ)
  183. planeYZ <- lapply(axonFiles, function(x) x$PROYECTION.LENGTH.YZ)
  184. # Get number of neurons per type
  185. numberTypes <- unlist(lapply(realLength, function(x) length(x)))
  186. # Store in data frames by plane
  187. xyData <- data.frame(id = 1:length(unlist(realLength)),
  188. LR = unlist(realLength), LP = unlist(planeXY), Plane = rep("XY", sum(numberTypes)),
  189. neurType = rep(c("Type1", "Type2", "Type3", "Type4"), numberTypes))
  190. xzData <- data.frame(id = 1:length(unlist(realLength)),
  191. LR = unlist(realLength), LP = unlist(planeXZ), Plane = rep("XZ", sum(numberTypes)),
  192. neurType = rep(c("Type1", "Type2", "Type3", "Type4"), numberTypes))
  193. yzData <- data.frame(id = 1:length(unlist(realLength)),
  194. LR = unlist(realLength), LP = unlist(planeYZ), Plane = rep("YZ", sum(numberTypes)),
  195. neurType = rep(c("Type1", "Type2", "Type3", "Type4"), numberTypes))
  196. # Merge into a single data frame and sort by id
  197. allData <- do.call("rbind", list(xyData, xzData, yzData))
  198. allData <- allData[order(allData$id),]
  199. # Add different number for every neuron-plane combination
  200. allData$NeurPlane <- c(rep(c(1,2,3), numberTypes[1]), rep(c(4,5,6), numberTypes[2]),
  201. rep(c(7,8,9), numberTypes[3]), rep(c(10,11,12), numberTypes[4]))
  202. allData$interact <-interaction(allData$neurType, allData$Plane, lex.order=T)
  203. # Create data frame w/o Type4
  204. dataNo4 <- allData[allData$neurType != "Type4", ]; dataNo4$neurType <- factor(dataNo4$neurType)
  205. # Data in short format
  206. dataShort <- data.frame(id = 1:length(unlist(realLength)),
  207. LR = unlist(realLength), XY = unlist(planeXY), XZ = unlist(planeXZ), YZ = unlist(planeYZ),
  208. neurType = rep(c("Type1", "Type2", "Type3", "Type4"), numberTypes))
  209. ```
  210. ##CV-Bootstrap Pred Interval
  211. This code is meant to be run independently for every neurType and in parallel. RUN AT CREMOTO
  212. ```{r}
  213. # Ensure reproducibility
  214. set.seed(123456)
  215. # Subset axon type
  216. dataCV <- allData[allData$neurType == "Type1", ]
  217. # List to store CV repetitions
  218. cvList <- list ()
  219. cvRand <- list()
  220. for (i in 1:10) {
  221. # Stratified random training/test sets
  222. stratDat <- stratified(dataCV, "Plane", 0.7, bothSets = TRUE, replace = FALSE)
  223. trainSet <- stratDat$SET1
  224. testSet <- stratDat$SET2
  225. # Store train and test set
  226. cvRand[[i]] <- list(trainSet, testSet)
  227. # OLS with training set
  228. my.reg <- lm(LR ~ LP - 1, trainSet)
  229. # Create adjusted residuals
  230. leverage <- influence(my.reg)$hat
  231. my.s.resid <- residuals(my.reg)/sqrt(1-leverage)
  232. my.s.resid <- my.s.resid - mean(my.s.resid)
  233. # List to store planes
  234. planeList <- list()
  235. planeRand <- list()
  236. # Assign test data to bootstrapping
  237. dataToBoot <- testSet
  238. # Bootstrapping on test set
  239. for (j in c("XY", "XZ", "YZ")) {
  240. # Subset plane
  241. dataSub <- dataToBoot[dataToBoot$Plane == j, ]
  242. # List to store errors
  243. epList <- list()
  244. epCount <- 1
  245. for (k in dataSub[ , "LP"]) {
  246. # Predict LR
  247. y.p <- coef(my.reg)*k
  248. # Do bootstrap with 100 replications
  249. ep.draws <- replicate(n=100, the.replication.type(reg = my.reg, s = my.s.resid, x_Np1 = k))
  250. # Store in list
  251. epList[[epCount]] <- ep.draws
  252. epCount <- epCount + 1
  253. print(epCount)
  254. }
  255. # Store in list
  256. planeList[[j]] <- epList
  257. print(paste("CV", i, "Plane", j, "finished"))
  258. }
  259. # Store in list
  260. cvList[[i]] <- planeList
  261. print(paste("CV", i, "finished"))
  262. }
  263. # names(axList) <- c("Type1", "Type2", "Type3", "Type4")
  264. # names(axRand) <- c("Type1", "Type2", "Type3", "Type4")
  265. #
  266. # saveRDS(list(distList, distRand), paste("stereo_bootPI_full_Type", axCount, "_boot20.rds", sep = ""))
  267. ```
  268. ####Check convergence
  269. ```{r}
  270. # Load object
  271. bootCheck <- readRDS("bootPI_full.rds")
  272. # First list is prediction error, second is random subset
  273. axList <- bootCheck[[1]]
  274. axRand <- bootCheck[[2]]
  275. # One example
  276. # OLS
  277. my.reg <- lm(LR ~ LP:neurType - 1, allData)
  278. # Estimate and store in list
  279. lrpList <- list()
  280. lrpCount <- 1
  281. axCount <- 1
  282. for (i in c("Type1", "Type2", "Type3", "Type4")) {
  283. for (j in c("XY", "XZ", "YZ")) {
  284. # Extract LP and LR values
  285. lp <- axRand[[i]][[j]]$LP
  286. lr <- axRand[[i]][[j]]$LR
  287. for (k in 1:length(lp)) {
  288. # Predicted LR
  289. lr.p <- coef(my.reg)[axCount]*lp[k]
  290. # Extract absolute errors
  291. ae <- axList[[i]][[j]][[k]]
  292. # Empty vectors for storage
  293. lrp.Low <- numeric()
  294. lrp.Up <- numeric()
  295. # Estimate PI95 for boot n=10...100
  296. for (l in 4:100) {
  297. lrp <- lr.p + quantile(ae[1:l], probs = c(0.025,0.975))
  298. # Store in vectors
  299. lrp.Low <- c(lrp.Low, lrp[1])
  300. lrp.Up <- c(lrp.Up, lrp[2])
  301. }
  302. # Store in data frame
  303. lrpFrame <- data.frame(nboot = 4:100, lowLim = lrp.Low, upLim = lrp.Up, LR = rep(lr[k], length(4:100)))
  304. # Store in list
  305. lrpList[[lrpCount]] <- lrpFrame
  306. lrpCount <- lrpCount + 1
  307. }
  308. }
  309. axCount <- axCount + 1
  310. }
  311. ```
  312. #####Reformat
  313. ```{r}
  314. # 1/3 random observations per plane, no replacement
  315. numExtract <- floor(1/3*numberTypes)
  316. checkFrame <- do.call(rbind, lrpList)
  317. checkFrame$neurType <- rep(c("Type1", "Type2", "Type3", "Type4"), times = numExtract*(length(4:100))*3)
  318. checkFrame$Plane <- c(rep(c("XY", "XZ", "YZ"), each = numExtract[1]*length(4:100)),
  319. rep(c("XY", "XZ", "YZ"), each = numExtract[2]*length(4:100)),
  320. rep(c("XY", "XZ", "YZ"), each = numExtract[3]*length(4:100)),
  321. rep(c("XY", "XZ", "YZ"), each = numExtract[4]*length(4:100)))
  322. lowMean <- aggregate(lowLim ~ neurType + Plane + nboot, checkFrame, mean)
  323. upMean <- aggregate(upLim ~ neurType + Plane + nboot, checkFrame, mean)
  324. LRmean <- aggregate(LR ~ neurType + Plane + nboot, checkFrame, mean)
  325. cumError <- with(lowMean[lowMean$neurType == "Type1" & lowMean$Plane == "XY", ],
  326. abs((lowLim[-1] - lowLim[-length(lowLim)])/lowLim[-length(lowLim)])*100)
  327. ```
  328. #####Plot
  329. ```{r, fig.height=20, fig.width=30}
  330. # Separate plots
  331. png(filename=paste("convergence_bootPI_separate.png", sep=""), type="cairo",
  332. units="in", width=30, height=20, pointsize=12,
  333. res=300)
  334. par(mfrow=c(4,6), cex = 1.1)
  335. for (i in c("Type1", "Type2", "Type3", "Type4")) {
  336. for (j in c("XY", "XZ", "YZ")) {
  337. plot(lowLim ~ nboot, lowMean[lowMean$neurType == i & lowMean$Plane == j, ], type = "l", lwd = 2)
  338. plot(upLim ~ nboot, upMean[upMean$neurType == i & upMean$Plane == j, ], type = "l", lwd = 2)
  339. }
  340. }
  341. dev.off()
  342. # Same plot
  343. png(filename=paste("convergence_bootPI_together.png", sep=""), type="cairo",
  344. units="in", width=15, height=20, pointsize=12,
  345. res=300)
  346. par(mfrow=c(4,3), cex = 1.1)
  347. for (i in c("Type1", "Type2", "Type3", "Type4")) {
  348. for (j in c("XY", "XZ", "YZ")) {
  349. minY <- min(lowMean[lowMean$neurType == i & lowMean$Plane == j, "lowLim"])
  350. maxY <- max(upMean[upMean$neurType == i & upMean$Plane == j, "upLim"])
  351. plot(lowLim ~ nboot, lowMean[lowMean$neurType == i & lowMean$Plane == j, ], type = "l", lwd = 2, ylim = c(minY, maxY))
  352. lines(upLim ~ nboot, upMean[upMean$neurType == i & upMean$Plane == j, ], type = "l", lwd = 2)
  353. # lines(LR ~ nboot, LRmean[LRmean$neurType == i & LRmean$Plane == j, ], type = "l", lwd = 2, col = "red")
  354. }
  355. }
  356. dev.off()
  357. # Diff plot
  358. png(filename=paste("convergence_bootPI_Diff.png", sep=""), type="cairo",
  359. units="in", width=15, height=20, pointsize=12,
  360. res=300)
  361. par(mfrow=c(4,3), cex = 1.1)
  362. for (i in c("Type1", "Type2", "Type3", "Type4")) {
  363. for (j in c("XY", "XZ", "YZ")) {
  364. minY <- min(lowMean[lowMean$neurType == i & lowMean$Plane == j, "lowLim"])
  365. maxY <- max(upMean[upMean$neurType == i & upMean$Plane == j, "upLim"])
  366. diflu <- upMean[upMean$neurType == i & upMean$Plane == j, "upLim"] - lowMean[lowMean$neurType == i & lowMean$Plane == j, "lowLim"]
  367. plot(4:100, diflu, type = "l", col = "red", lwd = 2)
  368. }
  369. }
  370. dev.off()
  371. # Same plot random intervals
  372. for (m in 1:4) {
  373. png(filename=paste("convergence_bootPI_together_example_rand", m, ".png", sep=""), type="cairo",
  374. units="in", width=15, height=20, pointsize=12,
  375. res=300)
  376. par(mfrow=c(4,3), cex = 1.1)
  377. for (i in c("Type1", "Type2", "Type3", "Type4")) {
  378. for (j in c("XY", "XZ", "YZ")) {
  379. subData <- checkFrame[checkFrame$neurType == i & checkFrame$Plane == j, ]
  380. randLR <- sample(subData[, "LR"], 1)
  381. # randLR <- sample(subData[which(subData$lowLim > subData$LR), "LR"], 1)
  382. minY <- min(subData[subData$LR == randLR, "lowLim"])
  383. maxY <- max(subData[subData$LR == randLR, "upLim"])
  384. plot(lowLim ~ nboot, subData[subData$LR == randLR, ], type = "l", lwd = 2, ylim = c(minY, maxY))
  385. lines(upLim ~ nboot, subData[subData$LR == randLR, ], type = "l", lwd = 2)
  386. lines(LR ~ nboot, subData[subData$LR == randLR, ], type = "l", lwd = 2, col = "red")
  387. }
  388. }
  389. dev.off()
  390. }
  391. # Separate with lines
  392. png(filename=paste("convergence_bootPI_separate_FULL.png", sep=""), type="cairo",
  393. units="in", width=30, height=20, pointsize=12,
  394. res=300)
  395. par(mfrow=c(4,6), cex = 1.1)
  396. for (i in c("Type1", "Type2", "Type3", "Type4")) {
  397. for (j in c("XY", "XZ", "YZ")) {
  398. plot(lowLim ~ nboot, checkFrame[checkFrame$neurType == i & checkFrame$Plane == j, ], type = "l", col = "lightgray")
  399. lines(lowLim ~ nboot, lowMean[lowMean$neurType == i & lowMean$Plane == j, ], type = "l", lwd = 2)
  400. plot(upLim ~ nboot, checkFrame[checkFrame$neurType == i & checkFrame$Plane == j, ], type = "l", col = "lightgray")
  401. lines(upLim ~ nboot, upMean[upMean$neurType == i & upMean$Plane == j, ], type = "l", lwd = 2)
  402. }
  403. }
  404. dev.off()
  405. ```
  406. ####Estimate porcentual errors
  407. ```{r, fig.width=8, fig.height=8}
  408. # # Load objects
  409. cvPI <- list()
  410. for (i in c("Type1", "Type2", "Type3", "Type4")) {
  411. # Whithin each Type, element 1 are erros and element 2 are random samples, training and test
  412. cvPI[[i]] <- readRDS(paste("ProyeccionAnalysis/proyec_boot10CVPI_", i, ".rds", sep = ""))
  413. }
  414. # Estimate and store in list
  415. peList <- list()
  416. peCount <- 1
  417. lrLength <- list()
  418. for (h in c("Type1", "Type2", "Type3", "Type4")) {
  419. for (i in 1:10) {
  420. for (j in c("XY", "XZ", "YZ")) {
  421. # Extract LR values
  422. lr <- cvPI[[h]][[2]][[i]][[2]]
  423. lr <- lr[lr$Plane == j, "LR"]
  424. for (k in 1:length(lr)) {
  425. # Divide absolute errors by LR to get porcentual errors
  426. pe <- cvPI[[h]][[1]][[i]][[j]][[k]] / lr[k]
  427. # Store in list
  428. peList[[peCount]] <- pe*100
  429. peCount <-peCount + 1
  430. }
  431. }
  432. }
  433. lrLength[h] <- length(lr)
  434. }
  435. lrLength <- unlist(lrLength)
  436. peFrame <- data.frame(pe = unlist(peList),
  437. neurType = rep(c("Type1", "Type2", "Type3", "Type4"), times = lrLength*10*3*100),
  438. Plane = c(rep(rep(c("XY", "XZ", "YZ"), each = lrLength[1]*100), 10),
  439. rep(rep(c("XY", "XZ", "YZ"), each = lrLength[2]*100), 10),
  440. rep(rep(c("XY", "XZ", "YZ"), each = lrLength[3]*100), 10),
  441. rep(rep(c("XY", "XZ", "YZ"), each = lrLength[4]*100), 10)))
  442. peFrame$interact <- interaction(peFrame$neurType, peFrame$Plane, lex.order=T)
  443. ```
  444. #####Plot density and histogram
  445. ```{r}
  446. # Plot
  447. # peGroup <- peFrame %>% group_by(interact)
  448. # peSample <- sample_n(peGroup, 100)
  449. setwd("ProyeccionAnalysis/")
  450. png(filename=paste("prederrorDistributions_ols_CV.png", sep=""), type="cairo",
  451. units="in", width=8, height=10, pointsize=12,
  452. res=300)
  453. (p <- ggplot(peFrame, aes(x = `pe`, y = `interact`, fill = ifelse(..x.. >=-5 & ..x.. <=5, ">= -5", "< 5"))) +
  454. scale_fill_manual(values = c("gray70", "coral2"), name = NULL) +
  455. geom_density_ridges_gradient(alpha = 0.8) +
  456. theme_ridges() +
  457. theme(legend.position = "none") +
  458. xlab("% Error")) +
  459. xlim(c(-40,40))
  460. dev.off()
  461. png(filename=paste("prederrorHistogram_ols_CV.png", sep=""), type="cairo",
  462. units="in", width=8, height=10, pointsize=12,
  463. res=300)
  464. # Histogram
  465. # nint 300 bootPI
  466. # nint 500 bootPI_full
  467. # nint 400 bootPI_full_replace
  468. (h <- histogram( ~ pe | Plane + neurType, peFrame, nint = 300, type ="density", xlim = c(-50,50)))
  469. dev.off()
  470. ```
  471. #####Plot errors vs cumul probability
  472. ```{r, fig.height=5, fig.width=20}
  473. # Inverse quantile
  474. quantInv <- function(distr, value) ecdf(distr)(value)
  475. #ecdf plot example
  476. # plot(ecdf(peFrame[peFrame$interact == "Type1.XY", "pe"]))
  477. # lines(ecdf(peFrame[peFrame$interact == "Type1.XZ", "pe"]), col = "red")
  478. # lines(ecdf(peFrame[peFrame$interact == "Type1.YZ", "pe"]), col = "blue")
  479. probList <- list()
  480. binProb <- 0.5
  481. probSeq <- seq(binProb, 100, binProb)
  482. for (i in unique(peFrame$interact)) {
  483. dataProb <- peFrame[peFrame$interact == i, "pe"]
  484. probVec <- numeric()
  485. for (j in probSeq) {
  486. errProb <- quantInv(dataProb, j) - quantInv(dataProb, -j)
  487. probVec <- c(probVec, errProb)
  488. }
  489. probList[[i]] <- probVec
  490. }
  491. # Plot with limited axis
  492. setwd("ProyeccionAnalysis/")
  493. png(filename=paste("errorProbabilityCum_bin", binProb, "_CV.png", sep=""), type="cairo",
  494. units="in", width=20, height=5, pointsize=12,
  495. res=300)
  496. par(mfrow = c(1,4))
  497. axCount <- 1
  498. for (i in seq(1,10,3)) {
  499. plot(probSeq, probList[[i]], type = "l", col = "green", lwd = 1.5, ylim = c(0,1), xlim = c(0.5, 40),
  500. ylab = "Cumulative probability", xlab = "% Error", main = paste("Axon Type", axCount))
  501. lines(probSeq, probList[[i + 1]], col = "red", lwd = 1.5)
  502. lines(probSeq, probList[[i + 2]], col = "blue", lwd = 1.5)
  503. abline(v = 5, lty = "dashed", col = "gray")
  504. abline(h = 0.95, lty = "dashed", col = "gray" )
  505. legend("bottomright", c("XY", "XZ", "YZ"), lwd = 1.5, col = c("green", "red", "blue"))
  506. axCount <- axCount + 1
  507. }
  508. dev.off()
  509. # Plot with free axis
  510. png(filename=paste("errorProbabilityCum_bin", binProb, "_freeAxis_CV.png", sep=""), type="cairo",
  511. units="in", width=20, height=5, pointsize=12,
  512. res=300)
  513. par(mfrow = c(1,4))
  514. axCount <- 1
  515. for (i in seq(1,10,3)) {
  516. plot(probSeq, probList[[i]], type = "l", col = "green", lwd = 1.5,
  517. ylab = "Cumulative probability", xlab = "% Error", main = paste("Axon Type", axCount))
  518. lines(probSeq, probList[[i + 1]], col = "red", lwd = 1.5)
  519. lines(probSeq, probList[[i + 2]], col = "blue", lwd = 1.5)
  520. abline(v = 5, lty = "dashed", col = "gray")
  521. abline(h = 0.95, lty = "dashed", col = "gray" )
  522. legend("bottomright", c("XY", "XZ", "YZ"), lwd = 1.5, col = c("green", "red", "blue"))
  523. axCount <- axCount + 1
  524. }
  525. dev.off()
  526. ```
  527. #####Random sample
  528. ```{r}
  529. size <- 0.05
  530. peRand <- peFrame %>% group_by(interact)
  531. peSample <- sample_frac(peRand, size)
  532. ```
  533. ######Plot density and histogram
  534. ```{r}
  535. setwd("ProyeccionAnalysis/")
  536. png(filename=paste("prederrorDistributions_ols_CV_rand", size, ".png", sep=""), type="cairo",
  537. units="in", width=8, height=10, pointsize=12,
  538. res=300)
  539. (p <- ggplot(peSample, aes(x = `pe`, y = `interact`, fill = ifelse(..x.. >=-5 & ..x.. <=5, ">= -5", "< 5"))) +
  540. scale_fill_manual(values = c("gray70", "coral2"), name = NULL) +
  541. geom_density_ridges_gradient(alpha = 0.8) +
  542. theme_ridges() +
  543. theme(legend.position = "none") +
  544. xlab("% Error")) +
  545. xlim(c(-40,40))
  546. dev.off()
  547. png(filename=paste("prederrorHistogram_ols_CV_rand", size, ".png", sep=""), type="cairo",
  548. units="in", width=8, height=10, pointsize=12,
  549. res=300)
  550. # Histogram
  551. # nint 300 bootPI
  552. # nint 500 bootPI_full
  553. # nint 400 bootPI_full_replace
  554. (h <- histogram( ~ pe | Plane + neurType, peSample, nint = 300, type ="density", xlim = c(-50,50)))
  555. dev.off()
  556. ```
  557. ######Plot errors vs cumul probability
  558. ```{r, fig.height=5, fig.width=20}
  559. # Inverse quantile
  560. quantInv <- function(distr, value) ecdf(distr)(value)
  561. peSample <- as.data.frame(peSample)
  562. probList <- list()
  563. binProb <- 0.5
  564. probSeq <- seq(binProb, 100, binProb)
  565. for (i in unique(peSample$interact)) {
  566. dataProb <- peSample[peSample$interact == i, "pe"]
  567. probVec <- numeric()
  568. for (j in probSeq) {
  569. errProb <- quantInv(dataProb, j) - quantInv(dataProb, -j)
  570. probVec <- c(probVec, errProb)
  571. }
  572. probList[[i]] <- probVec
  573. }
  574. # Plot with limited axis
  575. setwd("ProyeccionAnalysis/")
  576. png(filename=paste("errorProbabilityCum_bin", binProb, "_CV_rand", size, ".png", sep=""), type="cairo",
  577. units="in", width=20, height=5, pointsize=12,
  578. res=300)
  579. par(mfrow = c(1,4))
  580. axCount <- 1
  581. for (i in seq(1,10,3)) {
  582. plot(probSeq, probList[[i]], type = "l", col = "green", lwd = 1.5, ylim = c(0,1), xlim = c(0.5, 40),
  583. ylab = "Cumulative probability", xlab = "% Error", main = paste("Axon Type", axCount))
  584. lines(probSeq, probList[[i + 1]], col = "red", lwd = 1.5)
  585. lines(probSeq, probList[[i + 2]], col = "blue", lwd = 1.5)
  586. abline(v = 5, lty = "dashed", col = "gray")
  587. abline(h = 0.95, lty = "dashed", col = "gray" )
  588. legend("bottomright", c("XY", "XZ", "YZ"), lwd = 1.5, col = c("green", "red", "blue"))
  589. axCount <- axCount + 1
  590. }
  591. dev.off()
  592. # Plot with free axis
  593. png(filename=paste("errorProbabilityCum_bin", binProb, "_freeAxis_CV_rand", size, ".png", sep=""), type="cairo",
  594. units="in", width=20, height=5, pointsize=12,
  595. res=300)
  596. par(mfrow = c(1,4))
  597. axCount <- 1
  598. for (i in seq(1,10,3)) {
  599. plot(probSeq, probList[[i]], type = "l", col = "green", lwd = 1.5,
  600. ylab = "Cumulative probability", xlab = "% Error", main = paste("Axon Type", axCount))
  601. lines(probSeq, probList[[i + 1]], col = "red", lwd = 1.5)
  602. lines(probSeq, probList[[i + 2]], col = "blue", lwd = 1.5)
  603. abline(v = 5, lty = "dashed", col = "gray")
  604. abline(h = 0.95, lty = "dashed", col = "gray" )
  605. legend("bottomright", c("XY", "XZ", "YZ"), lwd = 1.5, col = c("green", "red", "blue"))
  606. axCount <- axCount + 1
  607. }
  608. dev.off()
  609. ```
  610. ####Estimate mean porcentual errors
  611. UNFINISHED
  612. ```{r, fig.width=8, fig.height=8}
  613. # mm <- peFrame[peFrame$Plane == "XY", "pe"]
  614. # mm.mat <- matrix(mm, 100, 270)
  615. # meanMat <- apply(mm.mat, 1, mean)
  616. ```