doc/FCnet_overview_of_package.R

In EBlini/FCnet: Analysis of Functional Connectivity matrices through elastic NETs

## ----setup, include=FALSE------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
knitr::opts_chunk$set(echo = TRUE)
options(width = 999)
Sys.setenv(LANG = "en")


## ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
library("FCnet") #for the main analysis routines
library("ggplot2") #beautiful depictions
library("gridExtra") #arrange plots
library("future.apply") #for speed and parallel computing


## ----eval= F-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
#  my_matrices= FCnet::loadFC()
#  

## ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
data("MeanFC") #the object MeanFC is now available

set.seed(1)

N_subs= 50 #number of participants/matrices
subjs_variability= 0.2 #variability (sd) between matrices

#this function creates several matrices by adding gaussian noise
#with sd "variability" to a reference matrix
m_start= simulateMat(mat = MeanFC, 
                     Nmat = N_subs, 
                     mat_variability = subjs_variability)


## ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
y= rnorm(N_subs, mean = 0, sd= 1) #the behavioral score

network1= 50:100 #network to perturb
network2= 200:250 #network to perturb

bias_multiplier= 0.6 #mean signal to inject (i.e. mean correlation with y)
bias_variability= 0.2 #variability - sd - around this mean signal

#this function biases a given network as a function 
#of a behavioral score y
m_bias= biasMat(matrices = m_start, 
                y = y, 
                network1 = network1, 
                network2 = network2,
                bias_multiplier = bias_multiplier,
                bias_variability = bias_variability,
                mat_variability = subjs_variability)
  

## ----fig.height=9--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
#this adds a grid on the plot corresponding to the
#network which was perturbed
grid_net= list(geom_rect(aes(xmin= min(network1), 
                         xmax= max(network1),
                         ymin= min(network2),
                         ymax= max(network2)),
                         size= 1.05, color= "black", 
                         fill= "transparent"))

ps= plotFC(m_start[[which.max(y)]], #we plot one specific matrix
           lim = c(-1, 1)) + ggtitle("Start") + grid_net

pb= plotFC(m_bias[[which.max(y)]], 
           lim = c(-1, 1)) + ggtitle("Bias") + grid_net

#arrange plots
grid.arrange(ps, pb)


## ----warning= F, message=F-----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
f_cor= function(lst, y){
  
  n= length(lst) 	   
  rc= dim(lst[[1]]) 	   
  ar1= array(unlist(lst), c(rc, n)) 	
  
  apply(ar1, c(1, 2), function(x)(cor(x, y))) 	         

}

cor_bias= f_cor(m_bias, y) #warnings due to vectors being identical (e.g. diagonal)

plotFC(cor_bias) + ggtitle("Correlation map") + grid_net


## ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
f_sd= function(lst){
  
  n= length(lst) 	   
  rc= dim(lst[[1]]) 	   
  ar1= array(unlist(lst), c(rc, n)) 	
  
  apply(ar1, c(1, 2), function(x)(sd(x))) 	         

}

bias_sd= f_sd(m_bias)

plotFC(bias_sd, limit = c(0.15, 0.25)) + ggtitle("Variability") + grid_net


## ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
rf_bias= reduce_featuresFC(FCmatrices = m_bias, 
                           Ncomp = 10) 


## ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
r_bias= FCnetLOO(y = y, 
                 x =  rf_bias, 
                 alpha= 0, #ridge regression
                 parallelLOO = F) 
r_bias$R2


## ----fig.width= 9, fig.height=9------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
plotFCnet(model = r_bias, 
          plot_labels = F) #avoid padding, otherwise useful

plotFCnet(model = r_bias,
          output = "coefficients") 


## ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
bp_bias= backprojectFCnet(coeffs = r_bias, 
                          reduce_features_object = rf_bias,
                          threshold = length(network1)*length(network2))

plotFC(bp_bias) + grid_net


## ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
bp_first= backprojectFCnet(coeffs = 1,  
                          reduce_features_object = rf_bias,
                          threshold = length(network1)*length(network2))

plotFC(bp_first) + grid_net


## ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
plan(multisession) #prepare for parallel computing
perm_bias= permutateLOO(y = y,
                        x= rf_bias, 
                        alpha = 0, 
                        parallelLOO = T, 
                        nperm = 100, 
                        model_R2 = r_bias)

knitr::kable(perm_bias$Summary)


## ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
sessionInfo()