OSCR

Structural Brain Correlates of Poor Reading Comprehension.

Code ↔ Paper

2 matches between paragraphs of the paper and lines of its authors' code, computed by the harvester (lexical-v1). Click a colored paragraph or line to see its counterpart.

The 2 matches
  1. [1] § MATERIALS AND METHODS › Group Classification › Cutoff classification model ↔ Code Files/Mahaffy_Masters_Real_DK__KM_20251223 (1).Rmd, lines 458–565 · score 0.66 · MatchIt, standard score, outlier, IQ, PC, age
  2. [2] § RESULTS › White Matter › Regression classification method ↔ Code Files/Mahaffy_Masters_Real_DK__KM_20251223 (1).Rmd, lines 1130–1237 · score 0.53 · left AF, right ILF, quadratic, linear, FA, residuals

Paper

Loaded from Europe PMC by your browser, not stored by OSCR: doi.org · Europe PMC

The paper is loaded when this pane is shown.

The authors' code

R Markdown · 1,237 lines · 68 KB · no license · 2 matches

  1. ---
  2. title: "Mahaffy_MA_Publication"
  3. author: "Kelly_Mahaffy"
  4. date: "`r Sys.Date()`"
  5. output: word_document
  6. ---
  8. # questions from Dan:
  9. # 1) Is it correct that no models are fit to the matched data?
  10. # 2) p-value adjustments: What should be grouped? (I don't have access to the folder with those outputs.)
  13. ```{r setup, include=FALSE}
  14. knitr::opts_chunk$set(echo = TRUE)
  15. library(psych)
  16. library(ggplot2)
  17. library(readxl)
  18. library(memisc)
  19. library(MatchIt)
  20. library(lm.beta)
  22. library(ggplot2)
  23. library(ggforce) # for geom_sina()
  24. library(gghalves)
  25. library(ggpubr)
  26. library(gridExtra) # for grid.arrange()
  28. # added by DK
  29. library(reshape2) # for melt()
  30. library(dplyr) # for rename()
  31. library(viridisLite) # for easier access to viridis color scale
  32. ```
  34. ```{r define dirs, include=FALSE}
  36. #super_dir <- "~/Mahaffy_PoorComprehenders_Pub"
  37. super_dir <- "/Users/dank/Documents/Haskins/Other people's projects/Kelly Mahaffy/MA_paper/Code_2024-12-17/"
  38. #
  39. # super_dir <- validate_dir(super_dir) # in case the final file separator character was omitted
  40. #
  41. data_dir <- paste0(super_dir, "Data_Files/")
  42. pvals_dir <- paste0(data_dir , "P_Values/" )
  44. dirs <- list()
  46. # file structure assumes this script is stored in a separate folder at the same level as Data_Files
  47. dirs[["code" ]] <- paste0(super_dir, "Code_Files/"); setwd(dirs[["code"]]) # set wd here!
  48. dirs[["data" ]] <- paste0(".." , .Platform$file.sep, "Data_Files", .Platform$file.sep)
  49. dirs[["input_brain" ]] <- paste0(dirs[["data"]], "Input_files" , .Platform$file.sep, "Brain_data", .Platform$file.sep)
  50. dirs[["input_pvals" ]] <- paste0(dirs[["data"]], "Input_files" , .Platform$file.sep, "P_values" , .Platform$file.sep)
  51. dirs[["input_other" ]] <- paste0(dirs[["data"]], "Input_files" , .Platform$file.sep, "Other" , .Platform$file.sep)
  52. dirs[["output_data" ]] <- paste0(dirs[["data"]], "Output_files", .Platform$file.sep, "Data" , .Platform$file.sep)
  53. dirs[["output_models"]] <- paste0(dirs[["data"]], "Output_files", .Platform$file.sep, "Models" , .Platform$file.sep)
  55. # create output directories that don't already exist
  56. for(next_output_dir in names(dirs)[grepl("^output_",names(dirs))]) {
  57. dir.create(dirs[[next_output_dir]], recursive=TRUE, showWarnings=FALSE)
  58. }
  60. output_file_prefix <- "Mahaffy_MA_"
  61. ```
  63. ```{r define functions, include=FALSE}
  65. # function which takes a filepath for a .csv, .xls or .xlsx file
  66. # and returns the data in that file
  67. read_data_file <- function(file_path, na="NA") {
  68. file_ext <- tail(strsplit(file_path, "\\.")[[1]], n=1)
  69. if(file_ext %in% c('xls', 'xlsx')) {
  70. file_data <- read_excel(file_path, na=na)
  71. } else if(file_ext %in% 'csv') {
  72. file_data <- read.csv(file_path, na=na)
  73. } else {
  74. error(paste0(file_path, ' file extension must be .csv, .xls, or .xlsx!'))
  75. }
  76. return(file_data)
  77. }
  79. # function which takes a list of data frames (data_list), a data frame to merge with the others (merge_with_others),
  80. # and (optionally) a column by which to join them (uses "id" by default).
  81. # returns a list with all data frames in the list merged with the to-be-merged data frame.
  82. merge_one_with_all <- function(data_list, merge_with_others, join_by="id") {
  83. merged_list <- sapply(data_list, function(next_file_data) {merge(next_file_data, merge_with_others, by=join_by)})
  84. return(merged_list)
  85. }
  87. # function which takes a list of data frames (data_list), the name of a data frame in that list to merge with the others (merge_with_others),
  88. # and (optionally) a column by which to join them (uses "id" by default).
  89. # returns a list with all data frames except the named one merged with the named data frame.
  90. # (builds on merge_one_with_all() by breaking out the to-be-merged data frame from the list, then passing them as separate args.)
  91. merge_one_with_others <- function(data_list, merge_with_others, join_by="id") {
  92. names_to_merge <- setdiff(names(data_list), merge_with_others)
  93. merged_list <- merge_one_with_all(data_list[names_to_merge], data_list[[merge_with_others]], join_by)
  94. return(merged_list)
  95. }
  97. # function which takes a list of data frames (e.g., brain_data[["grey"]]), a column to center (e.g., "cnr"), and the name of an output centered column ("cnr_c").
  98. # returns the list with a centered column in each data frame
  99. center_col_by_matter_type <- function(df_list, col_to_center, centered_col_name) {
  100. for(next_file_data in names(df_list)) {
  101. df_list[[next_file_data]][[centered_col_name]] <- scale(df_list[[next_file_data]][[col_to_center]], scale=F)
  102. }
  103. return(df_list)
  104. }
  106. ##KM attempt at scaling
  107. scale_col_by_matter_type <- function(df_list, col_to_scale, scaled_col_name) {
  108. for(next_file_data in names(df_list)) {
  109. df_list[[next_file_data]][[scaled_col_name]] <- scale(df_list[[next_file_data]][[col_to_scale]])
  110. }
  111. return(df_list)
  112. }
  114. # function that takes a list of data frames, a list of data frame names to write,
  115. # and a directory to write them in, and then writes those data frames as CSVs
  116. # in that dir. by default, output directory is dirs[["output_data"]]
  117. write_df_as_csv <- function(df_list, dfs_to_write, write_dir=dirs[["output_data"]]) {
  118. for(next_df_to_write in dfs_to_write) {
  119. write.csv(df_list[[next_df_to_write]], paste0(write_dir, next_df_to_write, ".csv"))
  120. }
  121. }
  123. # function that takes a data frame as input and a dependent variable (as a character),
  124. # then calls matchit() 2x and creates some summary output (which may not be output
  125. # in a usable way, though this could be easily changed). it returns a data frame
  126. # with the appropriate columns for matching.
  127. custom_matchit_calls <- function(df, dv) {
  128. matchit_formula <- formula(paste0(dv, " ~ age + sex + WISC_NVI_ss"))
  130. summary(m.out.1 <- matchit(matchit_formula, data=df, method=NULL , distance="glm" ))
  132. m.out.2 <- matchit(matchit_formula, data=df, method="optimal", distance="mahalanobis", replace=FALSE)
  134. summary(m.out.2)
  135. plot(m.out.2, type="qq")
  136. return(match.data(m.out.2))
  137. }
  139. # function to combine a PCPD and a PCTD dataset rather than doing this outside R
  140. merge_pcpd_pctd <- function(pcpd, pctd) {
  141. # not an rbind, because a given participant may appear in both, and we only want them to have one row.
  142. # for a given sub who appears in both, all values are identical (though the _Match col is called PC_PD_Match in one and PC_TD_Match in the other)
  144. # rename cols
  145. colnames(pctd)[colnames(pctd)=="PC_TD_Match"] <- "PC_TD"
  146. colnames(pcpd)[colnames(pcpd)=="PC_PD_Match"] <- "PC_PD"
  148. # create a df with one sub per row if they appear in at least one of the two matched data sets
  149. # remove X col, if present -- would need to be treated differently because vals would be different for PCPD and PCTD, and I don't think it's used later anyway
  151. merged_data <- unique.data.frame(rbind(
  152. as.data.frame(pctd[,setdiff(colnames(pctd), c("PC_TD", "X"))]),
  153. as.data.frame(pcpd[,setdiff(colnames(pcpd), c("PC_PD", "X"))])
  154. ))
  156. # merge PC_PD and PC_TD col headers
  157. merged_headers <- merge(pctd[,c("id","PC_TD")], pcpd[,c("id","PC_PD")], all.x=TRUE, all.y=TRUE)
  158. # add PC_All col headers: 1 if both !is.na(PC_TD) and !is.na(PC_PD), and 0 otherwise
  159. merged_headers[["PC_All"]] <- 1-apply(sapply(merged_headers[,c("PC_TD","PC_PD")], is.na), 1, sum)
  161. # merge headers and data
  162. merged_data <- merge(merged_headers, merged_data)
  164. # sort: PC_All (descending), PC_TD (descending), PC_PD (descending), id (ascending)
  165. merged_data <- merged_data[with(merged_data, order(-PC_All, -PC_TD, -PC_PD, id)),]
  167. # to do (maybe): add X col
  168. # NOTE: There is an error in the way this is handled in the previous code --
  169. # X is nested within matched dataset, so a given participant should have
  170. # either 1 or 2 values (potentially one each for PC_PD and PC_TD).
  171. # However, X is never used again, so it doesn't really matter. Safe to
  172. # just omit it for now.
  174. return(merged_data)
  175. }
  177. # function which takes various pieces of information about a set of lm() models to fit,
  178. # then writes the appropriate files with results and p-values, and returns
  179. # an object containing the fitted models and p-values
  180. # df = the data frame to be used in all models
  181. # dvs = a character vector of DVs (1 per model)
  182. # pred_pval = the string name of the variable for which to save the p-value
  183. # preds_other = a character vector of other predictors (not pred_pval) to include in the model
  184. # id_text = an identifying string to be prepended to the names of files that are written
  186. # df = get(next_model_df)[[paste0(output_file_prefix, next_model_info[["df"]])]]
  187. # dvs = dvs[[models_to_fit[next_model_idx,"dvs"]]]
  188. # pred_pval = next_model_info[["pred_pval"]]
  189. # preds_other = strsplit(next_model_info[["preds_other"]], " \\+ ")[[1]]
  190. # id_text = next_model_info[["id_text"]]
  192. fit_lms <- function(df, dvs, pred_pval, preds_other, id_text) {
  194. all_preds <- setdiff(c(pred_pval, preds_other), "") # exclude empty preds
  196. # convert data from wide to long
  197. df.long <- melt(df[,c(dvs,all_preds)], id.vars=all_preds, variable.name="region", value.name="dv")
  199. # construct formula -- will be the same for all models
  200. lm_formula_baseline <- formula(paste0("dv ~ ", paste(setdiff(preds_other, ""), collapse=" + ")))
  201. lm_formula <- formula(paste0("dv ~ ", paste( all_preds , collapse=" + ")))
  203. # fit all models
  204. all_baseline_lms <- lapply(dvs, function(x) {next_lm <- lm.beta(lm(lm_formula_baseline, df.long[df.long[["region"]]==x,]))})
  205. names(all_baseline_lms) <- dvs
  207. all_lms <- lapply(dvs, function(x) {next_lm <- lm.beta(lm(lm_formula, df.long[df.long[["region"]]==x,]))})
  208. names(all_lms) <- dvs
  210. # get r^2 vals
  211. rsq_vals <- do.call(rbind.data.frame, lapply(dvs, function(x) {
  212. rsq <- c(summary(all_lms[[x]])[["r.squared"]], summary(all_baseline_lms[[x]])[["r.squared"]]);
  213. rsq <- c(rsq, rsq[1]-rsq[2])}))
  214. rsq.headers <- c("rsq.linear", "rsq.intercept", "rsq.linear.partial")
  215. colnames(rsq_vals) <- rsq.headers
  217. #names(all_lms) <- paste0("lm", 1:length(all_lms)) # necessary for mtable() call to work
  219. # generate and save mtable
  220. # note: only works if the models are not standardized (due to extra column in model summary)
  221. # lm_mtable <- do.call(mtable, all_lms)
  222. # write.mtable(lm_mtable, file=paste0(dirs[["output_models"]], output_file_prefix, id_text, "_Results.csv"), format=c("delim"))
  224. # make our own output
  225. all_lms_for_output <- do.call(rbind.data.frame, lapply(1:length(all_lms), function(x) {data.frame(model_num=x, dv=dvs[[x]], effect=rownames(coef(summary(all_lms[[x]]))), coef(summary(all_lms[[x]])), rsq_vals[x,])}))
  226. rownames(all_lms_for_output) <- NULL
  227. colnames(all_lms_for_output) <- c("model_num","dv","effect","Estimate","Standardized","SE","t","p",rsq.headers)
  228. write.table(all_lms_for_output, file=paste0(dirs[["output_models"]], output_file_prefix, id_text, "_Results.csv"), sep=",")
  230. # extract and save p-values from specified predictor
  231. lm_pvals <- sapply(1:length(all_lms), function(next_lm_idx) {coef(summary(all_lms[[next_lm_idx]]))[pred_pval,"Pr(>|t|)"]})
  232. write.csv(lm_pvals, paste0(dirs[["output_models"]], output_file_prefix, id_text, "_PValue.csv"))
  235. objects_to_return <- list(model=all_lms, rsq.linear=data.frame(dv=dvs, rsq_vals), p=lm_pvals)
  237. # quad formula if needed
  238. unique_pred_vals <- unique(df.long[[pred_pval]])
  239. unique_pred_vals <- unique_pred_vals[!is.na(unique_pred_vals)] # exclude NAs
  240. pred_is_continuous <- length(unique_pred_vals) > 2
  241. if(pred_is_continuous) {
  242. pred_linear_and_quad <- as.data.frame(poly(df.long[[pred_pval]], 2))
  243. colnames(pred_linear_and_quad) <- paste0(pred_pval, c(".lin", ".quad"))
  244. df.long <- data.frame(df.long, pred_linear_and_quad)
  246. # construct quadratic formula
  247. lm_formula.quad <- formula(paste0("dv ~ ", paste(c(all_preds, paste0(pred_pval, ".quad")), collapse=" + ")))
  249. all_lms.quad <- lapply(dvs, function(x) {next_lm <- lm(lm_formula.quad, df.long[df.long[["region"]]==x,])})
  250. names(all_lms.quad) <- dvs
  252. rsq_vals.quad <- do.call(rbind.data.frame, lapply(1:length(all_lms), function(x) {
  253. rsq <- summary(all_lms.quad[[x]])[["r.squared"]]
  254. rsq <- c(rsq, rsq - rsq_vals[x,"rsq.linear"])
  255. }))
  256. rsq.quad.headers <- c("rsq.quad", "rsq.quad.partial")
  257. colnames(rsq_vals.quad) <- rsq.quad.headers
  259. # compare p-values for each one
  260. all_quad_stats <- do.call(rbind.data.frame, lapply(1:length(all_lms), function(x) {
  261. next_quad_comparison <- anova(all_lms[[x]], all_lms.quad[[x]])
  262. next_quad_comparison.stats <- data.frame(with(next_quad_comparison, data.frame(F=F, df=Df, p=`Pr(>F)`)[2,]), rsq_vals.quad[x,])
  263. return(next_quad_comparison.stats)
  264. }))
  266. objects_to_return[["quad_models" ]] <- all_lms.quad
  267. objects_to_return[["quad_results"]] <- data.frame(dv=dvs, all_quad_stats)
  269. }
  271. return(objects_to_return)
  272. }
  274. # function which takes a data frame containing brain data (df) and a data frame
  275. # (var_info) containing two cols: "match" and "store". for each row in var_info,
  276. # it averages across all values in columns named "slf[digit(s)]_x" (where x is
  277. # the value of "match") and stores the results in "slfavg_y" (where y is the
  278. # value of "store").
  280. create_slf_avg_vars <- function(df, var_info) {
  282. # average over appropriately named cols
  283. new_cols <- sapply(var_info[["match"]], function(next_string) {
  284. rowMeans(df[,grepl(paste0("^slf\\d+_", next_string, "$"), colnames(df))], na.rm=TRUE)
  285. })
  287. # add new cols to df and return
  288. df[,paste0("slfavg_", var_info[["store"]])] <- new_cols
  289. return(df)
  290. }
  293. create_slf_avg_vars_fake <- function(df, var_info) { # for printing only
  295. # average over appropriately named cols
  296. summary_strings <- sapply(1:nrow(var_info), function(next_row_idx) {
  297. next_match_string <- var_info[next_row_idx,"match"]
  298. cols_to_combine <- colnames(df)[grepl(paste0("^slf\\d+_", next_match_string, "$"), colnames(df))]
  299. next_store_string <- paste0("slfavg_", var_info[next_row_idx,"store"], " = mean(", paste(cols_to_combine, collapse=", "), ")")
  300. return(next_store_string)
  301. })
  303. return(summary_strings)
  304. }
  307. get_model_resids <- function(fitted_lm, pred_var) {
  308. outputs <- list()
  310. outputs[["fitted_lm.formula_split"]] <- as.character(formula(fitted_lm))[c(2,3)]
  311. outputs[["fitted_lm.dv"]] <- outputs[["fitted_lm.formula_split"]][1]
  312. outputs[["fitted_lm.preds"]] <- strsplit(outputs[["fitted_lm.formula_split"]][2], ' \\+ ')[[1]]
  314. outputs[["subset_lm.preds"]] <- setdiff(outputs[["fitted_lm.preds"]], pred_var)
  315. outputs[["subset_lm.formula"]] <- formula(paste0(outputs[["fitted_lm.dv"]], " ~ ", paste(outputs[["subset_lm.preds"]], collapse=" + ")))
  316. outputs[["subset_lm"]] <- lm(outputs[["subset_lm.formula"]], data=fitted_lm[["model"]])
  317. outputs[["subset_lm.resids"]] <- resid(outputs[["subset_lm"]])
  319. outputs[["pred_lm.formula"]] <- formula(paste0(pred_var, " ~ ", paste(outputs[["subset_lm.preds"]], collapse=" + ")))
  320. outputs[["pred_lm"]] <- lm(outputs[["pred_lm.formula"]], data=fitted_lm[["model"]])
  321. outputs[["pred.resids"]] <- resid(outputs[["pred_lm"]])
  323. return(outputs)
  324. }
  326. plot_grouped_vs_continuous <- function(fitted_lm, pred_var, num_quantiles=5, highlighted_quantiles=c(1,3,5), degree=c(1), linetype=NULL) {
  328. resid_models <- get_model_resids(fitted_lm, pred_var)
  330. # get plot information, colors, quantiles, etc.
  331. plot_data <- fitted_lm[["model"]]
  333. plot.y.label <- paste0(resid_models[["fitted_lm.dv"]], ".resid")
  334. plot_data[[plot.y.label]] <- resid_models[["subset_lm.resids"]]
  336. plot.x.label <- paste0(pred_var, ".resid")
  337. plot_data[[plot.x.label]] <- resid_models[["pred.resids"]]
  339. plot.quantiles <- quantile(plot_data[[plot.x.label]], probs=(1:(num_quantiles-1))/num_quantiles)
  340. plot_data[["quantile"]] <- as.factor(as.character(sapply(plot_data[[plot.x.label]], function(x) {1 + sum(x > plot.quantiles)})))
  342. # compute quantile boundaries as mean of mins/maxes of each quantile
  343. quantile_boundaries <- lapply(c(min,max), function(x) {aggregate(formula(paste0(plot.x.label, " ~ quantile")), data=plot_data, FUN=x)})
  344. quantile_boundaries <- data.frame(max=quantile_boundaries[[2]][1:(num_quantiles-1),plot.x.label], min=quantile_boundaries[[1]][2:num_quantiles,plot.x.label])
  345. quantile_boundaries[["mean"]] <- apply(quantile_boundaries, 1, mean)
  347. # create color map: grey for non-highlighted quantiles, viridis for highlighted quantiles
  348. plot_colors <- gray.colors(num_quantiles)
  349. plot_colors[highlighted_quantiles] <- viridis(num_quantiles)[highlighted_quantiles]
  351. # get means by quantile (there is SO a better way to do this...)
  352. summary_stats <- do.call(rbind.data.frame, lapply(c(plot.x.label, plot.y.label), function(y) {
  353. next_summary_stats <- do.call(cbind, lapply(c(mean,sd,length), function(x) {aggregate(formula(paste0(y, " ~ quantile")), data=plot_data, FUN=x)[,2]}))
  354. colnames(next_summary_stats) <- c("mean","sd","n")
  355. #next_summary_stats <- data.frame(var=rep(c("x","y"),each=num_quantiles), next_summary_stats)
  356. next_summary_stats <- data.frame(quantile=1:num_quantiles, next_summary_stats)
  357. next_summary_stats
  358. }))
  359. summary_stats <- data.frame(var=rep(c("x","y"), each=num_quantiles), summary_stats)
  360. summary_stats[["se" ]] <- with(summary_stats, sd/sqrt(n))
  361. summary_stats[["start"]] <- with(summary_stats, mean-se)
  362. summary_stats[["end" ]] <- with(summary_stats, mean+se)
  364. # oy
  365. summary_stats <- do.call(merge, lapply(c("x","y"), function(x) {
  366. subset_summary_stats <- summary_stats[summary_stats[["var"]]==x,]
  367. subset_col_names <- setdiff(colnames(subset_summary_stats), c("var","quantile"))
  368. colnames(subset_summary_stats)[colnames(subset_summary_stats) %in% subset_col_names] <- paste0(subset_col_names, ".", x)
  369. subset(subset_summary_stats, select=-var)
  370. }))
  373. # set figure x-limits at [range] + 10% on each end
  374. x.range <- range(plot_data[[plot.x.label]])
  375. x.range <- x.range + (diff(x.range) * .10) * c(-1,1)
  377. # make plot
  378. resid_plot <- ggplot(plot_data, aes_string(x=plot.x.label, y=plot.y.label, color="quantile")) +
  379. geom_point() +
  380. scale_color_manual(values=plot_colors) +
  381. theme_bw() +
  382. scale_x_continuous(limits=x.range)
  384. for(next_degree_idx in 1:length(degree)) {
  385. next_degree <- degree[next_degree_idx]
  387. if(next_degree > 1) {
  388. stat_smooth.formula <- formula(paste0("y ~ x + ", paste(paste0("I(x^", 2:next_degree, ")"), collapse=" + ")))
  389. } else {
  390. stat_smooth.formula <- formula("y ~ x")
  391. }
  393. if(is.null(linetype) || length(linetype) != length(degree)) {
  394. next_linetype <- next_degree
  395. } else {
  396. next_linetype <- linetype[next_degree_idx]
  397. }
  398. resid_plot <- resid_plot + stat_smooth(aes(group=1), color="black", linetype=next_linetype, method="lm", formula=stat_smooth.formula)
  399. }
  401. # add quantile boundaries
  402. for(next_quantile_boundary_idx in 1:nrow(quantile_boundaries)) {
  403. resid_plot <- resid_plot + geom_vline(xintercept=quantile_boundaries[next_quantile_boundary_idx,"mean"], color="black", linetype="dashed", alpha=0.2)
  404. }
  406. # just add plus signs rather than bars with +/- 1 SEM
  407. resid_plot <- resid_plot +
  408. geom_point(data=summary_stats, aes(x=mean.x, y=mean.y), color="black", size=5, shape=3)
  410. # for(next_quantile_mean_idx in (1:num_quantiles)) {
  411. # resid_plot <- resid_plot +
  412. # geom_segment(color=plot_colors[next_quantile_mean_idx], linewidth=1.5, # shows variability in X
  413. # x =summary_stats[next_quantile_mean_idx,"start.x"],
  414. # xend=summary_stats[next_quantile_mean_idx,"end.x" ],
  415. # y =summary_stats[next_quantile_mean_idx,"mean.y" ],
  416. # yend=summary_stats[next_quantile_mean_idx,"mean.y" ]) +
  417. # geom_segment(color=plot_colors[next_quantile_mean_idx], linewidth=1.5, # shows variability in Y
  418. # x =summary_stats[next_quantile_mean_idx,"mean.x" ],
  419. # xend=summary_stats[next_quantile_mean_idx,"mean.x" ],
  420. # y =summary_stats[next_quantile_mean_idx,"start.y"],
  421. # yend=summary_stats[next_quantile_mean_idx,"end.y" ])
  422. # }
  424. return(resid_plot)
  425. }
  427. ```
  429. ## Describe
  430. ```{r}
  431. Mahaffy_MA_Final <- read_data_file(paste0(dirs[["input_other"]], "Mahaffy_Masters_Final.xlsx"),na = "NA")
  433. #objects(Mahaffy_MA_Final)
  434. Mahaffy_MA_Final$age_round<-round(Mahaffy_MA_Final$age,1)
  435. describe(Mahaffy_MA_Final$age_round)
  436. sum(Mahaffy_MA_Final$sex)
  437. describe(Mahaffy_MA_Final$wiat_rc_stand)
  438. describe(Mahaffy_MA_Final$wiat_rc_raw)
  439. describe(Mahaffy_MA_Final$wiat_lcrv_ss)
  440. describe(Mahaffy_MA_Final$wiat_lcrv_raw)
  441. describe(Mahaffy_MA_Final$towre_pde_ss)
  442. describe(Mahaffy_MA_Final$towre_pde_raw)
  443. describe(Mahaffy_MA_Final$WISC_NVI_ss)
  444. ```
  446. ## Regression Model
  447. ```{r}
  448. BehavioralRegression<-lm(wiat_rc_raw ~ wiat_lcrv_raw + towre_pde_raw + age_round + sex + WISC_NVI_raw, data=Mahaffy_MA_Final)
  449. summary(BehavioralRegression)
  450. BehavioralRegression_Interactions<-lm(wiat_rc_raw ~ (wiat_lcrv_raw + towre_pde_raw + age_round + sex + WISC_NVI_raw)^2 , data=Mahaffy_MA_Final)
  451. summary(BehavioralRegression_Interactions)
  452. anova(BehavioralRegression, BehavioralRegression_Interactions)
  453. Mahaffy_MA_Final$Residuals_BehavioralRegression<-resid(BehavioralRegression_Interactions)
  454. plot(Mahaffy_MA_Final$Residuals_BehavioralRegression)
  455. hist(Mahaffy_MA_Final$Residuals_BehavioralRegression)
  456. ```
  458. ##Preparing Data
  459. ```{r}
  460. Mahaffy_MA_Final$age_c<-scale(Mahaffy_MA_Final$age_round, scale=FALSE)
  461. Mahaffy_MA_Final$word_c<-scale(Mahaffy_MA_Final$towre_pde_raw, scale=FALSE)
  462. Mahaffy_MA_Final$vocab_c<-scale(Mahaffy_MA_Final$wiat_lcrv_raw, scale=FALSE)
  463. Mahaffy_MA_Final$iq_c<-scale(Mahaffy_MA_Final$WISC_NVI_raw, scale=FALSE)
  467. #Import Brain Data
  469. brain_data <- list()
  471. #Grey Matter
  472. # grey matter files
  473. files.grey <- c("FS7_compiled_lh_area_2009", "FS7_compiled_lh_thickness_2009", "FS7_compiled_rh_area_2009", "FS7_compiled_rh_thicnkess_2009", "FS7_Compiled_Vol_Seg", "T1w_QC_Full")
  474. # read in files to list
  475. brain_data[["grey"]] <- lapply(files.grey, function(next_file) {read_data_file(paste0(dirs[["input_brain"]], next_file, ".xlsx"), na = "NA")})
  476. names(brain_data[["grey"]]) <- files.grey
  479. #White Matter
  480. # white matter files
  481. files.white <- c("WM_FA_MD_Tract_All", "WM_NODDI_All", "WM_QC_All")
  482. # read in files to list
  483. brain_data[["white"]] <- lapply(files.white, function(next_file) {read_data_file(paste0(dirs[["input_brain"]], next_file, ".xlsx"), na = "0")})
  484. names(brain_data[["white"]]) <- files.white
  487. #Combining Data
  489. brain_data[["grey" ]] <- merge_one_with_others(brain_data[["grey" ]], "T1w_QC_Full")
  490. brain_data[["white"]] <- merge_one_with_others(brain_data[["white"]], "WM_QC_All" )
  492. # rename dfs
  493. df_names <- c(
  494. LH_Area = "FS7_compiled_lh_area_2009",
  495. LH_Thickness = "FS7_compiled_lh_thickness_2009",
  496. RH_Area = "FS7_compiled_rh_area_2009",
  497. RH_Thickness = "FS7_compiled_rh_thicnkess_2009",
  498. Vol = "FS7_Compiled_Vol_Seg",
  499. Noddi = "WM_NODDI_All",
  500. MDFA = "WM_FA_MD_Tract_All"
  501. )
  502. names(brain_data[["grey" ]]) <- paste0(output_file_prefix, names(df_names)[match(names(brain_data[["grey" ]]), unname(df_names))])
  503. names(brain_data[["white"]]) <- paste0(output_file_prefix, names(df_names)[match(names(brain_data[["white"]]), unname(df_names))])
  506. #Centering Data
  508. # center CNR
  509. brain_data[["grey" ]] <- center_col_by_matter_type(brain_data[["grey" ]], "cnr" , "cnr_c")
  510. brain_data[["white"]] <- center_col_by_matter_type(brain_data[["white"]], "Avg _CNR", "cnr_c")
  513. #Scaling Data
  514. brain_data[["grey" ]] <- scale_col_by_matter_type(brain_data[["grey" ]], "cnr" , "cnr_c")
  515. brain_data[["white"]] <- scale_col_by_matter_type(brain_data[["white"]], "Avg _CNR", "cnr_c")
  517. #Combining Brain and Behavioral Data
  519. brain_data[["grey" ]] <- merge_one_with_all(brain_data[["grey" ]], Mahaffy_MA_Final)
  520. brain_data[["white"]] <- merge_one_with_all(brain_data[["white"]], Mahaffy_MA_Final)
  523. #Making spreadsheets for matching
  525. #Only using one GM metric here because I know the participants are the same across all GM data
  527. brain_data[["grey"]][["Mahaffy_GM_PC" ]] <- subset(brain_data[["grey"]][[paste0(output_file_prefix, "LH_Area")]], wiat_rc_stand <= 90 & towre_pde_ss >= 100)
  528. brain_data[["grey"]][["Mahaffy_GM_Possible_PD"]] <- subset(brain_data[["grey"]][[paste0(output_file_prefix, "LH_Area")]], wiat_rc_stand >= 100 & towre_pde_ss <= 90 & wiat_lcrv_ss >= 100)
  529. brain_data[["grey"]][["Mahaffy_GM_Possible_TD"]] <- subset(brain_data[["grey"]][[paste0(output_file_prefix, "LH_Area")]], wiat_rc_stand >= 100 & towre_pde_ss >= 100 & wiat_lcrv_ss >= 100)
  531. write_df_as_csv(brain_data[["grey"]], c("Mahaffy_GM_PC", "Mahaffy_GM_Possible_PD", "Mahaffy_GM_Possible_TD"))
  534. #Outside of R, I checked for and removed duplicate values, combined all there of the above into a single spreadsheet, and dummy coded both PC_PD and PC_TD variables, with poor comprehenders always being 1. Spreadsheet re-uploaded as Mahaffy_MatchGroups. We also identified one PC who was an outlier (a RC standard score more that 10 points from the next lowest SS) and created a second spreadsheet for group matching without this outlier. Spreadsheets were split into a PC_PD and PC_TD spreadsheet (with all poor comprehenders and possible PDs and TDs, respectively, for compliance with MatchIt! package used below.)
  536. # matched data files
  537. match_pool_files.grey <- c("Mahaffy_MatchGroups_PCPD", "Mahaffy_MatchGroups_PCTD", "Mahaffy_MatchGroups_NoOut_PCPD", "Mahaffy_MatchGroups_NoOut_PCTD")
  539. # read in files to list
  540. match_pool_data <- list()
  541. match_pool_data[["grey"]] <- lapply(match_pool_files.grey, function(next_file) {read_data_file(paste0(dirs[["input_other"]], next_file, ".xlsx"), na = "NA")})
  542. names(match_pool_data[["grey"]]) <- match_pool_files.grey
  545. #Repeating the Process for White Matter
  547. #Here, doing this process once for WM metrics from FSL and once for NODDI because the participants are not all the same
  549. brain_data[["white"]][["Mahaffy_WM_PC" ]] <- subset(brain_data[["white"]][[paste0(output_file_prefix, "MDFA" )]], wiat_rc_stand <= 90 & towre_pde_ss >= 100)
  550. brain_data[["white"]][["Mahaffy_WM_Possible_PD" ]] <- subset(brain_data[["white"]][[paste0(output_file_prefix, "MDFA" )]], wiat_rc_stand >= 100 & towre_pde_ss <= 90 & wiat_lcrv_ss >= 100)
  551. brain_data[["white"]][["Mahaffy_WM_Possible_TD" ]] <- subset(brain_data[["white"]][[paste0(output_file_prefix, "MDFA" )]], wiat_rc_stand >= 100 & towre_pde_ss >= 100 & wiat_lcrv_ss >= 100)
  553. brain_data[["white"]][["Mahaffy_Noddi_PC" ]] <- subset(brain_data[["white"]][[paste0(output_file_prefix, "Noddi")]], wiat_rc_stand <= 90 & towre_pde_ss >= 100)
  554. brain_data[["white"]][["Mahaffy_Noddi_Possible_PD"]] <- subset(brain_data[["white"]][[paste0(output_file_prefix, "Noddi")]], wiat_rc_stand >= 100 & towre_pde_ss <= 90 & wiat_lcrv_ss >= 100)
  555. brain_data[["white"]][["Mahaffy_Noddi_Possible_TD"]] <- subset(brain_data[["white"]][[paste0(output_file_prefix, "Noddi")]], wiat_rc_stand >= 100 & towre_pde_ss >= 100 & wiat_lcrv_ss >= 100)
  557. write_df_as_csv(brain_data[["white"]], c("Mahaffy_WM_PC" , "Mahaffy_WM_Possible_PD" , "Mahaffy_WM_Possible_TD" ,
  558. "Mahaffy_Noddi_PC", "Mahaffy_Noddi_Possible_PD", "Mahaffy_Noddi_Possible_TD"))
  561. match_pool_files.white <- c("Mahaffy_WM_MatchGroups_PCPD", "Mahaffy_WM_MatchGroups_PCTD", "Mahaffy_Noddi_MatchGroups_PCPD", "Mahaffy_Noddi_MatchGroups_PCTD")
  563. match_pool_data[["white"]] <- lapply(match_pool_files.white, function(next_file) {read_data_file(paste0(dirs[["input_other"]], next_file, ".csv"), na = "NA")})
  564. names(match_pool_data[["white"]]) <- match_pool_files.white
  565. ```
  568. ##Matching For Cutoff
  569. ```{r}
  571. # create a data frame with all parameters needed, then call custom_matchit_calls() via an apply function to match everything at once
  572. data_to_match <- data.frame(rbind(
  573. c("grey" , "Mahaffy_MatchGroups_PCPD" , "PC_PD_Match", "Matched_PCPD_GM" )
  574. , c("grey" , "Mahaffy_MatchGroups_PCTD" , "PC_TD_Match", "Matched_PCTD_GM" )
  575. , c("grey" , "Mahaffy_MatchGroups_NoOut_PCPD", "PC_PD_Match", "Matched_PCPD_NoOut")
  576. , c("grey" , "Mahaffy_MatchGroups_NoOut_PCTD", "PC_TD_Match", "Matched_PCTD_NoOut")
  577. , c("white", "Mahaffy_WM_MatchGroups_PCPD" , "PC_PD_Match", "Matched_PCPD_WM" )
  578. , c("white", "Mahaffy_WM_MatchGroups_PCTD" , "PC_TD_Match", "Matched_PCTD_WM" )
  579. , c("white", "Mahaffy_Noddi_MatchGroups_PCPD", "PC_PD_Match", "Matched_PCPD_Noddi")
  580. , c("white", "Mahaffy_Noddi_MatchGroups_PCTD", "PC_TD_Match", "Matched_PCTD_Noddi")
  581. ))
  582. names(data_to_match) <- c("matter_type", "df_name", "dv", "id_text")
  584. # run all models in data_to_match
  585. matched_group_data <- lapply(1:nrow(data_to_match), function(next_data_to_match_idx) {
  586. next_match_info <- data_to_match[next_data_to_match_idx,]
  587. next_matched_data <- custom_matchit_calls(match_pool_data[[next_match_info[["matter_type"]]]][[next_match_info[["df_name"]]]], next_match_info[["dv"]])
  588. return(next_matched_data)
  589. })
  590. names(matched_group_data) <- paste0(output_file_prefix, data_to_match[["id_text"]])
  592. # write CSVs
  593. write_df_as_csv(matched_group_data, data_to_match[["id_text"]])
  597. #Outside of R, I combined each of these into a single spreadsheet with all contrasts of interest, which resulted in three spreadsheets- one for GM, one for WM, and one for Noddi
  599. ## we can do that in R!
  601. # define the four file groupings
  602. matched_group_data.groupings <- c("GM","WM","Noddi","NoOut")
  603. # for each grouping, combine the two dfs
  604. matched_group_data.merged <- lapply(matched_group_data.groupings, function(next_data) {
  605. # algorithmically construct the names of the fields to pull
  606. ordered_pcpd_pctd <- c("PCPD","PCTD")
  607. matched_data_names <- paste0(output_file_prefix, "Matched_", ordered_pcpd_pctd, "_", next_data)
  608. names(matched_data_names) <- ordered_pcpd_pctd
  610. # run custom function which combines the data
  611. merged_data <- merge_pcpd_pctd(pcpd=matched_group_data[[matched_data_names["PCPD"]]], pctd=matched_group_data[[matched_data_names["PCTD"]]])
  612. return(merged_data)
  613. })
  614. # name the fields in the resulting list
  615. names(matched_group_data.merged) <- paste0(output_file_prefix, "Matched_", matched_group_data.groupings)
  618. #Writing a CSV to prepare groups for mixed regression (save the lowest, middle, and highest 20% of residual scores with the lowest group being dummy coded as 1 and all others as 0) and create a list to merge with other spreadsheets like for previous. Done separately for GM, WM, and Noddi.
  620. write_df_as_csv(c(brain_data[["grey"]], brain_data[["white"]]), paste0(output_file_prefix, c("LH_Area", "MDFA", "Noddi"))) # mix of GM & WM
  623. # reading these in from the output directory (?) because they were just created
  624. brain_data[["grey" ]][[paste0(output_file_prefix, "LH_Area")]] <- read_data_file(paste0(dirs[["output_data"]], output_file_prefix, "LH_Area.csv"), na = "NA")
  625. brain_data[["white"]][[paste0(output_file_prefix, "MDFA" )]] <- read_data_file(paste0(dirs[["output_data"]], output_file_prefix, "MDFA.csv" ), na = "NA")
  626. brain_data[["white"]][[paste0(output_file_prefix, "Noddi" )]] <- read_data_file(paste0(dirs[["output_data"]], output_file_prefix, "Noddi.csv" ), na = "NA")
  629. #Reading in lists to merge so that group contrasts can be merged with other data files
  631. merged_data_files <- data.frame(rbind( # keep all info about each file together
  632. c("GM_Matched_ToMerge.xlsx" , "GM_Cutoff_Matched_Merge", "_Cutoff" , 1)
  633. , c("Mixed_Merge.csv" , "Mixed_Merge" , "" , 2)
  634. , c("Matched_NoOut_ToMerge.csv", "Matched_NoOut_Merge" , "_Cutoff_NoOut", 1)
  635. ))
  636. names(merged_data_files) <- c("filename", "field_name", "suffix", "merge_order")
  637. merged_data_files[["merge_order"]] <- as.numeric(merged_data_files[["merge_order"]])
  639. data_for_merging <- lapply(merged_data_files[["filename"]],
  640. function(next_filename) {read_data_file(paste0(dirs[["input_other"]], output_file_prefix, next_filename), na="NA")}
  641. )
  642. names(data_for_merging) <- paste0(output_file_prefix, merged_data_files[["field_name"]])
  645. ## comprehensively merge brain data files with matching files
  646. data_for_analysis <- list()
  648. for(df_1_name in paste0(output_file_prefix, c("LH_Thickness","RH_Thickness","LH_Area","RH_Area","Vol"))) {
  649. df_1 <- brain_data[["grey"]][[df_1_name]]
  651. for(df_2_idx in 1:nrow(merged_data_files)) {
  652. df_2 <- data_for_merging[[paste0(output_file_prefix, merged_data_files[df_2_idx,"field_name"])]]
  654. next_merge_order <- merged_data_files[df_2_idx,"merge_order"]
  656. if(next_merge_order==1) {
  657. next_merged_df <- merge(df_2, df_1, by="id")
  658. } else if(next_merge_order==2) {
  659. next_merged_df <- merge(df_1, df_2, by="id")
  660. } else {
  661. next_merged_df <- NULL
  662. }
  664. data_for_analysis[[paste0(df_1_name, merged_data_files[df_2_idx,"suffix"])]] <- next_merged_df
  665. }
  666. }
  669. #WM
  671. #Creating SLF Avg variables
  673. # store info about col names to match ("match" col) and the name of the column
  674. # in which averages should be stored ("stored")
  675. slf_avg_var_ids <- list()
  678. slf_avg_var_ids[["wm"]] <- data.frame(rbind(
  679. c("l_FA" , "l_FA" )
  680. , c("l_MD" , "l_MD" )
  681. , c("r_FA" , "r_FA" )
  682. , c("r_MD" , "r_MD" )
  683. ))
  684. slf_avg_var_ids[["noddi"]] <- data.frame(rbind(
  685. c("l_Neurite_density", "l_Neurite_density" )
  686. , c("l_ODI" , "l_ODI")
  687. , c("r_Neurite_density", "r_Neurite_density" )
  688. , c("r_ODI" , "r_ODI")
  689. ))
  690. for(next_df_name in names(slf_avg_var_ids)) {
  691. colnames(slf_avg_var_ids[[next_df_name]]) <- c("match","store")
  692. }
  694. ##MD and FA
  695. brain_data[["white"]][[paste0(output_file_prefix, "MDFA")]] <- create_slf_avg_vars(brain_data[["white"]][[paste0(output_file_prefix, "MDFA")]] , slf_avg_var_ids[["wm"]] )
  696. matched_group_data.merged[[paste0(output_file_prefix, "Matched_WM")]] <- create_slf_avg_vars(matched_group_data.merged[[paste0(output_file_prefix, "Matched_WM")]] , slf_avg_var_ids[["wm"]] )
  698. ##ND and ODI
  699. brain_data[["white"]][[paste0(output_file_prefix, "Noddi")]] <- create_slf_avg_vars(brain_data[["white"]][[paste0(output_file_prefix, "Noddi")]] , slf_avg_var_ids[["noddi"]])
  700. matched_group_data.merged[[paste0(output_file_prefix, "Matched_Noddi")]] <- create_slf_avg_vars(matched_group_data.merged[[paste0(output_file_prefix, "Matched_Noddi")]], slf_avg_var_ids[["noddi"]])
  703. ## FOR DEBUGGING PURPOSES:
  704. ## just printing out the col names being averaged here:
  705. slf_summary_str <- c()
  706. slf_summary_str <- c(slf_summary_str, paste0(paste(c("brain_data", "white", paste0(output_file_prefix, "MDFA")), collapse="$"), "$",
  707. create_slf_avg_vars_fake(brain_data[["white"]][[paste0(output_file_prefix, "MDFA")]], slf_avg_var_ids[["wm"]])))
  708. slf_summary_str <- c(slf_summary_str, paste0(paste(c("matched_group_data.merged", paste0(output_file_prefix, "Matched_WM")), collapse="$"), "$",
  709. create_slf_avg_vars_fake(matched_group_data.merged[[paste0(output_file_prefix, "Matched_WM")]], slf_avg_var_ids[["wm"]])))
  710. slf_summary_str <- c(slf_summary_str, paste0(paste(c("brain_data", "white", paste0(output_file_prefix, "Noddi")), collapse="$"), "$",
  711. create_slf_avg_vars_fake(brain_data[["white"]][[paste0(output_file_prefix, "Noddi")]], slf_avg_var_ids[["noddi"]])))
  712. slf_summary_str <- c(slf_summary_str, paste0(paste(c("matched_group_data.merged", paste0(output_file_prefix, "Matched_Noddi")), collapse="$"), "$",
  713. create_slf_avg_vars_fake(matched_group_data.merged[[paste0(output_file_prefix, "Matched_Noddi")]], slf_avg_var_ids[["noddi"]])))
  715. slf_summary_str <- do.call(rbind.data.frame, strsplit(slf_summary_str, "="))
  716. colnames(slf_summary_str) <- c("store","avg")
  718. # slf_summary_str[c(2,4,1,3,9,11,10,12),] # order of rowMeans averaging
  719. # slf_summary_str[c(1,3,2,4,9,11,10,12),] # order of mod1-8
  722. data_for_merging[[paste0(output_file_prefix, "MDFA_Mixed_Merge")]] <- read_data_file(paste0(dirs[["input_other"]], output_file_prefix, "MDFA_Mixed_ToMerge.csv"), na = "NA")
  723. data_for_analysis[[paste0(output_file_prefix, "MDFA" )]] <- merge(brain_data[["white"]][[paste0(output_file_prefix, "MDFA" )]], data_for_merging[[paste0(output_file_prefix, "MDFA_Mixed_Merge")]])
  724. data_for_analysis[[paste0(output_file_prefix, "Noddi")]] <- merge(brain_data[["white"]][[paste0(output_file_prefix, "Noddi")]], data_for_merging[[paste0(output_file_prefix, "MDFA_Mixed_Merge")]])
  725. ```
  728. ## Preparing Brain Data
  729. ```{r}
  731. dvs <- list()
  733. # Most of these (all except Vol) can be determined algorithmically from column names in data files
  735. dvs[["lh_area" ]] <- data_for_analysis[[paste0(output_file_prefix, "LH_Area" )]] %>% colnames %>% .[grepl("^lh_" , .)]
  737. dvs[["rh_area" ]] <- data_for_analysis[[paste0(output_file_prefix, "RH_Area" )]] %>% colnames %>% .[grepl("^rh_" , .)]
  739. dvs[["lh_thick"]] <- data_for_analysis[[paste0(output_file_prefix, "LH_Thickness")]] %>% colnames %>% .[grepl("^lh_.*_thickness$", .)]
  741. dvs[["rh_thick"]] <- data_for_analysis[[paste0(output_file_prefix, "RH_Thickness")]] %>% colnames %>% .[grepl("^rh_.*_thickness$", .)]
  743. # can't determine this one algorithmically
  744. dvs[["vol" ]] <- c("Brain_Stem", "CC_Anterior", "CC_Central", "CC_Mid_Anterior", "CC_Mid_Posterior", "CC_Posterior", "Left_Accumbens_area" , "Left_Amygdala", "Left_Caudate", "Left_Cerebellum_Cortex" , "Left_Cerebellum_White_Matter", "Left_choroid_plexus", "Left_Hippocampus", "Left_Pallidum" , "Left_Putamen" , "Left_Thalamus", "Optic_Chiasm", "Right_Accumbens_area" , "Right_Amygdala" , "Right_Caudate" , "Right_Cerebellum_Cortex" , "Right_Cerebellum_White_Matter", "Right_choroid_plexus", "Right_Hippocampus" , "Right_Lateral_Ventricle" , "Right_Pallidum", "Right_Putamen", "Right_Thalamus", "TotalGrayVol")
  746. dvs[["mdfa" ]] <- data_for_analysis[[paste0(output_file_prefix, "MDFA" )]] %>% colnames %>% .[(grepl("_FA$" , .) | grepl("_MD$" , .)) & !grepl("^slf[0-9]", .)]
  748. dvs[["noddi" ]] <- data_for_analysis[[paste0(output_file_prefix, "Noddi" )]] %>% colnames %>% .[(grepl("_Neurite_density$", .) | grepl("_ODI$", .)) & !grepl("^slf[0-9]", .)]
  750. ```
  752. ```{r All Models}
  754. # define predictors which are used as pred_others
  755. cnr_pred <- "cnr_c"
  756. control_preds <- c("age_c","sex","iq_c")
  757. all_preds <- paste(c(control_preds, cnr_pred), collapse=" + ")
  759. # create a data frame with all parameters needed, then call fit.lm() via an apply function to call everything at once!
  760. models_to_fit <- data.frame(rbind(
  761. # df dvs pred_pval pred_others id_text
  763. ## grey matter
  764. c("LH_Area" , "lh_area" , "Residuals_BehavioralRegression", cnr_pred , "LH_Area_Regression" )
  765. , c("Matched_GM" , "lh_area" , "PC_TD" , all_preds , "LH_Area_PC_TD" )
  766. , c("LH_Area_Cutoff_NoOut" , "lh_area" , "PC_TD_NoOut" , all_preds , "LH_Area_PC_TD_NoOut" )
  767. , c("Matched_GM" , "lh_area" , "PC_PD" , all_preds , "LH_Area_PC_PD" )
  768. , c("LH_Area_Cutoff_NoOut" , "lh_area" , "PC_PD_NoOut" , all_preds , "LH_Area_PC_PD_NoOut" )
  769. , c("Matched_GM" , "lh_area" , "PC_All" , all_preds , "LH_Area_PC_All" )
  770. , c("LH_Area_Cutoff_NoOut" , "lh_area" , "PC_All_NoOut" , all_preds , "LH_Area_PC_All_NoOut" )
  771. , c("LH_Area" , "lh_area" , "UPC_UGC" , all_preds , "LH_Area_UPC_UGC" )
  772. , c("LH_Area" , "lh_area" , "UPC_EC" , all_preds , "LH_Area_UPC_EC" )
  773. , c("LH_Area" , "lh_area" , "UPC_All" , all_preds , "LH_Area_UPC_All" )
  775. , c("RH_Area" , "rh_area" , "Residuals_BehavioralRegression", cnr_pred , "RH_Area_Regression" )
  776. , c("RH_Area_Cutoff" , "rh_area" , "PC_TD" , all_preds , "RH_Area_PC_TD" )
  777. , c("RH_Area_Cutoff_NoOut" , "rh_area" , "PC_TD_NoOut" , all_preds , "RH_Area_PC_TD_NoOut" )
  778. , c("RH_Area_Cutoff" , "rh_area" , "PC_PD" , all_preds , "RH_Area_PC_PD" )
  779. , c("RH_Area_Cutoff_NoOut" , "rh_area" , "PC_PD_NoOut" , all_preds , "RH_Area_PC_PD_NoOut" )
  780. , c("RH_Area_Cutoff" , "rh_area" , "PC_All" , all_preds , "RH_Area_PC_All" )
  781. , c("RH_Area_Cutoff_NoOut" , "rh_area" , "PC_All_NoOut" , all_preds , "RH_Area_PC_All_NoOut" )
  782. , c("RH_Area" , "rh_area" , "UPC_UGC" , all_preds , "RH_Area_UPC_UGC" )
  783. , c("RH_Area" , "rh_area" , "UPC_EC" , all_preds , "RH_Area_UPC_EC" )
  784. , c("RH_Area" , "rh_area" , "UPC_All" , all_preds , "RH_Area_UPC_All" )
  786. , c("LH_Thickness" , "lh_thick", "Residuals_BehavioralRegression", cnr_pred , "LH_Thick_Regression" )
  787. , c("LH_Thickness_Cutoff" , "lh_thick", "PC_TD" , all_preds , "LH_Thick_PC_TD" )
  788. , c("LH_Thickness_Cutoff_NoOut", "lh_thick", "PC_TD_NoOut" , all_preds , "LH_Thick_PC_TD_NoOut" )
  789. , c("LH_Thickness_Cutoff" , "lh_thick", "PC_PD" , all_preds , "LH_Thick_PC_PD" )
  790. , c("LH_Thickness_Cutoff_NoOut", "lh_thick", "PC_PD_NoOut" , all_preds , "LH_Thick_PC_PD_NoOut" )
  791. , c("LH_Thickness_Cutoff" , "lh_thick", "PC_All" , all_preds , "LH_Thick_PC_All" )
  792. , c("LH_Thickness_Cutoff_NoOut", "lh_thick", "PC_All_NoOut" , all_preds , "LH_Thick_PC_All_NoOut")
  793. , c("LH_Thickness" , "lh_thick", "UPC_UGC" , all_preds , "LH_Thick_UPC_UGC" )
  794. , c("LH_Thickness" , "lh_thick", "UPC_EC" , all_preds , "LH_Thick_UPC_EC" )
  795. , c("LH_Thickness" , "lh_thick", "UPC_All" , all_preds , "LH_Thick_UPC_All" )
  797. , c("RH_Thickness" , "rh_thick", "Residuals_BehavioralRegression", cnr_pred , "RH_Thick_Regression" )
  798. , c("RH_Thickness_Cutoff" , "rh_thick", "PC_TD" , all_preds , "RH_Thick_PC_TD" )
  799. , c("RH_Thickness_Cutoff_NoOut", "rh_thick", "PC_TD_NoOut" , all_preds , "RH_Thick_PC_TD_NoOut" )
  800. , c("RH_Thickness_Cutoff" , "rh_thick", "PC_PD" , all_preds , "RH_Thick_PC_PD" )
  801. , c("RH_Thickness_Cutoff_NoOut", "rh_thick", "PC_PD_NoOut" , all_preds , "RH_Thick_PC_PD_NoOut" )
  802. , c("RH_Thickness_Cutoff" , "rh_thick", "PC_All" , all_preds , "RH_Thick_PC_All" )
  803. , c("RH_Thickness_Cutoff_NoOut", "rh_thick", "PC_All_NoOut" , all_preds , "RH_Thick_PC_All_NoOut")
  804. , c("RH_Thickness" , "rh_thick", "UPC_UGC" , all_preds , "RH_Thick_UPC_UGC" )
  805. , c("RH_Thickness" , "rh_thick", "UPC_EC" , all_preds , "RH_Thick_UPC_EC" )
  806. , c("RH_Thickness" , "rh_thick", "UPC_All" , all_preds , "RH_Thick_UPC_All" )
  808. , c("Vol" , "vol" , "Residuals_BehavioralRegression", cnr_pred , "Vol_Regression" )
  809. , c("Vol_Cutoff" , "vol" , "PC_TD" , all_preds , "Vol_PC_TD" )
  810. , c("Vol_Cutoff_NoOut" , "vol" , "PC_TD_NoOut" , all_preds , "Vol_PC_TD_NoOut" )
  811. , c("Vol_Cutoff" , "vol" , "PC_PD" , all_preds , "Vol_PC_PD" )
  812. , c("Vol_Cutoff_NoOut" , "vol" , "PC_PD_NoOut" , all_preds , "Vol_PC_PD_NoOut" )
  813. , c("Vol_Cutoff" , "vol" , "PC_All" , all_preds , "Vol_PC_All" )
  814. , c("Vol_Cutoff_NoOut" , "vol" , "PC_All_NoOut" , all_preds , "Vol_PC_All_NoOut" )
  815. , c("Vol" , "vol" , "UPC_UGC" , all_preds , "Vol_UPC_UGC" )
  816. , c("Vol" , "vol" , "UPC_EC" , all_preds , "Vol_UPC_EC" )
  817. , c("Vol" , "vol" , "UPC_All" , all_preds , "Vol_UPC_All" )
  819. ## white matter
  820. , c("Matched_WM" , "mdfa" , "PC_PD" , all_preds , "MDFA_PC_PD" )
  821. , c("Matched_WM" , "mdfa" , "PC_TD" , all_preds , "MDFA_PC_TD" )
  822. , c("Matched_WM" , "mdfa" , "PC_All" , all_preds , "MDFA_PC_All" ) # NOT THE SAME DATA AS KELLY'S
  823. , c("MDFA" , "mdfa" , "UPC_EC" , all_preds , "MDFA_UPC_EC" )
  824. , c("MDFA" , "mdfa" , "UPC_UGC" , all_preds , "MDFA_UPC_UGC" )
  825. , c("MDFA" , "mdfa" , "UPC_All" , all_preds , "MDFA_UPC_All" )
  826. , c("MDFA" , "mdfa" , "Residuals_BehavioralRegression", cnr_pred , "MDFA_Regression" )
  828. ## NODDI
  829. , c("Matched_Noddi" , "noddi" , "PC_PD" , all_preds , "Noddi_PC_PD" )
  830. , c("Matched_Noddi" , "noddi" , "PC_TD" , all_preds , "Noddi_PC_TD" )
  831. , c("Matched_Noddi" , "noddi" , "PC_All" , all_preds , "Noddi_PC_All" )
  832. , c("Noddi" , "noddi" , "UPC_EC" , all_preds , "Noddi_UPC_EC" )
  833. , c("Noddi" , "noddi" , "UPC_UGC" , all_preds , "Noddi_UPC_UGC" )
  834. , c("Noddi" , "noddi" , "UPC_All" , all_preds , "Noddi_UPC_All" )
  835. , c("Noddi" , "noddi" , "Residuals_BehavioralRegression", cnr_pred , "Noddi_Regression" )
  836. ))
  837. names(models_to_fit) <- c("df", "dvs", "pred_pval", "preds_other", "id_text")
  839. ## COMMENTED OUT FOR DEBUGGING
  840. # # run all models in models_to_fit
  841. # all_models <- lapply(1:nrow(models_to_fit), function(next_model_idx) {
  842. # # extract info about next model
  843. # next_model_info <- models_to_fit[next_model_idx,]
  844. # # if next df name contains the word "Matched", look for it in matched_group_data.merged; otherwise, look for it in data_for_analysis
  845. # next_model_df <- ifelse(grepl("Matched", next_model_info[["df"]]), "matched_group_data.merged", "data_for_analysis")
  846. # # fit next set of models
  847. # next_fit_models <- fit_lms(
  848. # df = get(next_model_df)[[paste0(output_file_prefix, next_model_info[["df"]])]],
  849. # dvs = dvs[[models_to_fit[next_model_idx,"dvs"]]],
  850. # pred_pval = next_model_info[["pred_pval"]],
  851. # preds_other = strsplit(next_model_info[["preds_other"]], " \\+ ")[[1]],
  852. # id_text = next_model_info[["id_text"]]
  853. # )
  854. # return(next_fit_models)
  855. # })
  856. # names(all_models) <- models_to_fit[["id_text"]]
  859. # ENABLED FOR DEBUGGING
  860. all_models <- list()
  861. for(next_model_idx in 1:nrow(models_to_fit)) {
  862. print(paste0(next_model_idx, " / ", nrow(models_to_fit), "..."))
  863. # extract info about next model
  864. next_model_info <- models_to_fit[next_model_idx,]
  865. # if next df name contains the word "Matched", look for it in matched_group_data.merged; otherwise, look for it in data_for_analysis
  866. next_model_df <- ifelse(grepl("Matched", next_model_info[["df"]]), "matched_group_data.merged", "data_for_analysis")
  867. # fit next set of models
  868. next_fit_models <- fit_lms(
  869. df = get(next_model_df)[[paste0(output_file_prefix, next_model_info[["df"]])]],
  870. dvs = dvs[[models_to_fit[next_model_idx,"dvs"]]],
  871. pred_pval = next_model_info[["pred_pval"]],
  872. preds_other = strsplit(next_model_info[["preds_other"]], " \\+ ")[[1]],
  873. id_text = next_model_info[["id_text"]]
  874. )
  875. all_models[[next_model_info[["id_text"]]]] <- next_fit_models
  876. }
  880. # extract key features of linear results into a single df
  881. all_linear_results <- do.call(rbind.data.frame, lapply(names(all_models), function(model) {
  882. data.frame(model, all_models[[model]][["rsq.linear"]], p=all_models[[model]][["p"]])
  883. }))
  886. # extract all quadratic results into a single df
  887. models_with_quad_results_binary <- sapply(names(all_models), function(x) {"quad_results" %in% names(all_models[[x]])})
  888. all_quad_results <- do.call(rbind.data.frame, lapply(names(all_models)[models_with_quad_results_binary], function(model) {
  889. data.frame(model, all_models[[model]][["quad_results"]])
  890. }))
  892. all_rsq_results <- do.call(rbind.data.frame, lapply(names(all_models), function(model) {
  893. data.frame(model, all_models[[model]][["rsq.linear"]])
  894. }))
  895. all_rsq_results <- merge(all_rsq_results, all_quad_results[,colnames(all_quad_results) %in% colnames(all_rsq_results) | grepl("^rsq\\.",colnames(all_quad_results))], all.x=TRUE)
  898. #all_quad_results # contains all p-values from comparing linear & quadratic models
  900. # shows proportion of p-values <= .05 separately for each model
  901. sapply(unique(all_quad_results[["model"]]), function(model) {mean(all_quad_results[all_quad_results[["model"]]==model,"p"] <= .05)})
  903. # figure of same
  904. all_quad_results[["model"]] <- factor(all_quad_results[["model"]], unique(all_quad_results[["model"]]))
  905. all_quad_results[["p_thresh"]] <- as.factor(ifelse(all_quad_results[["p"]] <= .05, 1, 0))
  907. # for MDFA: subdivide based on l vs. r (or neither), and MD vs. FA
  908. # for NODDI: subdivide based on l vs. r (or neither), and ODI vs. NeuriteDensity
  909. all_quad_results.for_plotting <- all_quad_results
  910. all_quad_results.for_plotting[["model_group"]] <- gsub("^LH_", "l_", gsub("^RH_", "r_", gsub("_Regression$","",all_quad_results.for_plotting[["model"]])))
  911. # reorder measure and hemisphere designations...
  912. all_quad_results.for_plotting[["model_group"]] <- sapply(all_quad_results.for_plotting[["model_group"]], function(x) {paste(rev(strsplit(x, "_")[[1]]), collapse="_")})
  914. # oy!
  915. mdfa_split <- sapply(all_quad_results.for_plotting[all_quad_results.for_plotting[["model_group"]]=="MDFA","dv"], function(x) {strsplit(x, "_")})
  916. mdfa_num_fields <- sapply(mdfa_split, function(x) {length(x)})
  917. mdfa_groups <- sapply(1:length(mdfa_split), function(y) {paste(mdfa_split[[y]][mdfa_num_fields[y]:2], collapse="_")})
  918. all_quad_results.for_plotting[all_quad_results.for_plotting[["model_group"]]=="MDFA","model_group"] <- mdfa_groups
  920. noddi_split <- sapply(gsub("Neurite_density","NeuriteDensity",all_quad_results.for_plotting[all_quad_results.for_plotting[["model_group"]]=="Noddi","dv"]), function(x) {strsplit(x, "_")})
  921. noddi_num_fields <- sapply(noddi_split, function(x) {length(x)})
  922. noddi_groups <- sapply(1:length(noddi_split), function(y) {paste(noddi_split[[y]][noddi_num_fields[y]:2], collapse="_")})
  923. all_quad_results.for_plotting[all_quad_results.for_plotting[["model_group"]]=="Noddi","model_group"] <- noddi_groups
  925. ordered_model_groups <- c("Area_l", "Area_r", "Thick_l", "Thick_r", "Vol", "FA_l", "FA", "FA_r", "MD_l", "MD", "MD_r", "NeuriteDensity_l", "NeuriteDensity", "NeuriteDensity_r", "ODI_l", "ODI", "ODI_r")
  926. all_quad_results.for_plotting[["model_group"]] <- factor(all_quad_results.for_plotting[["model_group"]], ordered_model_groups)
  927. # determine where to add vertical lines for visual breaks in the more complex figure
  928. model_group_prefixes <- sapply(strsplit(ordered_model_groups, "_"), function(x) {x[[1]]})
  929. model_group_breaks <- which(model_group_prefixes[2:length(model_group_prefixes)]!=model_group_prefixes[1:(length(model_group_prefixes)-1)])+0.5
  931. # plot p-values by model (simpler)
  932. quad_p.plot <- ggplot(all_quad_results, aes(x=model, y=p, color=p_thresh)) +
  933. ggforce::geom_sina(size=1, scale=F, method="density", maxwidth=0.3, position=position_identity()) +
  934. scale_color_manual(values=c("black","red")) +
  935. geom_hline(yintercept=0.05, linetype="dashed", color="black") +
  936. theme_bw() +
  937. xlab("Model") +
  938. ylab("p")
  939. quartz(width=11, height=5); quad_p.plot
  941. # plot p-values by model, measure, hemisphere... (more complex)
  942. quad_p.for_plotting.plot <- ggplot(all_quad_results.for_plotting, aes(x=model_group, y=p, color=p_thresh)) +
  943. ggforce::geom_sina(size=1, scale=F, method="density", maxwidth=0.2, position=position_identity()) +
  944. scale_color_manual(values=c("black","red")) +
  945. geom_hline(yintercept=0.05, linetype="dashed", color="black") +
  946. geom_vline(xintercept=model_group_breaks, linetype="dashed", color="black") +
  947. theme_bw() +
  948. xlab("Model") +
  949. ylab("p")
  950. quartz(width=18, height=5); quad_p.for_plotting.plot
  954. ```
  957. ## How to access stored models
  958. ```{r Access model objects}
  960. # all_models is a list with 64 objects:
  961. # LH_Area_Regression, LH_Area_PCTD, ... RH_Area_Regression, RH_Area_PCTD, ...
  962. # LH_Thick_Regression, LH_Thick_PCTD, ... RH_Thick_Regression, RH_Thick_PCTD, ...
  963. # Vol_Regression, Vol_PCTD, ... MDFA_PCPD, MDFA_PCTD, ... Noddi_PCPD, Noddi_PCTD, ...
  964. # One of those can be accessed via (e.g.) all_models[["LH_Area_Regression"]]
  965. # That is itself a list with 3 objects:
  966. # "model" (a list of 75 fitted models); can be accessed via all_models[["LH_Area_Regression"]][["model"]]
  967. # "p" (a vector of 75 p-values, 1 per model); can be accessed via all_models[["LH_Area_Regression"]][["p"]]
  968. # "quad_results" (a data frame of results from quadratic tests); can be accessed via all_models[["LH_Area_Regression"]][["quad_results"]]
  969. # Unfortunately, the names of the model fields had to be changed to work with mtable(). So, at the moment, they're all "lm1", "lm2", ... "lm75".
  970. # You can see this via names(all_models[["LH_Area_Regression"]][["model"]])
  971. # So, in order to match those up with names which are actually useful, you can grab the list of dependent variables which they correspond to.
  972. # Those can be found in dvs[["lh_area"]] which returns a list of 75 area names. This means that, for example,
  973. # all_models[["LH_Area_Regression"]][["model"]][[1]] (or, equivalently, all_models[["LH_Area_Regression"]][["model"]][["lm1"]] )
  974. # will return the LH area regression model for dvs[["lh_area"]][[1]] , which is "lh_G_and_S_frontomargin_area" .
  975. # You can see the list of different sets of DVs via names(dvs) ,
  976. # and the list of everything that was matched up for model fitting is in models_to_fit
  978. # We can also make an enormous df of model p-values, as follows...
  979. all_model_results.flat <- suppressWarnings(lapply(1:nrow(models_to_fit), function(idx) {data.frame(subset(models_to_fit[idx,], select=-dvs), dv=dvs[[models_to_fit[idx,"dvs"]]], p=all_models[[idx]][["p"]])})) # warnings regard row names
  980. all_model_results.flat <- do.call(rbind.data.frame, all_model_results.flat)
  982. # all_model_results.flat now has everything in one place!
  984. all_rsq_results <- merge(all_model_results.flat, all_rsq_results, by.x=c("id_text","dv"), by.y=c("model","dv"))
  986. ```
  990. #FDR Correction for Multiple Comparisons
  991. ```{r MultipleComparisonsCorrection}
  993. dv_pcorr_groupings <- unique.data.frame(rbind(all_linear_results[,c("model","dv")],all_quad_results[,c("model","dv")])) # quad should be a subset, but just in case...
  995. model_suffixes <- c(
  996. "Regression"
  997. , "UPC_UGC"
  998. , "UPC_EC"
  999. , "UPC_All"
  1000. , "PC_PD"
  1001. , "PC_TD"
  1002. , "PC_All"
  1003. , "PC_PD_NoOut"
  1004. , "PC_TD_NoOut"
  1005. , "PC_All_NoOut"
  1006. )
  1008. dv_pcorr_groupings[["model_type"]] <- sapply(dv_pcorr_groupings[["model"]], function(x) {
  1009. match_idx <- which(sapply(model_suffixes, function(y) { grepl(paste0("_", y, "$"), x) }))
  1010. ifelse(length(match_idx)==1, return(model_suffixes[match_idx]), stop(paste0("Not exactly 1 match! : ", x, " has ", length(match_idx), " matches")))
  1011. })
  1013. dv_pcorr_groupings[["dv_type"]] <- sapply(strsplit(dv_pcorr_groupings[["dv"]], "_"), function(x) {x[length(x)]})
  1014. dv_pcorr_groupings[["dv_type"]][!dv_pcorr_groupings[["dv_type"]] %in% c("area","thickness","MD","FA","ODI","density")] <- "volume"
  1015. dv_pcorr_groupings[["pcorr_group_id"]] <- with(dv_pcorr_groupings, paste0(model_type, "_", dv_type))
  1017. # count how many models are in each group
  1018. unique_model_groups <- unique.data.frame(dv_pcorr_groupings[,c("model_type","dv_type","pcorr_group_id")])
  1019. unique_model_groups[["num_models_in_group"]] <- sapply(unique_model_groups[["pcorr_group_id"]], function(x) {sum(dv_pcorr_groupings[["pcorr_group_id"]]==x)})
  1021. # show # models per group and # of groups
  1022. dplyr::rename(aggregate(pcorr_group_id ~ dv_type + num_models_in_group, data=unique_model_groups, FUN=length), num_groups="pcorr_group_id")
  1025. ```
  1027. ```{r Explore sig. quadratic models}
  1029. # quad correction
  1030. all_quad_results <- merge(all_quad_results, dv_pcorr_groupings[,c("model","dv","pcorr_group_id")])
  1032. quad_summary <- do.call(rbind.data.frame, lapply(unique(all_quad_results[["pcorr_group_id"]]), function(pcorr_group_id) {
  1033. matching_rows <- all_quad_results[all_quad_results[["pcorr_group_id"]]==pcorr_group_id,]
  1034. c(nrow(matching_rows), sum(matching_rows[["p"]] < .05))}))
  1035. colnames(quad_summary) <- c("n","sig")
  1037. # correct p-values in groups
  1038. for(next_group in unique(all_quad_results[["pcorr_group_id"]])) {
  1039. temp_p_list <- all_quad_results[all_quad_results[["pcorr_group_id"]]==next_group,"p"]
  1040. all_quad_results[all_quad_results[["pcorr_group_id"]]==next_group,"p.corr"] <- p.adjust(temp_p_list, method="fdr")
  1041. all_quad_results[all_quad_results[["pcorr_group_id"]]==next_group,"n.corr"] <- length(temp_p_list)
  1042. }
  1044. sig_quad_results <- all_quad_results[all_quad_results[["p.corr"]] < .05,] # sig. post-correction
  1046. # plot post-correction sig. p-values by model
  1047. all_sig_quad_plots <- list()
  1048. for(next_sig_quad_result_idx in 1:nrow(sig_quad_results)) {
  1049. next_sig_quad_model_info <- sig_quad_results[next_sig_quad_result_idx,]
  1052. next_sig_quad_model <- with(next_sig_quad_model_info, all_models[[as.character(model)]][["model"]][[dv]]) # (1) this is the linear model, not the quad model; (2) do I need this?
  1053. #next_sig_quad_model <- with(next_sig_quad_model_info, all_models[[as.character(model)]][["quad_models"]][[as.character(dv)]]) # oy
  1055. #next_sig_quad_plot <- plot_grouped_vs_continuous(fitted_lm=next_sig_quad_model, pred_var=c("Residuals_BehavioralRegression", "Residuals_BehavioralRegression.quad"))
  1056. next_sig_quad_plot <- plot_grouped_vs_continuous(fitted_lm=next_sig_quad_model, pred_var="Residuals_BehavioralRegression", degree=c(1,2))
  1058. next_sig_quad_plot <- next_sig_quad_plot + ggtitle(next_sig_quad_model_info[["dv"]])
  1060. all_sig_quad_plots[[next_sig_quad_result_idx]] <- next_sig_quad_plot
  1061. }
  1063. all_sig_quad_plots.together <- grid.arrange(grobs=all_sig_quad_plots, ncol=5)
  1064. ggsave(filename=paste0(super_dir, "all_sig_quad_plots.png"), plot=all_sig_quad_plots.together, width=30, height=15, units="in")
  1068. all_sig_quad_plots.quad_models <- list()
  1069. for(next_sig_quad_result_idx in 1:nrow(sig_quad_results)) {
  1070. next_sig_quad_model_info <- sig_quad_results[next_sig_quad_result_idx,]
  1071. # next_sig_quad_model <- with(next_sig_quad_model_info, all_models[[as.character(model)]][["model"]][[dv]]) # (1) this is the linear model, not the quad model; (2) do I need this?
  1072. next_sig_quad_model <- with(next_sig_quad_model_info, all_models[[as.character(model)]][["quad_models"]][[as.character(dv)]]) # oy
  1073. all_sig_quad_plots.quad_models[[next_sig_quad_result_idx]] <- next_sig_quad_model
  1074. }
  1076. all_sig_quad_plots.quad_models <- lapply(1:length(all_sig_quad_plots.quad_models), function(x) {coef(summary(all_sig_quad_plots.quad_models[[x]]))})
  1077. sig_quad_results[["quad_coeff"]] <- sapply(1:length(all_sig_quad_plots.quad_models), function(x) {all_sig_quad_plots.quad_models[[x]][grepl("\\.quad$", rownames(all_sig_quad_plots.quad_models[[x]])),"Estimate"]})
  1079. sig_quad_results
  1081. ```
  1083. ```{r Explore sig. linear models}
  1085. # linear correction
  1086. all_linear_results <- merge(all_linear_results, dv_pcorr_groupings[,c("model","dv","pcorr_group_id")])
  1088. linear_summary <- do.call(rbind.data.frame, lapply(unique(all_linear_results[["pcorr_group_id"]]), function(pcorr_group_id) {
  1089. matching_rows <- all_linear_results[all_linear_results[["pcorr_group_id"]]==pcorr_group_id,]
  1090. c(nrow(matching_rows), sum(matching_rows[["p"]] < .05))}))
  1091. colnames(linear_summary) <- c("n","sig")
  1093. # correct p-values in groups
  1094. for(next_group in unique(all_linear_results[["pcorr_group_id"]])) {
  1095. temp_p_list <- all_linear_results[all_linear_results[["pcorr_group_id"]]==next_group,"p"]
  1096. all_linear_results[all_linear_results[["pcorr_group_id"]]==next_group,"p.corr"] <- p.adjust(temp_p_list, method="fdr")
  1097. all_linear_results[all_linear_results[["pcorr_group_id"]]==next_group,"n.corr"] <- length(temp_p_list)
  1098. }
  1101. sig_linear_results <- all_linear_results[all_linear_results[["p.corr"]] < .05,] # sig. post-correction
  1103. # get p-values for quadratic plots for these same models, to identify models that only show linear but not quadratic components...
  1104. compare_linear_and_quad <- merge(rename(sig_linear_results, p.linear="p", p.corr.linear="p.corr", n.corr.linear="n.corr"), subset(rename(all_quad_results, p.quad="p", p.corr.quad="p.corr", n.corr.quad="n.corr"), select=-c(F,df,p_thresh)), by.all=c("model","dv","group"), all.x=TRUE) # lol
  1105. # remove models for which a quadratic model was not fit
  1106. compare_linear_and_quad <- compare_linear_and_quad[!is.na(compare_linear_and_quad[["rsq.quad"]]),] # only 6!
  1107. # remove models with sig. quad fits
  1108. compare_linear_and_quad <- compare_linear_and_quad[compare_linear_and_quad[["p.corr.quad"]] > .05,]
  1111. # plot quadratic fits for these models (not a typo!)
  1112. all_sig_linear_nonsig_quad_plots <- list()
  1113. for(next_result_idx in 1:nrow(compare_linear_and_quad)) {
  1114. next_model_info <- compare_linear_and_quad[next_result_idx,]
  1116. next_model <- with(next_model_info, all_models[[as.character(model)]][["model"]][[dv]])
  1118. next_plot <- plot_grouped_vs_continuous(fitted_lm=next_model, pred_var="Residuals_BehavioralRegression", degree=c(1,2))
  1119. next_plot <- next_plot + ggtitle(next_model_info[["dv"]])
  1121. all_sig_linear_nonsig_quad_plots[[next_result_idx]] <- next_plot
  1122. }
  1124. all_sig_linear_nonsig_quad_plots <- grid.arrange(grobs=all_sig_linear_nonsig_quad_plots, ncol=3)
  1125. ggsave(filename=paste0(super_dir, "all_sig_linear_nonsig_quad_plots.png"), plot=all_sig_linear_nonsig_quad_plots, width=18, height=10, units="in")
  1127. ```
  1130. ```{r Figs for paper}
  1132. # added 12/23/25
  1134. # First: Sinaplots for the AF, SLF, and R ILF for the UPC vs EC contrast (that one has the largest effect sizes) with mean lines by group (all MD)
  1136. regions_for_sinaplots <- c(af_l_MD="MD - Left AF", af_r_MD="MD - Right AF", ilf_r_MD="MD - Right ILF", slfavg_l_MD="MD - Left SLF", slfavg_r_MD="MD - Right SLF")
  1138. models_for_sinaplots <- data.frame(model="MDFA_UPC_EC", dv=names(regions_for_sinaplots), desc=unname(regions_for_sinaplots), pred.crit="UPC_EC")
  1140. # overly complicated function for computing multiple stats on a given df
  1141. agg_df <- function(char_funs=c("min","max","mean","sd"), df, x.col, y.col) {
  1142. # uses aggregate() to apply each of the functions listed (in character format) in char_funs
  1143. # to formula(y.col ~ [all of the vars in x.cols]) in data frame df, then returns the merged result
  1144. resids_by_group <- Reduce(function(x,y) merge(x, y, by="group", all=TRUE), lapply(char_funs, function(fun) {agg <- aggregate(formula(paste0(y.col, " ~ ", paste(x.col, collapse=" + "))), FUN=eval(parse(text=fun)), data=df); colnames(agg) <- c(x.col,fun); return(agg) })) # !!!
  1145. return(resids_by_group)
  1146. }
  1148. group_colors <- viridis(4)[2:3] # group colors are defined here -- first for UPC, second for EC
  1149. names(group_colors) <- c("UPC","EC")
  1151. x_margin <- 0.15
  1153. all_sinaplots <- list()
  1155. set.seed(123) # sinaplots involve random jitter; setting the seed for reproducibility
  1157. for(next_model_for_sinaplot_idx in 1:nrow(models_for_sinaplots)) {
  1158. next_model <- with(models_for_sinaplots[next_model_for_sinaplot_idx,], all_models[[model]][["model"]][[dv]]) # extract appropriate model
  1160. # fit subset models to get residualized predictor & DV
  1161. next_model.resids <- get_model_resids(next_model, models_for_sinaplots[next_model_for_sinaplot_idx,"pred.crit"])
  1163. # combine df with residualized DV
  1164. next_model.for_plotting <- data.frame(next_model[["model"]], dv.resid=next_model.resids[["subset_lm.resids"]], pred.resids=next_model.resids[["pred.resids"]])
  1166. # this hard-codes in something which is earlier assumed to be a variable... the key col may not be UPC_EC
  1167. next_model.for_plotting[["group"]] <- as.factor(with(next_model.for_plotting, ifelse(UPC_EC==0, "UPC", "EC"))) # it's either 0 or 1
  1168. next_model.for_plotting[["group"]] <- factor(next_model.for_plotting[["group"]], c("UPC", "EC"))
  1170. # compute means for separate plotting
  1171. next_model.for_plotting.means <- dplyr::rename(agg_df(char_funs=c("mean","sd","length"), df=next_model.for_plotting, x.col=c("group"), y="dv.resid"), n="length")
  1172. next_model.for_plotting.means[["se"]] <- with(next_model.for_plotting.means, sd/sqrt(n))
  1173. next_model.for_plotting.means <- next_model.for_plotting.means[order(next_model.for_plotting.means[["group"]]),] # order (important!)
  1175. next_plot <- ggplot(next_model.for_plotting, aes(x=group, y=dv.resid, color=group)) +
  1176. ggforce::geom_sina(size=1, scale=F, method="density", maxwidth=0.3, position=position_identity()) +
  1177. scale_color_manual(values=group_colors) +
  1178. theme_bw() +
  1179. xlab("Group") +
  1180. ylab(paste0(models_for_sinaplots[next_model_for_sinaplot_idx,"desc"], " (resid.)")) +
  1181. theme(panel.grid.major.x=element_blank(), legend.position="none") # no vertical lines
  1183. # add means
  1184. for(next_group_mean_idx in 1:nrow(next_model.for_plotting.means)) {
  1185. next_plot <- next_plot +
  1186. geom_segment(color="black", linewidth=0.5, # horizontal line at y-mean
  1187. x = with(next_model.for_plotting.means[next_group_mean_idx,], next_group_mean_idx-x_margin),
  1188. xend = with(next_model.for_plotting.means[next_group_mean_idx,], next_group_mean_idx+x_margin),
  1189. y = with(next_model.for_plotting.means[next_group_mean_idx,], mean ),
  1190. yend = with(next_model.for_plotting.means[next_group_mean_idx,], mean )) +
  1191. geom_segment(color="black", linewidth=1.5, # vertical line to show +/- 1 SE variability
  1192. x = next_group_mean_idx, xend=next_group_mean_idx,
  1193. y = with(next_model.for_plotting.means[next_group_mean_idx,], mean-se),
  1194. yend = with(next_model.for_plotting.means[next_group_mean_idx,], mean+se))
  1195. }
  1197. # store plot in list
  1198. all_sinaplots[[next_model_for_sinaplot_idx]] <- next_plot
  1199. }
  1201. # combine plots and save
  1202. all_sinaplots.together <- grid.arrange(grobs=all_sinaplots, ncol=3)
  1203. ggsave(filename=paste0(super_dir, "all_sinaplots.png"), plot=all_sinaplots.together, width=9, height=6, units="in")
  1207. # Second: Scatterplots with a nonlinear fit line for thel AF (af_l_fa) and the right SLF (slf_r_Neurite_density) with colors from the groups
  1209. # note: the second one is not sig. (p = .051)
  1210. models_for_quadplots <- data.frame(rbind(c("MDFA_Regression" , "af_l_FA" , "FA - Left AF" , "Residuals_BehavioralRegression"),
  1211. c("Noddi_Regression", "slfavg_r_Neurite_density", "NODDI - Right SLF", "Residuals_BehavioralRegression")))
  1212. colnames(models_for_quadplots) <- c("model","dv","desc","pred.crit")
  1215. all_quad_results[(all_quad_results[["dv"]] %in% models_for_quadplots[["dv"]]),] # show info for these two models
  1218. # generate plots for these quadratic models
  1219. all_paper_quad_plots <- list()
  1220. for(next_paper_quad_result_idx in 1:nrow(models_for_quadplots)) {
  1221. next_paper_quad_model_info <- models_for_quadplots[next_paper_quad_result_idx,]
  1223. next_paper_quad_model <- with(next_paper_quad_model_info, all_models[[as.character(model)]][["model"]][[dv]])
  1224. next_paper_quad_plot <- plot_grouped_vs_continuous(fitted_lm=next_paper_quad_model, pred_var=next_paper_quad_model_info[["pred.crit"]], degree=c(2), linetype=c(1)) # to show both linear AND quad, change degree to c(1,2)
  1226. next_paper_quad_plot <- next_paper_quad_plot +
  1227. xlab("Comprehension regression residuals (resid.)") +
  1228. ylab(paste0(next_paper_quad_model_info[["desc"]], " (resid.)"))
  1231. all_paper_quad_plots[[next_paper_quad_result_idx]] <- next_paper_quad_plot
  1232. }
  1234. all_paper_quad_plots.together <- grid.arrange(grobs=all_paper_quad_plots, ncol=2)
  1235. ggsave(filename=paste0(super_dir, "all_paper_quad_plots.png"), plot=all_paper_quad_plots.together, width=12, height=5, units="in")
  1237. ```

