Graph theory identifies altered prefrontal microcircuit organization in Shank3 mice, a mouse Model of autism.

1. ---
2. title: "K group fold ML winsorized w/grid tuning"
3. output: html_notebook
4. ---
6. This script is for Model 1, predicting genotype. Hyperparameter tuning will be employed for the mtry parameter.
8. ```{r, warning=FALSE, message=FALSE}
9. library(tidyverse)
10. library(caret)
11. library(randomForest)
12. library(beepr)
13. library(tictoc)
14. prl <- read.csv("winsorizedOutliersRemoved.csv", stringsAsFactors = T)
15. ```
17. # Pre-processing
19. The following features will be used: area H, area SB, area SN, normDiff_H_SB, normDiff_H_SN, Lg_SB, Cl_SB, Cl_SN, Cg_H, Cg_SB, Cg_SN.
21. ```{r}
22. # Select variables for ML.
23. mlVars <- prl %>%
24. select(-c(totalNeurons, Session))
25. ```
27. Convert ID variable to factor.
29. ```{r}
30. # Convert ID to factor.
31. mlVars$ID <- as.factor(mlVars$ID)
32. ```
34. We will be performing nested cross-validation. An outer loop will use leave two-out grouped cross-validation to determine the model error. An inner loop will use leave-one-out grouped CV for model training and validation.
36. # Functions for cross-validation
38. We want to create testing and validation sets that contain mice from both the WT and SHANK3 genotypes. To start, there are 11 WT mice and 9 SHANK3 mice. For the outer loop, we want to select 2 mice from each genotype. For the inner loop, we want to select 1 mouse from each genotype. For both loops, we want to iterate over all possible combinations of selecting mice for testing/validation from each genotype.
40. From a technical perspective, we want to create cv splits for tuning by index, similar to the output seen from caret's groupKfold(). This requires us to create a list of indices indicating the training set.
42. The following function will identify all possible combinations of testing/validation mice, based on how many mice will be selected from each genotype, and split the data for each combination into lists of indices indicating an observation's place in either the training or testing set.
44. ```{r}
45. # Create CV folds, sampling mice for validation.
46. cvFolds <- function(df, numMice){
48. # Subset WT mice and SHANK3 mice.
49. wt <- df %>% filter(gt == "WT") %>% pull(ID) %>% unique()
50. shank3 <- df %>% filter(gt == "SHANK3") %>% pull(ID) %>% unique()
52. # Define all possible combinations given the number of mice
53. # to select for each genotype.
54. # All ways to selecting WT mice.
55. wtTest <- combn(wt, numMice, simplify = FALSE)
56. # All ways to select SHANK3 mice.
57. shank3Test <- combn(shank3, numMice, simplify = FALSE)
58. # All ways to combine WT and SHANK3 mice.
59. # Output is a dataframe with two columns: WT & SHANK3.
60. validationMice <- expand.grid(wt = wtTest,
61. shank3 = shank3Test)
63. # Convert to list.
64. testMiceList <- lapply(1:nrow(validationMice),
65. function(x) c(unlist(validationMice[x, 1]),
66. unlist(validationMice[x, 2])))
68. # Get row index for sampling.
69. folds <- lapply(1:length(testMiceList), function(x) df %>%
70. rowid_to_column("index") %>%
71. filter(!ID %in% testMiceList[[x]]) %>%
72. pull(index))
74. # Name folds to match groupKfold()
75. names(folds) <- paste0("Fold", 1:length(testMiceList))
77. return(folds)
79. }
81. ```
83. ```{r}
84. # Start timer for tracking how long the code runs.
85. tic()
86. ```
88. ```{r}
89. # Create outer loop folds.
90. outerFolds <- cvFolds(df = mlVars, numMice = 2)
91. # Define number of outer fold iterations.
92. it <- length(outerFolds)
94. # Create list for storing results.
95. resultsList <- vector(mode = 'list', length = it)
96. # Create list for optimal variables.
97. optVariablesList <- vector(mode = 'list', length = it)
98. # Create list for storing confusion matrix from test results.
99. confMatList <- vector(mode = 'list', length = it)
100. # Create list for storing ROC curve information.
101. rocList <- vector(mode = 'list', length = it)
102. ```
103. Run nested CV procedure:
105. ```{r}
106. set.seed(33)
108. for(i in 1:it){
110. # Subset outer loop training data.
111. train <- mlVars[outerFolds[[i]], ]
112. # Subset outer loop testing data.
113. test <- mlVars[-outerFolds[[i]], ]
115. # Specify inner loop folds.
116. folds <- cvFolds(df = train, numMice = 1)
118. # Specify details for RFE.
119. ctrl <- rfeControl(functions = rfFuncs,
120. method = "cv",
121. index = folds,
122. number = length(folds),
123. saveDetails = T)
125. # Train.
126. rfFit <- rfe(gt ~ ., data = train %>% select(-c(ID)),
127. sizes = c(1:(ncol(train) - 2)), rfeControl = ctrl)
129. # Store results.
130. resultsList[[i]] <- cbind(rfFit$results, cv = i)
131. # Store optimal variables.
132. optVariablesList[[i]] <- rfFit$optVariables
134. # Subset training data to only include gt and optimal variables from rfe.
135. optTrain <- train %>% select(c(gt, rfFit$optVariables))
137. # Find max mtry.
138. maxMtry <- (ncol(optTrain) - 1) %>% sqrt() %>% ceiling()
140. # Re-fit model with optimal variables.
141. rfFitOpt <- train(gt ~ ., data = optTrain,
142. tuneGrid = data.frame(mtry = c(1:maxMtry)),
143. trControl = trainControl(index = folds,
144. method = "cv",
145. number = length(folds)))
147. # Predict test data.
148. preds <- predict(rfFitOpt,
149. newdata = test %>%
150. select(rfFit$optVariables))
152. # ROC information.
153. # Grab probabilities.
154. probs <- predict(rfFitOpt,
155. newdata = test %>%
156. select(rfFit$optVariables), type = "prob")
157. # Create data.frame of probabilities, true class, and iteration.
158. rocDat <- data.frame(trueClass = test$gt, probSHANK3 = probs$SHANK3,
159. cv = i)
160. # Add to rocList.
161. rocList[[i]] <- rocDat
163. # Performance table.
164. confusion <- confusionMatrix(preds, test$gt)
165. tableStats <- data.frame(tp = confusion$table[1, 1],
166. fn = confusion$table[2, 1],
167. fp = confusion$table[1, 2],
168. tn = confusion$table[2, 2])
170. # Store performance results.
171. confMatList[[i]] <- c(confusion$overall, confusion$byClass, tableStats)
173. print(i)
175. }
176. ```
178. ```{r}
179. # Audible cue letting user know cv is finished
180. beep('coin')
181. ```
183. ```{r}
184. # Stop timer and return time it took to run code.
185. toc()
186. ```
189. ```{r}
190. # Turn confusion matrix list into a data.frame.
191. confDF <- bind_rows(lapply(confMatList, as.data.frame.list))
192. # Turn results list into a data.frame.
193. resultsDF <- bind_rows(lapply(resultsList, as.data.frame.list))
194. # Turn optimal variables into a data.frame.
195. optVarsDF <- bind_rows(lapply(optVariablesList, as.data.frame.list))
196. # Turn rocList into a data.rame
197. rocDF <- bind_rows(lapply(rocList, as.data.frame.list))
198. ```
200. Save dataframes.
202. ```{r}
203. write.csv(confDF, "PrLwinsConfMatrix.csv", row.names = F)
204. write.csv(resultsDF, "PrLwinsResults.csv", row.names = F)
205. write.csv(optVarsDF, "PrLwinsOptVars.csv", row.names = F)
206. write.csv(rocDF, "PrLwinsROC.csv", row.names = F)
207. ```

PrLrfRFE.Rmd, under CC-BY-4.0 · at the source