# Calculating estimated AR parametersp <-ncol(similar_mvts[[1]])# Function for summarize estimated AR(1) parameters for participantsar_summary <-function(x) {summary(diag(matrix(ar(similar_mvts[[x]], order.max =1, aic =FALSE)$ar, nrow = p))) %>% broom::tidy() %>% flextable::flextable()}map(.x =1:nCandidates, .f = ar_summary)
Individual-Level Analysis
We demonstrate the utility of our proposed inferential methods by analyzing resting-state functional magnetic resonance imaging data from the Autism Brain Imaging Data Exchange (ABIDE) initiative (Craddock et al., 2013). We analyze imaging data from a single site provided by Stanford University that was preprocessed using the Data Processing Assistant for Resting-State fMRI pipeline. Within this site, we consider data from female participants diagnosed with autism spectrum disorder who were approximately 10 years old at the time of data acquisition. The Automated Anatomical Labeling atlas was used to extract data for 116 regions of interest, but we focus our analysis on ten brain regions in the Default Mode Network. See Craddock et al. (2013) for more details on data acquisition and preprocessing steps.
Since the resting-state fMRI data consists of 175 volumes for 10 regions of interest, and the average autoregressive parameter estimated for the data for each participant is about 0.40, the simulation setting that most aligns with these settings has \(N=100\) observations, \(p=10\) variables, and \(\phi = 0.4\) for the autoregressive parameter.
Partial correlation plot
Code
# Vector of brain regionsregions <- stringr::str_remove_all(roiNames$old, pattern ="[.]")fullRegions <-str_replace_all(str_replace_all(roiNames$new, pattern =" R", replacement =" Right"), pattern =" L", replacement =" Left")# Function for creating plots of group-level networksanalyzeNetwork <-function(myGraph, rois = regions) { p <-ncol(myGraph[[1]])# Creating "network" objects myNet <- network::network(myGraph >0, directed =FALSE,vertex.attr =list(rois),vertex.attrnames =list("ROI"))# Setting edge weightsset.edge.value(myNet, "weight", myGraph)# Creating network data frames for plottingroiCoords <-ggnetwork(myNet, layout ="fruchtermanreingold") %>% dplyr::distinct(x, y, ROI) %>% dplyr::select(x, y) %>%as.matrix()ggnetDF <-ggnetwork(myNet, layout = roiCoords)# Plotting networksggNet <- ggnetDF %>%ggplot(aes(x = x, y = y, xend = xend, yend = yend))# Plotting darkest edges firstfor(w insort(ggnetDF$weight, decreasing =FALSE)) { ggNet <- ggNet +geom_edges(data = ggnetDF %>% dplyr::filter(weight == w), aes(color = weight), linetype ="solid", linewidth =1.5)}# Adding nodes ggNet <- ggNet +geom_nodes(aes(x, y), color ="steelblue", size =11) +geom_nodetext(aes(label = ROI),size =2.5, color ="white") +scale_color_gradient2(high ="black", mid ="white", low ="white",limits=c(0,1), breaks =seq(0, 1, by =0.20), midpoint =0.50,guide =guide_colorbar(frame.colour ="black", ticks.colour ="black")) +labs(color ="Proportion of individuals with an edge present:") +theme_facet() +theme(legend.position ="bottom")return(ggNet)}# Function for calculating matrix of p-values based on covariance matrices & pc estimatescalcPvalueMat <-function(covMat, pcorrs) { se2s <-diag(covMat) pc2s <- pcorrs^2 zs <- pc2s / se2s chiCritical <-qchisq(p =1-0.05, df =1, ncp =0) q <-length(pcorrs) p <- (1+sqrt(1+8*q)) /2 pvals <-pchisq(q= zs, ncp =0, df =1, lower.tail =FALSE) pvalMat <-diag(rep(0, p)) pvalMat[upper.tri(pvalMat, diag =FALSE)] <- pvalsreturn(pvalMat +t(pvalMat))}library(ggcorrplot)library(grid)p <-ncol(similar_mvts[[1]])q <-choose(p, 2)triInds <-upper.tri(diag(p))# Creating list of vectors of partial correlationspartialCorrs <-map(1:nCandidates, .f =function(part_num) { pcCov::corrMat_cpp(tsData = similar_mvts[[part_num]], partial =TRUE)[triInds]})# Obtaining list of matrices of p-values based on Wald testpvalueMats <-map(1:nCandidates, .f =function(part_num) { pMat <-calcPvalueMat(covMat = asymp_covs[[part_num]], pcorrs = partialCorrs[[part_num]])colnames(pMat) <-rownames(pMat) <- regionsreturn(pMat) })# Creating common plot functionmatrixPlot <-function(myMatrix, myTitle) { ggcorrplot::ggcorrplot(myMatrix, method ="circle", type ="upper") +scale_fill_gradient2(high ="#0072B2", low ="white",limits=c(0, 10), breaks =seq(0, 10, by =2)) +labs(x ="", y ="", fill =expression("-log"[10]*"(p-value)"), title = myTitle) +theme_bw() +theme(legend.key.height =unit(1, "cm"), text =element_text(face ="bold"),axis.text.y =element_text(size =6, face ="bold"),axis.text.x =element_text(size =6, face ="bold"),plot.title =element_text(size =15, face ="bold"))}# Plots for partial correlationsmyMats <-map(pvalueMats, .f =function(pmat) { pmat <--log(pmat, base =10)# Thresholding for plottingfor(r in1:nrow(pmat)) {for(column in1:ncol(pmat)) { pmat[r, column] <-min(pmat[r, column], 10) }}return(pmat) })myCorrGGs <-map(myMats, .f =function(myMat) { myGG <-matrixPlot(myMat, myTitle ="Individual Wald Tests of Partial Correlations") myGG[["plot_env"]][["p"]][["layers"]][[1]][["mapping"]][["size"]] <-1return(myGG)})# Saving plotsmap(1:nCandidates, .f =function(p_num) {if(!file.exists(paste0("correlation_plot_", p_num, ".png"))) {ggsave(filename =paste0("correlation_plot_", p_num, ".png"), plot = myCorrGGs[[p_num]], units ="mm",width =1344/8, height =960/8) }})
Results of Wald tests of individual partial correlations being 0 or not for each participant.
# Determining common statistically significant correlationscommonSigs <-Reduce(map(.x =1:nCandidates, .f =function(x) {p.adjust(pvalueMats[[x]][triInds], method ="holm", n = nCandidates*length(pvalueMats[[x]][triInds])) <0.05}), f ="*")for(s in1:sum(commonSigs)) {print(roiNames[regionInds[as.logical(commonSigs), ][s, ], ])}
After applying a multiplicity adjustment as proposed by Holm (1979), the partial correlation between the right and left cingulum posterior regions, and the partial correlation between the right and left precuneus regions were statistically significant for all participants at the \(5\%\) significance level.
References
Cameron Craddock, Yassine Benhajali, Carlton Chu, Francois Chouinard, Alan Evans, András Jakab, Budhachandra Singh Khundrakpam, John David Lewis, Qingyang Li, Michael Milham, Chaogan Yan, Pierre Bellec (2013). The Neuro Bureau Preprocessing Initiative: open sharing of preprocessed neuroimaging data and derivatives. In Neuroinformatics 2013, Stockholm, Sweden.
Holm, S. (1979). A simple sequentially rejective multiple test procedure. Scandinavian Journal of Statistics, 6, 65–70. https://www.jstor.org/stable/4615733.