Mahaffy_Masters_Real_DK__KM_20251223 (1).Rmd, no license · at the source

Overview

Authors: Kelly Mahaffy1,2, Nabin Koirala2,3,4, Daniel Kleinman2, Nicole Landi1,2,3
ORCID iDs: Kelly Mahaffy, Nabin Koirala, Daniel Kleinman, Nicole Landi
  1. Department of Psychological Sciences, University of Connecticut, Storrs, CT, USA
  2. Child Study Center, Yale School of Medicine, New Haven, CT, USA
  3. Brain Imaging Research Core, University of Connecticut, Storrs, CT, USA
  4. The Nathan S. Klein Institute for Psychiatric Research, Orangeburg, NY, USA
Institutions: University of Connecticut (United States); Yale University (United States); Nathan Kline Institute for Psychiatric Research (United States)
Journal: Neurobiology of language (Cambridge, Mass.), volume 7, article NOL.a.265
Dates: received 7 May 2025; accepted 14 April 2026; published online 1 July 2026
Type: Research article · Language: English
License: CC BY
Identifiers: DOI 10.1162/nol.a.265 · PMID 42492006 · PMCID PMC13379302 · OpenAlex W4412638000
Open access: gold, a free copy (OpenAlex)
Preprint: doi:10.1101/2025.07.23.666369
Status: code verified
Categories: structural MRI / diffusion (modality), human (organism), cognitive (subfield)
Methods: Statistics, Connectivity, fMRI & imaging, Machine learning, Preprocessing
Keywords: educational neuroscience, NODDI (neurite orientation dispersion and density imaging), poor comprehension, reading comprehension
Topic: Reading and Literacy Development (Developmental and Educational Psychology, Psychology), according to OpenAlex
Funding: NIDCD NIH HHS (T32 DC017703)
Citations: not cited yet (Europe PMC); 102 references in the paper

