Rubio_etal_2021/Code/Statistical Analysis/dataforPepo.Rmd

---
title: "Axonal Length"
author: "Sergio D?ez Hermano"
date: '`r format(Sys.Date(),"%e de %B, %Y")`'
output:
  html_document:
      highlight: tango
      toc: yes
      toc_depth: 4
      toc_float:
        collapsed: no
---

```{r setup, include=FALSE}
require(knitr)
# include this code chunk as-is to set options
opts_chunk$set(comment = NA, prompt = FALSE, fig.height=5, fig.width=5, dpi=300, fig.align = "center", echo = TRUE, message = FALSE, warning = FALSE, cache=FALSE)
Sys.setlocale("LC_TIME", "C")
```

##Libraries and functions

```{r}
suppressMessages(library(R.matlab))  #import mat files
```


##PROJECTIONS

###Cumulative probability

####Format and save

```{r}
# Load rds
probList <- readRDS("ProyeccionAnalysis/cumProb_5CV.rds")
binProb <- 0.5
probSeq <- seq(binProb, 100, binProb)

# Store as long format data frame
probFrame <- data.frame(axonType = rep(c("Type1", "Type2", "Type3", "Type4"), each = 3*200),
                        Plane = rep(rep(c("XY", "XZ", "YZ"), each = 200),4),
                        error = rep(probSeq, 12),
                        cumulProb = unlist(probList))

# Save as csv
write.csv(probFrame, "dataPepo/projection_error_cumulativeProbability.csv", sep=",", row.names = FALSE)
```

###Error density

####Estimation

```{r}
# Load 5-fold CV
axList <- readRDS("ProyeccionAnalysis/modelBased_projections_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:length(lr)) {
        
        # Divide errors by LR to get porcentual errors
        pe <- lp[k, ] / lr[k]
        
        # Store in list
        peList[[peCount]] <- pe*100
        peCount <-peCount + 1
      }
      
    }
    
  }
  
  lrLength[h] <- length(lr)
  
}

# 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)

# Find mean bandwidth
bwFrame <- aggregate(pe ~ neurType + Plane, peFrame, function(x) density(x)$bw)

# Estimate density
densList <- by(peFrame$pe, peFrame$interact, function(x) density(x, bw = mean(bwFrame$pe)))
```

####Format and save

```{r}
# Extract from list
densXY <- lapply(densList, function(x) cbind(error = x$x, density = x$y))
densFrame <- data.frame(axonType = rep(c("Type1", "Type2", "Type3", "Type4"), each = 512*3),
                        Plane = rep(rep(c("XY", "XZ", "YZ"), each = 512), 4))

# Store in long format data frame
densFrame <- cbind(densFrame, do.call(rbind, densXY))

# Save density as csv
write.csv(densFrame, "dataPepo/projection_error_Densities.csv", sep=",", row.names = FALSE)

# Save original data as csv
oriFrame <- cbind(axonType = peFrame[, "neurType"], Plane = peFrame[, "Plane"], error = peFrame[, "pe"], data.frame(abserror = abs(peFrame$pe)))
write.csv(oriFrame, "dataPepo/projection_error_OriginalData.csv", sep=",", row.names = FALSE)
```

###Mean errors

####Estimate, format and save

```{r}
# Mean absolute error
meanErr <- aggregate(abs(pe) ~ Plane + neurType, peFrame, mean)

# SD absolute error
sdErr <- aggregate(abs(pe) ~ Plane + neurType, peFrame, sd)

# SE absolute error
seErr <- aggregate(abs(pe) ~ Plane + neurType, peFrame, sjstats::se)

# Common data frame
meansdErr <- data.frame(axonType = meanErr$neurType, 
                        Plane = meanErr$Plane, 
                        mean = meanErr$`abs(pe)`, 
                        sd = sdErr$`abs(pe)`,
                        se = seErr$`abs(pe)`,
                        ci95 = 1.96*seErr$`abs(pe)`)

# Save as csv
write.csv(meansdErr, "dataPepo/projection_error_MeanSD.csv", sep=",", row.names = FALSE)
```

###Five-number summary