Abstract

Poor comprehenders (PCs) have typical word reading skill and intelligence but poorer than expected reading comprehension. While the prevalence of PCs is similar to that of poor decoders (individuals who struggle fluently converting written text into spoken language), less is known about the neurobiological substrates of poor comprehension. Extant studies have found small differences in gray matter volume between PCs and poor or typically reading peers. However, a detailed quantification of cortical morphometric features and white matter integrity remains unexplored. Data from 2,100 children (1,200 with imaging data), aged 8–16 years were analyzed to determine if there is a distinct neuroanatomy associated with poor reading comprehension across three common methods for classifying PCs. We computed gray matter volume, cortical thickness, and surface area, as well as white matter measures (mean diffusivity, fractional anisotropy, neurite orientation, and neurite density) for PCs and compared these measures with those of poor decoders and typical readers. Results revealed small but widespread white matter differences, but no gray matter differences between PCs and other readers. PCs showed decreased white matter integrity (increased mean diffusivity, decreased neurite density) in tracts previously associated with reading, including the superior longitudinal fasciculus and inferior longitudinal fasciculus, and in tracts that have been associated with cognitive performance such as the uncinate fasciculus. These results suggest that diffuse structural connectivity differences may underlie reading comprehension weaknesses in the face of intact decoding skills. This is consistent with the behavioral profile of PCs who exhibit a broad pattern of subclinical impairments in language and integrative cognitive processes.

Reproduced under the paper's license (CC BY), from the paper cited above.

Repository

Its files are read in the Code ↔ Paper reader above, with 2 matches between paragraphs and lines of code.

OSF 2ng57

License: none: the authors keep all their rights
State: the link answers, verified on 27 September 2026
Evidence: files inventoried
Languages: R (1)
Size: 14 files, 1 script
Software Heritage: not checked
Found in: “DATA AVAILABILITY STATEMENT”
Holds: 1 notebook
Not found: README, license file, CITATION.cff, environment file, tests, continuous integration, documentation
Tools: ggplot2 (1 file), ggpubr (1 file), psych (1 file), reshape2 (1 file), tidyverse (1 file)
Availability: 1 check, the latest on 27 September 2026: the link answers (HTTP 200)
  • 27 September 2026: the link answers (HTTP 200)
1 file
At the source: osf.io/2ng57/overview

The paper's code and data availability statement is in the Data section.

Tracing map

Proposed by the machine: these links were found in the paper and verified at the source, without human review. The map will receive a Zenodo DOI once one of the paper's authors has validated it with their ORCID.

What the map holds:

  • 1 repository of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
  • 1 script, each with its path and the digest of its content;
  • 2 matches between paragraphs of the paper and lines of the code (method lexical-v1);
  • neither the text of the paper nor the code itself.