####Format and save
```{r}
# Store boxplot data
boxData <- boxplot(abs(pe) ~ Plane + neurType, peFrame)
fiveNum <- boxData$stats

# Reformat as data frame
fiveNum.frame <- data.frame(axonType = rep(c("Type1", "Type2", "Type3", "Type4"), each = 15),
                            Plane = rep(rep(c("XY", "XZ", "YZ"), each = 5), 4),
                            quantileType = rep(c("lowlim", "p25", "p50", "p75", "uplim"), 12),
                            quantileDat = as.vector(fiveNum))

# Save
write.csv(fiveNum.frame, "dataPepo/projection_error_boxplot.csv", sep=",", row.names = FALSE)
```


##STEREO PLANES

###Load data

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

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

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

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

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

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

###Mean errors

####Estimate, format and save

```{r}
# Estimate mean residual error
meanFrame <- aggregate(error ~ dist + step + neurType, allData, mean)
n <- 200
grid <- expand.grid(dist = seq(3, 30, length = n), step = seq(70, 150, length = n))

# Smooth lm
for (i in unique(meanFrame$neurType)) {
  
  axlm <- lm(error ~ dist * step, data = meanFrame[meanFrame$neurType == i, ])
  print(summary(axlm))
  axfit <- cbind(grid, prediction = predict(axlm, grid))
  
  # Save as csv
  write.csv(axfit, paste("dataPepo/stereoPlanes_meanError_smooth_", i, ".csv", sep=""))
  
}
```

###Error density

####Estimation

```{r}
# Find mean bandwidth
bwFrame <- aggregate(error2 ~ neurType + dist + step, allData, function(x) density(x)$bw)

# Estimate density using mean bw per axonType
densList <- list()

for (i in c("Type1", "Type2", "Type3", "Type4")) {
  
  axData <- allData[allData$neurType == i, ]
  densList[[i]] <- by(axData$error2, axData$interact, function(x) density(x, 
                                                                           bw = ceiling(mean(bwFrame[bwFrame$neurType == i, "error2"]))))
  
}

# Restore format
densList <- unlist(densList, recursive = FALSE)
```


####Format and save

```{r}
# Extract from list
densXY <- lapply(densList, function(x) cbind(error = x$x, density = x$y))
densFrame <- data.frame(axonType = rep(c("Type1", "Type2", "Type3", "Type4"), each = 512*90),
                        step = rep(rep(unique(allData$step), each = 512*10), 4), 
                        step = rep(rep(unique(allData$dist), each = 512), 4*9))

# Store in long format data frame
densFrame <- cbind(densFrame, do.call(rbind, densXY))

# Save density as csv
write.csv(densFrame, "dataPepo/stereoPlanes_error_Densities.csv", sep=",", row.names = FALSE)

# Save original data as csv
oriFrame <- cbind(axonType = allData[, "neurType"], dist = allData[, "dist"], step = allData[, "step"],
                  error = allData[, "error2"], abserror = allData[, "error"])
write.csv(oriFrame, "dataPepo/stereoPlanes_error_OriginalData.csv", sep=",", row.names = FALSE)
```

###Cumul probability

####Estimate, format and save

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


n <- 200
grid <- expand.grid(dist = seq(3, 30, length = n), step = seq(70, 150, length = n))

# Smooth lm
for (h in c(5, 10)) {
  
  for (i in unique(probFrame$neurType)) {
    
    axlm <- lm(prob ~ dist * step, data = probFrame[probFrame$error == h & probFrame$neurType == i, ])
    print(summary(axlm))
    axfit <- cbind(grid, prediction = predict(axlm, grid))
    
    # Save as csv
    write.csv(axfit, paste("dataPepo/stereoPlanes_cumulProb_error_", h, "_smooth_", i, ".csv", sep=""))
    
  }
  
  
}