Its JSON (tracing-map.json) is deposited on Zenodo with its DOI once the map is validated.

Data

No dataset and no data link were found in the paper.

Data availability statement

All data and code used for this project, along with project preregistration, are publicly available in an OSF repository: https://osf.io/2ng57/overview.

Reproduced under the paper's license (CC BY), from the paper cited above.

Versions

The history of this record: each version stored by the harvester or made by a correction of its authors or of the maintainers of its code, and what changed in its facts. The texts of the paper (its abstract, its availability statements) are not part of it; versions that changed only those are not listed.

Version 1, 27 September 2026: the first record

Recorded: type, language, journal, volume, pages, dates, 4 authors, 4 keywords, 1 funder, 96 references.

Cite

This paper

Mahaffy, K., Koirala, N., Kleinman, D., & Landi, N. (2026). Structural Brain Correlates of Poor Reading Comprehension. Neurobiology of language (Cambridge, Mass.), 7, NOL.a.265. https://doi.org/10.1162/nol.a.265

BibTeX

@article{mahaffy2026structural,
author = {Mahaffy, Kelly and Koirala, Nabin and Kleinman, Daniel and Landi, Nicole},
title = {{Structural Brain Correlates of Poor Reading Comprehension}},
journal = {Neurobiology of language (Cambridge, Mass.)},
year = {2026},
month = jul,
volume = {7},
pages = {NOL.a.265},
publisher = {MIT Press},
issn = {2641-4368},
doi = {10.1162/nol.a.265},
url = {https://doi.org/10.1162/nol.a.265},
pmid = {42492006},
pmcid = {PMC13379302}
}

RIS

TY - JOUR
AU - Mahaffy, Kelly
AU - Koirala, Nabin
AU - Kleinman, Daniel
AU - Landi, Nicole
TI - Structural Brain Correlates of Poor Reading Comprehension
T2 - Neurobiology of language (Cambridge, Mass.)
J2 - Neurobiol Lang (Camb)
PY - 2026
DA - 2026/07/01
VL - 7
SP - NOL.a.265
SN - 2641-4368
PB - MIT Press
DO - 10.1162/nol.a.265
UR - https://doi.org/10.1162/nol.a.265
LA - en
ER -

CSL-JSON

{
"id": "10.1162/nol.a.265",
"type": "article-journal",
"title": "Structural Brain Correlates of Poor Reading Comprehension",
"container-title": "Neurobiology of language (Cambridge, Mass.)",
"author": [
{
"family": "Mahaffy",
"given": "Kelly"
},
{
"family": "Koirala",
"given": "Nabin"
},
{
"family": "Kleinman",
"given": "Daniel"
},
{
"family": "Landi",
"given": "Nicole"
}
],
"container-title-short": "Neurobiol Lang (Camb)",
"volume": "7",
"page": "NOL.a.265",
"DOI": "10.1162/nol.a.265",
"PMID": "42492006",
"PMCID": "PMC13379302",
"ISSN": "2641-4368",
"publisher": "MIT Press",
"URL": "https://doi.org/10.1162/nol.a.265",
"language": "en",
"issued": {
"date-parts": [
[
2026,
7,
1
]
]
}
}

The tracing map gets a citation of its own once an author has validated it and it has a DOI.

Similar papers

The papers with a page that share the most with this one: the tools found in their code, their categories, datasets, cited references and authors, the rarest counting most.

[1] doi:10.1162/imag.a.1303
Composite reaction time and variability correlate with whole-brain white-matter characteristics.
Journal: Imaging neuroscience (Cambridge, Mass.)
In common: structural MRI / diffusion, 14 references
[2] doi:10.3389/frai.2026.1771088 [code]
Few-shot deployment of pretrained MRI transformers in brain imaging tasks.
Journal: Frontiers in artificial intelligence
In common: structural MRI / diffusion, 13 references
[3] doi:10.3389/fpain.2026.1850836 [code]
Diffusion tensor imaging in chronic tension-type headache.
Journal: Frontiers in pain research (Lausanne, Switzerland)
In common: structural MRI / diffusion, 11 references
[4] doi:10.1038/s41593-026-02359-0 [code]
The cross-site reproducibility of MRI morphometric phenotypes in psychiatric disorders.
Journal: Nature neuroscience
In common: structural MRI / diffusion, 9 references
[5] doi:10.1038/s41467-026-73262-2 [code]
Robust but independent sex differences in human brain function, structure, and behavior.
Journal: Nature communications
In common: reshape2, ggpubr, ggplot2, 1 other tool, structural MRI / diffusion, 4 references
[6] doi:10.1016/j.nicl.2026.104012 [code]
Structural-functional multilayer brain network properties and outcome of combined repetitive transcranial magnetic stimulation and psychotherapy for obsessive-compulsive disorder.
Journal: NeuroImage. Clinical
In common: psych, reshape2, ggpubr, 2 other tools, structural MRI / diffusion, 2 references
[7] doi:10.7554/elife.108109 [code]
Multimodal MRI marker of cognition explains the association between cognition and mental health in the UK Biobank.
Journal: eLife
In common: psych, ggplot2, structural MRI / diffusion, cognitive, 5 references
[8] doi:10.1162/netn.a.536 [code]
Increased sensitivity in identifying language-related functional connectivity using jackknife resampling analyses.
Journal: Network neuroscience (Cambridge, Mass.)
In common: psych, reshape2, ggpubr, 2 other tools, 2 references
[9] doi:10.1038/s41593-026-02363-4 [code]
Cortical thickness changes precede high levels of amyloid by at least 7 years.
Journal: Nature neuroscience
In common: ggpubr, ggplot2, tidyverse, structural MRI / diffusion, 4 references
[10] doi:10.1007/s00415-026-14034-2
Similarities and differences in the late-onset GM2 gangliosidoses: Tay-Sachs and Sandhoff diseases.
Journal: Journal of neurology
In common: structural MRI / diffusion, 7 references

Contribute

The authors of this paper can claim it, correct its record and validate its tracing map, and the maintainers of its code (its owner, or a public member of its organization) correct what it says of their repository; anyone signed in can ask for its removal. Every request goes to OSCR's own machine, which answers it; your account page follows them.

Sign in with ORCID to claim this paper as one of its authors, correct its record or validate its tracing map: when the paper's metadata lists your ORCID iD, you are recognized at once. Maintainers of its code: sign in with GitHub, then claim the repository on your account page.

Request its removal

To ask OSCR to remove this record, the copies of its authors' scripts or its tracing map, use the removal request page: signed in, you say who you are, what to remove and why, then review and confirm the request. Published rules decide every request (how).

Discussion, reproductions, activity

Discussion: questions and error reports about this paper and its code, from signed-in readers and its authors. It opens with sign-in.

Reproductions: reports from readers who ran the authors' code: what they reproduced, with which environment, commit and data. It opens with sign-in.

Activity: what happens around this paper: new versions of its record, its map's validation, discussions and reproductions. It opens with sign-in.