```

##STEREO SPHERES

###Load data

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

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

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

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

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

# Store in data frames
allData <- data.frame(id = rep(1:sum(numberTypes), each = 81),
                      diams = rep(rep(diams, each = 9), sum(numberTypes)), 
                      step = rep(rep(step, 9), sum(numberTypes)),
                      error = abs(errorVector),
                      error2 = errorVector,
                      LRe = lpVector,
                      LR = rep(lrVector, each = 81))

# Remove all LR = 0 (duplicated Type3)
allData <- allData[allData$LR != 0, ]

# Get number of neurons per type
numberTypes  <- numberTypes[-3]
names(numberTypes) <- c("Type1", "Type2", "Type3", "Type4")

# Reassign id and neurType
allData$id <- rep(1:sum(numberTypes), each = 81)
allData$neurType <- rep(c("Type1", "Type2", "Type3", "Type4"), times = 81*numberTypes)

allData$errorRaw <- allData$LR - allData$LRe
allData$interact <- interaction(allData$step, allData$diams, lex.order = T)
allData$interact2 <- interaction(allData$neurType, allData$step, allData$diams, lex.order = T)
```

###Mean errors

####Estimate, format and save

```{r}
# Estimate mean residual error
meanFrame <- aggregate(error ~ diams + step + neurType, allData, mean)
n <- 200
grid <- expand.grid(diams = seq(10, 50, length = n), step = seq(70, 150, length = n))

# Smooth lm
for (i in unique(meanFrame$neurType)) {
  
  axlm <- lm(error ~ diams * step, data = meanFrame[meanFrame$neurType == i, ])
  print(summary(axlm))
  axfit <- cbind(grid, prediction = predict(axlm, grid))
  
  # Save as csv
  write.csv(axfit, paste("dataPepo/stereoSpheres_meanError_smooth_", i, ".csv", sep=""))
  
}
```

###Error density

####Estimation

```{r}
# Find mean bandwidth
bwFrame <- aggregate(error2 ~ neurType + diams + step, allData, function(x) density(x)$bw)

# Estimate density using mean bw per axonType
densList <- list()

for (i in c("Type1", "Type2", "Type3", "Type4")) {
  
  axData <- allData[allData$neurType == i, ]
  densList[[i]] <- by(axData$error2, axData$interact, function(x) density(x, 
                                                                           bw = ceiling(mean(bwFrame[bwFrame$neurType == i, "error2"]))))
  
}

# Restore format
densList <- unlist(densList, recursive = FALSE)
```


####Format and save

```{r}
# Extract from list
densXY <- lapply(densList, function(x) cbind(error = x$x, density = x$y))
densFrame <- data.frame(axonType = rep(c("Type1", "Type2", "Type3", "Type4"), each = 512*81),
                        step = rep(rep(unique(allData$step), each = 512*9), 4), 
                        step = rep(rep(unique(allData$diams), each = 512), 4*9))

# Store in long format data frame
densFrame <- cbind(densFrame, do.call(rbind, densXY))

# Save density as csv
write.csv(densFrame, "dataPepo/stereoSpheres_error_Densities.csv", sep=",", row.names = FALSE)

# Save original data as csv
oriFrame <- cbind(axonType = allData[, "neurType"], diams = allData[, "diams"], step = allData[, "step"],
                  error = allData[, "error2"], abserror = allData[, "error"])
write.csv(oriFrame, "dataPepo/stereoSpheres_error_OriginalData.csv", sep=",", row.names = FALSE)
```

###Cumul probability

####Estimate, format and save

```{r}
# Reformat as dataframe
binProb <- 0.5
probSeq <- seq(binProb, 100, binProb)
axProb <- readRDS("EsferasAnalysis/errorProbs_RMSE.rds")
probFrame <- data.frame(error = rep(probSeq, 81*4), 
                        neurType = rep(c("Type1", "Type2", "Type3", "Type4"), each = length(probSeq)*81),
                        diams = rep(rep(unique(allData$diams), each = length(probSeq)*9), 4),
                        step = rep(rep(unique(allData$step), each = length(probSeq)), 9*4),
                        prob = unlist(axProb))


n <- 200
grid <- expand.grid(diams = seq(10, 50, length = n), step = seq(70, 150, length = n))

# Smooth lm
for (h in c(5, 10)) {
  
  for (i in unique(probFrame$neurType)) {
    
    axlm <- lm(prob ~ diams * step, data = probFrame[probFrame$error == h & probFrame$neurType == i, ])
    axfit <- cbind(grid, prediction = predict(axlm, grid))
    
    # Save as csv
    write.csv(axfit, paste("dataPepo/stereoSpheres_cumulProb_error_", h, "_smooth_", i, ".csv", sep=""))
    
  }
  
  
}

```