Descriptions and example scripts for running network preprocessing and analysis functions contained in this repository. See paper for details when it comes out.
#From original location to server
rsync -avz -O --omit-dir-times --no-perms --include="7*" --include="7*/Learn?_PEpriorD.feat" --include="7*/mprage.nii" --include="7*/Rest?.nii" --include="7*/Learn?_PEpriorD.feat/filtered_func_data.nii.gz" --include="7*/Learn?_PEpriorD.feat/reg" --include="7*/Learn?_PEpriorD.feat/mc" --include="7*/Learn?_PEpriorD.feat/mc/*" --include="7*/Learn?_PEpriorD.feat/reg/*" --exclude="*" --exclude="*/*" --exclude="*/*/*" rgerraty@lovelace.psych.columbia.edu:/data/engine/juliet/adoles/ /danl/Harmon_dynCon/
#Just taking Learn & Rest preprocessed data to run community detection on elvis
rsync --dry-run -avz -O --omit-dir-times --no-perms --include="4*" --include="4*/Learn?_PEprior.feat" --include="4*/Rest" --include="4*/Learn?_PEprior.feat/36par+spikes.feat" --include="4*/Rest/Rest?.feat" --include="41*/Rest/Rest?.feat/36par+spikes.feat" --include="4*/Learn?_PEprior.feat/36par+spikes.feat/H-O_rois" --include="4*/Learn?_PEprior.feat/36par+spikes.feat/H-O_rois" --include="4*/Learn?_PEprior.feat/36par+spikes.feat/H-O_rois/conn_cells.mat" --include="4*/Rest/Rest?.feat/36par+spikes.feat/H-O_rois" --include="4*/Rest/Rest?.feat/36par+spikes.feat/H-O_rois/conn_cells.mat" --exclude="*" --exclude="*/*" --exclude="*/*/*" charmon@lovelace.psych.columbia.edu:/data/engine/rgerraty/learn_dyncon/ /danl/Harmon_dynCon/
%Remove --dry_run when ready to run
#From server to Habanero
rsync -avz -O --omit-dir-times --no-perms --include="7*" --include="7*/Learn?_PEpriorD.feat" --include="7*/Learn?_PEpriorD.feat/*" --include="7*/Learn?_PEpriorD.feat/reg" --include="7*/Learn?_PEpriorD.feat/mc" --include="7*/Learn?_PEpriorD.feat/mc/*" --include="7*/Learn?_PEpriorD.feat/reg/*" --exclude="*" --exclude="*/*" --exclude="*/*/*" charmon@lux.psych.columbia.edu:/danl/Harmon_dynCon/ /rigel/psych/users/cmh2228/dynCon/Prepare folder structure. Run fsl_anat to prepare anatomical scans. Run initial preprocessing of resting state data.
#lux
/danl/Harmon_dynCon/scripts/0.set_up_scripts
./0.1preprepFolderStructure.sh
#habanero
/rigel/psych/users/cmh2228/dynCon/scripts/
for n in 7* ; do sbatch fslAnat.sh $n ; done
for n in 7* ; do sbatch restFeat.sh $n ; done
#lux
/danl/Harmon_dynCon/scripts/1.preprocessing
./restReg.sh
Because of the known effect of motion on measures of connectivity, we followed up standard preprocessing in FSL with an extended nuisance regression. Affine transformation parameters from motion correction, CSF, white matter, and whole-brain signals are regressed against preprocessed 4D data, along with the squares, derivatives, and squared derivatives of these confounds. See Satterthwaite et al 2013 for details. Bash code:
For Task
#do this on lux
/danl/Harmon_dynCon/scripts/1.preprocessing
./1.0extract_confts
#do this on habanero
for i in /rigel/psych/users/cmh2228/dynCon/7*/Learn*/filtered_func_data.nii.gz;
do
subdir=$(dirname $i);
sbatch --export=arg1=$i,arg2=$subdir/36par+spikes.txt,arg3=/rigel/psych/users/cmh2228/dynCon/rl_flexibility/conf_reg_design.fsf,arg4=/rigel/psych/app/fsl/data,arg5=$subdir/conf_reg_design.fsf /rigel/psych/users/cmh2228/dynCon/scripts/1.preprocessing/1st_level_conf.sub.sh;
doneFor Rest
rsync -avz -O --omit-dir-times --no-perms --include="7*" --include="7*/Rest" --include="7*/Rest/Rest?.feat" --include="7*/Rest/Rest?.feat/*" --include="7*/Rest/Rest?.feat" --include="7*/Rest/Rest?.feat/logs" --include="7*/Rest/Rest?.feat" --include="7*/Rest/Rest?.feat/logs/*" --include="7*/Rest/Rest?.feat" --include="7*/Rest/Rest?.feat/mc" --include="7*/Rest/Rest?.feat" --include="7*/Rest/Rest?.feat/mc/*" --exclude="*" --exclude="*/*" --exclude="*/*/*" --exclude="*/*/*/*" cmh2228@habanero.rcs.columbia.edu:/rigel/psych/users/cmh2228/dynCon/ /danl/Harmon_dynCon/
#must register data before doing:
#do this on lux
/danl/Harmon_dynCon/scripts/1.preprocessing
./1.0extract_confts_rest
#do this on habanero
rsync -avz -O --omit-dir-times --no-perms --include="7*" --include="7*/Rest" --include="7*/Rest/Rest?.feat" --include="7*/Rest/Rest?.feat/*" --include="7*/Rest/Rest?.feat" --include="7*/Rest/Rest?.feat/logs" --include="7*/Rest/Rest?.feat" --include="7*/Rest/Rest?.feat/logs/*" --include="7*/Rest/Rest?.feat" --include="7*/Rest/Rest?.feat/mc" --include="7*/Rest/Rest?.feat" --include="7*/Rest/Rest?.feat/mc/*" --include="7*/Rest/Rest?.feat" --include="7*/Rest/Rest?.feat/reg" --include="7*/Rest/Rest?.feat" --include="7*/Rest/Rest?.feat/reg/*" --exclude="*" --exclude="*/*" --exclude="*/*/*" --exclude="*/*/*/*" charmon@lux.psych.columbia.edu:/danl/Harmon_dynCon/ /rigel/psych/users/cmh2228/dynCon/
rsync -avz -O --omit-dir-times --no-perms --include="7*" --include="7*/Rest" --include="7*/Rest/Rest?.feat" --include="7*/Rest/Rest?.feat/*" --include="7*/Rest/Rest?.feat" --include="7*/Rest/Rest?.feat/reg" --include="7*/Rest/Rest?.feat/reg/*" --exclude="*" --exclude="*/*" --exclude="*/*/*" --exclude="*/*/*/*" charmon@lux.psych.columbia.edu:/danl/Harmon_dynCon/ /rigel/psych/users/cmh2228/dynCon/
for i in /rigel/psych/users/cmh2228/dynCon/7*/Rest/Rest?.feat/filtered_func_data.nii.gz;
do
subdir=$(dirname $i);
sbatch --export=arg1=$i,arg2=$subdir/36par+spikes.txt,arg3=/rigel/psych/users/cmh2228/dynCon/rl_flexibility/conf_reg_design.fsf,arg4=/rigel/psych/app/fsl/data,arg5=$subdir/conf_reg_design.fsf /rigel/psych/users/cmh2228/dynCon/scripts/1.preprocessing/1st_level_conf.sub.sh;
done
After nuisance regression has been run, the residual timeseries needs to be transformed into standard space (in this case, MNI). Make sure fsl_anat has been run on each structural image first. The following bash code was used to perform these transformations:
for Task
#fnirt has already been run, just applying transformation
for i in /rigel/psych/users/cmh2228/dynCon/7*/Learn?_PEpriorD.feat;
do sbatch --export=arg1=$i,arg2=/rigel/psych/app/fsl/data /rigel/psych/users/cmh2228/dynCon/scripts/non_lin_reg_sub.sh;
done
for Rest
#fnirt has already been run, just applying transformation
for i in /rigel/psych/users/cmh2228/dynCon/7*/Rest/Rest?.feat;
do sbatch --export=arg1=$i,arg2=/rigel/psych/app/fsl/data /rigel/psych/users/cmh2228/dynCon/scripts/non_lin_reg_sub.sh;
done#Moving newly created extended preprocesses folders (36par+spikes.feat) to lux
#Practice - dry run displays the files that will be copied
rsync --dry-run -avz --include="7*" --include="7*/Learn?_PEpriorD.feat" --include="7*/Learn?_PEpriorD.feat/36*.feat" --include="7*/Learn?_PEpriorD.feat/36*.feat/stats" --include="7*/Learn?_PEpriorD.feat/36*.feat/logs/" --include="7*/Learn?_PEpriorD.feat/36*.feat/stats/*" --include="7*/Learn?_PEpriorD.feat/36*.feat/logs/*" --exclude="*" --exclude="*/*" --exclude="*/*/*" --exclude="*/*/*/*" --exclude="*/*/*/*/*" cmh2228@habanero.rcs.columbia.edu:/rigel/psych/users/cmh2228/dynCon/ /danl/Harmon_dynCon/
#Task
rsync -avz -O --omit-dir-times --no-perms --include="7*" --include="7*/Learn?_PEpriorD.feat" --include="7*/Learn?_PEpriorD.feat/36*.feat" --include="7*/Learn?_PEpriorD.feat/36*.feat/stats" --include="7*/Learn?_PEpriorD.feat/36*.feat/logs/" --include="7*/Learn?_PEpriorD.feat/36*.feat/stats/*" --include="7*/Learn?_PEpriorD.feat/36*.feat/logs/*" --exclude="*" --exclude="*/*" --exclude="*/*/*" --exclude="*/*/*/*" --exclude="*/*/*/*/*" cmh2228@habanero.rcs.columbia.edu:/rigel/psych/users/cmh2228/dynCon/ /danl/Harmon_dynCon/
#Rest
rsync -avz -O --omit-dir-times --no-perms --include="7*" --include="7*/Rest" --include="7*/Rest/Rest?.feat" --include="7*/Rest/Rest?.feat/36*.feat" --include="7*/Rest/Rest?.feat/36*.feat/*" --include="7*/Rest/Rest?.feat/36*.feat/stats" --include="7*/Rest/Rest?.feat/36*.feat/logs/" --include="7*/Rest/Rest?.feat/36*.feat/stats/*" --include="7*/Rest/Rest?.feat/36*.feat/logs/*" --exclude="*" --exclude="*/*" --exclude="*/*/*" --exclude="*/*/*/*" --exclude="*/*/*/*/*" --exclude="*/*/*/*/*/*" cmh2228@habanero.rcs.columbia.edu:/rigel/psych/users/cmh2228/dynCon/ /danl/Harmon_dynCon/
Once the preprocessed images have been registered, we extract mean timecourses for each Harvard-Oxford ROI, using the function extract_ROIs.sh. The output of this function is a timecourse for each ROI in the specified input folder, as well as a .txt file containing all of the ROIs. The bash code used to run this function on each learning block for each subject is below:
#scp Harvard-Oxford_ROIs folder from rgerraty
for i in /danl/Harmon_dynCon/7*/Learn?_PEprior.feat/36par+spikes.feat/;
do
#extract timeseries (mean or 1st eigenvector, see function) data from each ROI in ~/Harvard-Oxford_ROIs/
~/GitHub/rl_flexibility/extract_ROIs.sh $i/stats/res4d_std.nii.gz /danl/Harmon_dynCon/Harvard-Oxford_ROIs/ $i/H-O_rois/;
done
#Also can be run by running script ./extract_time_course.shfor Rest
nohup ./extract_time_course_rest.sh &
Connectivity between pairs of ROIs was measured by average magnitude squared coherence in the .06-.12 Hz band, computed in MATLAB. The code below calls a function for creating a coherence matrices in a specified frequency range for specified time windows (in this case 25 TRs, or 50 s). These are saved as a .mat file for multi-slice community detection.
for Task
screen -r
matlab -nosplash -nojvm
%Alternatively(and slower) matlab -nosplash -nodesktop
addpath ~/GitHub/rl_flexibility
%read in all subject/run ROI timeseries directories
[a,b]=system('ls -d /danl/Harmon_dynCon/7*/Learn?_PEpriorD.feat/36par+spikes.feat/H-O_rois');
%do separately above section followed by next section
c=strread(b,'%s');
for i=1:size(c,1)
%calculate coherence per time window from concatenated ROI file
filename=char(strcat(c(i),'/all_rois.txt'))
%need to specify filename, window length in TR, sampling rate, bandpass
conn_cell=coherence_by_block(filename,25,.5,.06,.12);
save(char(strcat(c(i),'/conn_cells')),'conn_cell')
end
Subject 713 has 2 missing ROIs (column 13 & 61) - this code removes those columns in the all_rois.txt file before running community detection Rerun calculating the coherence matrices having removed these two columns from all 4 Learning blocks for this subject After rerunning add columns and rows back in
addpath ~/GitHub/rl_flexibility
%read in all subject/run ROI timeseries directories
[a,b]=system('ls -d /danl/Harmon_dynCon/713/Learn4_PEpriorD.feat/36par+spikes.feat/H-O_rois');
%do separately above section followed by next section
c=strread(b,'%s');
incomp_rois=load(char(strcat(c(4),'/all_rois.txt')))
[row, col]=find(~incomp_rois)
column_missing=unique(col)
%This is one way to remove missing ROIs for the file that has 0's incomp_rois( :, ~any(incomp_rois,1) ) = []
%Below is to get rid of excluded coher matrices for all runs for this subject
%Creating new ROI txt files excluding NA columns
[a,b]=system('ls -d /danl/Harmon_dynCon/713/Learn?_PEpriorD.feat/36par+spikes.feat/H-O_rois');
c=strread(b,'%s');
for i=1:size(c,1)
incomp_rois=load(char(strcat(c(i),'/all_rois.txt')))
incomp_rois(:,[13 61]) = []
save(char(strcat(c(1),'/all_rois_incomp.txt')),'incomp_rois', '-ascii')
end
%calculating coherence matrices excluding NA rois
[a,b]=system('ls -d /danl/Harmon_dynCon/713/Learn?_PEpriorD.feat/36par+spikes.feat/H-O_rois');
%do separately above section followed by next section
c=strread(b,'%s');
for i=1:size(c,1)
%calculate coherence per time window from concatenated ROI file
filename=char(strcat(c(i),'/all_rois_incomp.txt'))
%need to specify filename, window length in TR, sampling rate, bandpass
conn_cell=coherence_by_block(filename,25,.5,.06,.12);
save(char(strcat(c(i),'/conn_cells')),'conn_cell')
end
for Rest
screen -r
matlab -nosplash -nojvm
%Alternatively(and slower) matlab -nosplash -nodesktop
addpath ~/GitHub/rl_flexibility
%read in all subject/run ROI timeseries directories
[a,b]=system('ls -d /danl/Harmon_dynCon/7*/Rest/Rest?.feat/36par+spikes.feat/H-O_rois');
%do separately above section followed by next section
c=strread(b,'%s');
for i=1:size(c,1)
%calculate coherence per time window from concatenated ROI file
filename=char(strcat(c(i),'/all_rois.txt'))
%need to specify filename, window length in TR, sampling rate, bandpass
conn_cell=coherence_by_block(filename,25,.5,.06,.12);
save(char(strcat(c(i),'/conn_cells')),'conn_cell')
end
Input coherence matrix for each block. Also need number of blocks, resolution and coupling parameters. In Matlab
Because Subject 713 has 2 missing ROIs - add them back to conn_cells here
forTask
#Do separately for 70* 71* and 72* by using:
screen
control^ a d
screen -r
matlab -nosplash -nodisplay
%need multi-slice, flexibility codes not yet on GitHub for network_diags to run
addpath ~/GitHub/rl_flexibility
addpath ~/GitHub/rl_flexibility/GenLouvain/
addpath ~/GitHub/rl_flexibility/Bassett_Code/
%read in data
[a,b]=system('ls -d /danl/Harmon_dynCon/7*/Learn?_PEpriorD.feat/36par+spikes.feat/H-O_rois/');
[a,b]=system('ls -d /danl/Harmon_dynCon/72*/Learn?_PEpriorD.feat/36par+spikes.feat/H-O_rois/');
[a,b]=system('ls -d /danl/Harmon_dynCon/713/Learn?_PEpriorD.feat/36par+spikes.feat/H-O_rois/');
%do the above first then do the below
c=sort(strread(b,'%s'));
%concatenate runs for each subject
numruns=4
k=1;
for j=1:size(c,1)/numruns
c(k)
conn_cell_cat=[];
for i=1:numruns
load(strcat(char(c(k-1+i)),'/conn_cells'))
conn_cell_cat=cat(3,conn_cell_cat,conn_cell)
end
%network_diags code:
%runs multi-slice community detection
%gives flexibility for each run
%also allegiance matrix (not using yet)
%need to specify number of blocks, simulations, coupling, resolution
[a_mat,flex]=network_diags(conn_cell_cat,4,500,1,1.1813) %This will be ammended *see below* for saving comm detect
save(char(strcat(c(k),'/../../../a_mat')),'a_mat')
save(char(strcat(c(k),'/../../../flex')),'flex')
%save(char(strcat(c(k),'/../../../prom')),'prom') %See below - the next 3 lines are added seperately in re-run comm detect
%save(char(strcat(c(k),'/../../../S_tmp')),'S_tmp')
%save(char(strcat(c(k),'/../../../Q_tmp')),'Q_tmp')
k=k+numruns;
endforRest
screen
control^ a d
screen -r
matlab -nosplash -nodisplay
%need multi-slice, flexibility codes not yet on GitHub for network_diags to run
addpath ~/GitHub/rl_flexibility
addpath ~/GitHub/rl_flexibility/GenLouvain/
addpath ~/GitHub/rl_flexibility/Bassett_Code/
%read in data
[a,b]=system('ls -d /danl/Harmon_dynCon/7*/Rest/Rest?.feat/36par+spikes.feat/H-O_rois/');
%do the above first then do the below
c=sort(strread(b,'%s'));
%concatenate runs for each subject
numruns=1
k=1;
for j=1:size(c,1)/numruns
c(k)
conn_cell=[];
for i=1:numruns
load(strcat(char(c(k-1+i)),'/conn_cells'))
conn_cell=cat(3,conn_cell,conn_cell)
end
[a_mat,flex]=network_diags(conn_cell,1,500,1,1.1813)
save(char(strcat(c(k),'/../../a_mat')),'a_mat')
save(char(strcat(c(k),'/../../flex')),'flex')
k=k+numruns;
end
To add back in the missing columns/rows for subject 713 before pulling all subj flexibility stats #Should be updated for a_mat.mat & flex.mat
%do this from a script in matlab in the Learn folder seperately to index properly - to add back NA columns before concatenating
load('/danl/Harmon_dynCon/713/flex.mat')
flexTEMP=flex
missing=[13 61; 61 13]
for i=1:length(missing)
%FOR inserting rows
len=size(flexTEMP,i);
A = [flexTEMP(1:missing(i)-1,:)];
B = NaN(1,4);
C = [flexTEMP(missing(i):end,:)];
flexTEMP=[A; B; C];
end
save('flex', 'flexTEMP')
load('a_mat.mat')
a_matTEMP=a_mat;
a_matTEST=zeros(110,110,32);
missing=[13 61; 61 13];
for k=1:size(a_mat,3)
a_matTEMP=a_mat(:,:,k)
missing=[13 61; 61 13]
for i=1:length(missing)
%For inserting rows
len=size(a_matTEMP,i);
A = [a_matTEMP(1:missing(i)-1,:)]';
B = NaN(1,len)';
C = [a_matTEMP(missing(i):end,:)]';
a_matTEMP1=[A B C];
%For inserting columns
len=size(a_matTEMP1,2)
A = [a_matTEMP1(1:missing(i)-1,:)];
B = NaN(1,len);
C = [a_matTEMP1(missing(i):end,:)];
a_matTEMP=[A; B; C];
end
a_matTEST(:,:,k)=a_matTEMP;
end
Re-Run multi-slice community detection and flexibility statistics to save Community Detection Numbers
Input coherence matrix for each block. Also need number of blocks, resolution and coupling parameters. In Matlab
%forTask
%need multi-slice, flexibility codes not yet on GitHub for network_diags to run
addpath ~/GitHub/rl_flexibility
addpath ~/GitHub/rl_flexibility/GenLouvain/
addpath ~/GitHub/rl_flexibility/Bassett_Code/
%read in data - to make faster break up into different screens
[a,b]=system('ls -d /danl/Harmon_dynCon/7*/Learn?_PEpriorD.feat/36par+spikes.feat/H-O_rois/');
[a,b]=system('ls -d /danl/Harmon_dynCon/71*/Learn?_PEpriorD.feat/36par+spikes.feat/H-O_rois/');
[a,b]=system('ls -d /danl/Harmon_dynCon/72*/Learn?_PEpriorD.feat/36par+spikes.feat/H-O_rois/');
[a,b]=system('ls -d /danl/Harmon_dynCon/4*/Learn?_PEprior.feat/36par+spikes.feat/H-O_rois/');
[a,b]=system('ls -d /danl/Harmon_dynCon/41*/Learn?_PEprior.feat/36par+spikes.feat/H-O_rois/');
[a,b]=system('ls -d /danl/Harmon_dynCon/42*/Learn?_PEprior.feat/36par+spikes.feat/H-O_rois/');
%do the above first then do the below
c=sort(strread(b,'%s'));
sim=500;
maxComm = zeros(1,sim)';
%concatenate runs for each subject
numruns=4
k=1;
for j=1:size(c,1)/numruns
c(k)
conn_cell_cat=[];
for i=1:numruns
load(strcat(char(c(k-1+i)),'/conn_cells'))
conn_cell_cat=cat(3,conn_cell_cat,conn_cell)
end
%network_diags code:
%runs multi-slice community detection
%gives flexibility for each run
%also allegiance matrix (not using yet)
%need to specify number of blocks, simulations, coupling, resolution
[a_mat,flex,prom,S_tmp,Q_tmp]=network_diags(conn_cell_cat,4,500,1,1.1813)
save(char(strcat(c(k),'/../../../a_mat_rerun1')),'a_mat') %Remove _rerun1 when doing for the first time
save(char(strcat(c(k),'/../../../flex_rerun1')),'flex') %Remove _rerun1 when doing for the first time
save(char(strcat(c(k),'/../../../prom')),'prom')
save(char(strcat(c(k),'/../../../S_tmp')),'S_tmp')
save(char(strcat(c(k),'/../../../Q_tmp')),'Q_tmp')
for i=1:size(S_tmp,3)
maxComm(i)=max(unique(S_tmp(:,:,i)));
end
save(char(strcat(c(k),'/../../../maxComm')),'maxComm')
k=k+numruns;
end
%To get spreadsheet of max Communities
%Adolescents
cd /danl/Harmon_dynCon
[a,b]=system('ls -d /danl/Harmon_dynCon/7*');
c=sort(strread(b,'%s'));
commNum=zeros(size(c,1), 32);
for i=1:size(c)
load(strcat(char(c(i)),'/S_tmp.mat'))
commNum(i,:)=mean(max(S_tmp),3)
csvwrite('AdolesCommNum.csv', commNum)
end
%Adults
cd /danl/Harmon_dynCon
[a,b]=system('ls -d /danl/Harmon_dynCon/4*/');
c=sort(strread(b,'%s'));
commNum=zeros(size(c,1), 32);
for i=1:size(c)
load(strcat(char(c(i)),'/S_tmp.mat'))
commNum(i,:)=mean(max(S_tmp),3)
csvwrite('AdultsCommNum.csv', commNum)
end
%forRest
%need multi-slice, flexibility codes not yet on GitHub for network_diags to run
addpath ~/GitHub/rl_flexibility
addpath ~/GitHub/rl_flexibility/GenLouvain/
addpath ~/GitHub/rl_flexibility/Bassett_Code/
%read in data
[a,b]=system('ls -d /danl/Harmon_dynCon/7*/Rest/Rest?.feat/36par+spikes.feat/H-O_rois/');
[a,b]=system('ls -d /danl/Harmon_dynCon/71*/Rest/Rest?.feat/36par+spikes.feat/H-O_rois/');
[a,b]=system('ls -d /danl/Harmon_dynCon/72*/Rest/Rest?.feat/36par+spikes.feat/H-O_rois/');
[a,b]=system('ls -d /danl/Harmon_dynCon/4*/Rest/Rest?.feat/36par+spikes.feat/H-O_rois/');
[a,b]=system('ls -d /danl/Harmon_dynCon/41*/Rest/Rest?.feat/36par+spikes.feat/H-O_rois/');
[a,b]=system('ls -d /danl/Harmon_dynCon/42*/Rest/Rest?.feat/36par+spikes.feat/H-O_rois/');
%do the above first then do the below
c=sort(strread(b,'%s'));
sim=500;
maxComm = zeros(1,sim)';
%concatenate runs for each subject
numruns=1
k=1;
for j=1:size(c,1)/numruns
c(k)
conn_cell_cat=[];
for i=1:numruns
load(strcat(char(c(k-1+i)),'/conn_cells'))
conn_cell_cat=cat(3,conn_cell_cat,conn_cell) %This does nothing for resting state b/c numruns=1
end
%network_diags code:
%runs multi-slice community detection
%gives flexibility for each run
%also allegiance matrix (not using yet)
%need to specify number of blocks, simulations, coupling, resolution
[a_mat,flex, prom, S_tmp, Q_tmp]=network_diags(conn_cell_cat,1,500,1,1.1813)
save(char(strcat(c(k),'/../../a_mat_rerun1')),'a_mat')
save(char(strcat(c(k),'/../../flex_rerun1')),'flex')
save(char(strcat(c(k),'/../../prom')),'prom')
save(char(strcat(c(k),'/../../S_tmp')),'S_tmp')
save(char(strcat(c(k),'/../../Q_tmp')),'Q_tmp')
for i=1:size(S_tmp,3)
maxComm(i)=max(unique(S_tmp(:,:,i)));
end
save(char(strcat(c(k),'/../../maxComm')),'maxComm')
k=k+numruns;
end
%To get spreadsheet of max Communities
%Adolescents
cd /danl/Harmon_dynCon
[a,b]=system('ls -d /danl/Harmon_dynCon/7*/Rest/Rest1.feat/');
c=sort(strread(b,'%s'));
commNum=zeros(size(c,1), 6);
for i=1:size(c)
load(strcat(char(c(i)),'/S_tmp.mat'))
commNum(i,:)=mean(max(S_tmp),3)
csvwrite('AdolesCommNumRest1.csv', commNum)
end
[a,b]=system('ls -d /danl/Harmon_dynCon/7*/Rest/Rest2.feat/');
c=sort(strread(b,'%s'));
commNum=zeros(size(c,1), 6);
for i=1:size(c)
load(strcat(char(c(i)),'/S_tmp.mat'))
commNum(i,:)=mean(max(S_tmp),3)
csvwrite('AdolesCommNumRest2.csv', commNum)
end
%Adults
cd /danl/Harmon_dynCon
[a,b]=system('ls -d /danl/Harmon_dynCon/4*/Rest/Rest1.feat/');
c=sort(strread(b,'%s'));
commNum=zeros(size(c,1), 6);
for i=1:size(c)
load(strcat(char(c(i)),'/S_tmp.mat'))
commNum(i,:)=mean(max(S_tmp),3)
csvwrite('AdultsCommNumRest1.csv', commNum)
end
[a,b]=system('ls -d /danl/Harmon_dynCon/4*/Rest/Rest2.feat/');
c=sort(strread(b,'%s'));
commNum=zeros(size(c,1), 6);
for i=1:size(c)
load(strcat(char(c(i)),'/S_tmp.mat'))
commNum(i,:)=mean(max(S_tmp),3)
csvwrite('AdultsCommNumRest2.csv', commNum)
end
AdultsCommRest <- read.csv("AdultsCommNumRest.csv")
To add back in the missing columns/rows for subject 713 before pulling all subj promiscuity stats #Should be updated for a_mat.mat & flex.mat
%For promiscuity
load('/danl/Harmon_dynCon/713/prom.mat')
promTEMP=prom
missing=[13 61; 61 13]
for i=1:length(missing)
%FOR inserting rows
len=size(promTEMP,i);
A = [promTEMP(1:missing(i)-1,:)];
B = NaN(1,4);
C = [promTEMP(missing(i):end,:)];
promTEMP=[A; B; C];
end
save('/danl/Harmon_dynCon/713/prom', 'promTEMP')
prom=promTEMP
save('/danl/Harmon_dynCon/713/prom', 'prom')
%For allegiance
load('/danl/Harmon_dynCon/713/a_mat_rerun1.mat')
a_matTEMP=a_mat;
a_matTEST=zeros(110,110,32);
missing=[13 61; 61 13];
for k=1:size(a_mat,3)
a_matTEMP=a_mat(:,:,k)
missing=[13 61; 61 13]
for i=1:length(missing)
%For inserting rows
len=size(a_matTEMP,i);
A = [a_matTEMP(1:missing(i)-1,:)]';
B = NaN(1,len)';
C = [a_matTEMP(missing(i):end,:)]';
a_matTEMP1=[A B C];
%For inserting columns
len=size(a_matTEMP1,2)
A = [a_matTEMP1(1:missing(i)-1,:)];
B = NaN(1,len);
C = [a_matTEMP1(missing(i):end,:)];
a_matTEMP=[A; B; C];
end
a_matTEST(:,:,k)=a_matTEMP;
end
save('/danl/Harmon_dynCon/713/a_mat_rerun1', 'a_matTEST')
a_mat=a_matTEST
save('/danl/Harmon_dynCon/713/a_mat_rerun1', 'a_mat')
%For flexibility
load('/danl/Harmon_dynCon/713/flex_rerun1.mat') %Renamed to '/danl/Harmon_dynCon/713/flex_rerun1_108.mat'
flexTEMP=flex
missing=[13 61; 61 13]
for i=1:length(missing)
%FOR inserting rows
len=size(flexTEMP,i);
A = [flexTEMP(1:missing(i)-1,:)];
B = NaN(1,4);
C = [flexTEMP(missing(i):end,:)];
flexTEMP=[A; B; C];
end
save('/danl/Harmon_dynCon/713/flex_rerun1', 'flexTEMP')
flex=flexTEMP;
save('/danl/Harmon_dynCon/713/flex_rerun1', 'flex')
For plotting and preparing for heirarchical models. Matlab.
forTask
%ADOLESCENTS
%load data and concatenate flexibility statistics
%[a,b]=system('ls -d /danl/Harmon_dynCon/7*/flex.mat');
[a,b]=system('ls -d /danl/Harmon_dynCon/7*/flex_rerun1.mat');
%do separately above and below
c=sort(strread(b,'%s'));
flex_cat=[];
for j=1:size(c,1)
load(char(c(j)))
flex_cat=cat(3,flex_cat,flex)
end
plot(squeeze(mean(flex_cat)))
block=repmat([1:4]',25,1);
sub=repmat([1:25]',1,4)'
sub=sub(:);
%reshape whole-brain average flexibility
meanflex=squeeze(mean(flex_cat));
meanflex=meanflex(:);
%get striatal average flexibility
%check to make sure indices are correct
[trash,roi_names]=system('ls /danl/Harmon_dynCon/Harvard-Oxford_ROIs/*nii.gz | xargs -n1 basename');
roi_names=strread(roi_names,'%s');
%str_ind=[49,51,54,104,106,109]; %FOR adults only
str_ind=[6,8,11,66,68,71];
roi_names(str_ind)
roi_names(hipp_ind)
strflex=squeeze(mean(flex_cat(str_ind,:,:)));
strflex=strflex(:);
plot(squeeze(mean(flex_cat(str_ind,:,:))))
flex_allrois=[];
for i=1:size(flex_cat,3)
flex_allrois=[flex_allrois;flex_cat(:,:,i)'];
end
%dlmwrite('/danl/Harmon_dynCon/flex_allrois.csv',flex_allrois)
dlmwrite('/danl/Harmon_dynCon/flex_allrois03.04.2019.csv',flex_allrois) %This is flex_rerun1
forRest
%Adolescents
%load data and concatenate flexibility statistics
%[a,b]=system('ls -d /danl/Harmon_dynCon/7*/Rest/Rest1.feat/flex.mat');
[a,b]=system('ls -d /danl/Harmon_dynCon/7*/Rest/Rest1.feat/flex_rerun1.mat');
%do separately above and below
c=sort(strread(b,'%s'));
flex_cat=[];
for j=1:size(c,1)
load(char(c(j)))
flex_cat=cat(3,flex_cat,flex);
end
%get flexibility scores for each ROI for each run for whole-brain search
%can prob do this more effeciently but this is easier to see, harder to botch
flex_allrois=[];
for i=1:size(flex_cat,3)
flex_allrois=[flex_allrois;flex_cat(:,:,i)'];
end
dlmwrite('/danl/Harmon_dynCon/flex_allroisRest1.03.04.2019.csv',flex_allrois)
%Rest 1 and Rest 2 are seperate
%load data and concatenate flexibility statistics
%[a,b]=system('ls -d /danl/Harmon_dynCon/7*/Rest/Rest2.feat/flex.mat');
[a,b]=system('ls -d /danl/Harmon_dynCon/7*/Rest/Rest2.feat/flex_rerun1.mat');
%do separately above and below
c=sort(strread(b,'%s'));
flex_cat=[];
for j=1:size(c,1)
load(char(c(j)))
flex_cat=cat(3,flex_cat,flex);
end
%get flexibility scores for each ROI for each run for whole-brain search
%can prob do this more effeciently but this is easier to see, harder to botch
flex_allrois=[];
for i=1:size(flex_cat,3)
flex_allrois=[flex_allrois;flex_cat(:,:,i)'];
end
dlmwrite('/danl/Harmon_dynCon/flex_allroisRest2.03.04.2019.csv',flex_allrois)
%Adults
%load data and concatenate flexibility statistics
%[a,b]=system('ls -d /danl/Harmon_dynCon/4*/Rest/Rest1.feat/flex.mat');
[a,b]=system('ls -d /danl/Harmon_dynCon/4*/Rest/Rest1.feat/flex_rerun1.mat');
%do separately above and below
c=sort(strread(b,'%s'));
flex_cat=[];
for j=1:size(c,1)
load(char(c(j)))
flex_cat=cat(3,flex_cat,flex);
end
%get flexibility scores for each ROI for each run for whole-brain search
%can prob do this more effeciently but this is easier to see, harder to botch
flex_allrois=[];
for i=1:size(flex_cat,3)
flex_allrois=[flex_allrois;flex_cat(:,:,i)'];
end
dlmwrite('/danl/Harmon_dynCon/flex_allroisRest1Adults.03.04.2019.csv',flex_allrois)
%Rest 1 and Rest 2 are seperate
%load data and concatenate flexibility statistics
%[a,b]=system('ls -d /danl/Harmon_dynCon/4*/Rest/Rest2.feat/flex.mat');
[a,b]=system('ls -d /danl/Harmon_dynCon/4*/Rest/Rest2.feat/flex_rerun1.mat');
%do separately above and below
c=sort(strread(b,'%s'));
flex_cat=[];
for j=1:size(c,1)
load(char(c(j)))
flex_cat=cat(3,flex_cat,flex);
end
%get flexibility scores for each ROI for each run for whole-brain search
%can prob do this more effeciently but this is easier to see, harder to botch
flex_allrois=[];
for i=1:size(flex_cat,3)
flex_allrois=[flex_allrois;flex_cat(:,:,i)'];
end
dlmwrite('/danl/Harmon_dynCon/flex_allroisRest2Adults.03.04.2019.csv',flex_allrois)
For plotting and preparing for heirarchical models. Matlab.
forTask
%load data and concatenate flexibility statistics
[a,b]=system('ls -d /danl/Harmon_dynCon/7*/prom.mat');
%do separately above and below
c=sort(strread(b,'%s'));
prom_cat=[];
for j=1:size(c,1)
load(char(c(j)))
prom_cat=cat(3,prom_cat,prom)
end
plot(squeeze(mean(prom_cat)))
saveas(gcf,'meanProm25.png')
saveas(gcf,'meanProm25.fig')
block=repmat([1:4]',25,1);
sub=repmat([1:25]',1,4)'
sub=sub(:);
%reshape whole-brain average flexibility
meanprom=squeeze(mean(prom_cat));
meanprom=meanprom(:);
%get striatal average flexibility
%check to make sure indices are correct
[trash,roi_names]=system('ls /danl/Harmon_dynCon/Harvard-Oxford_ROIs/*nii.gz | xargs -n1 basename');
%Do this seperately
roi_names=strread(roi_names,'%s');
%str_ind=[49,51,54,104,106,109]; %For adults only
str_ind=[6,8,11,66,68,71]; %For adoles only
roi_names(str_ind)
strprom=squeeze(mean(prom_cat(str_ind,:,:)));
strprom=strflex(:);
plot(squeeze(mean(prom_cat(str_ind,:,:))))
saveas(gcf,'strProm25.png')
saveas(gcf,'strProm25.fig')
%Hippocampus
hipp_ind=[9, 69]; %For adoels; When sorted use: hipp_ind=[52, 107];
roi_names(hipp_ind)
%roi_names(hipp_ind)
left_hipp=[9];
right_hipp=[69];
hipprom=squeeze(mean(prom_cat(hipp_ind,:,:)));
hipprom=hipprom(:);
hippromR=squeeze(mean(prom_cat(right_hipp,:,:)));
hippromR=hippromR(:);
hippromL=squeeze(mean(prom_cat(right_hipp,:,:)));
hippromL=hippromL(:);
plot(squeeze((prom_cat(right_hipp,:,:))))
plot(squeeze((prom_cat(left_hipp,:,:))))
plot(squeeze(mean(prom_cat(hipp_ind,:,:))))
saveas(gcf,'hippProm25.png')
saveas(gcf,'hippProm25.fig')
%write out csv for modeling in R
promdata=[sub block meanprom strprom]
dlmwrite('/danl/Harmon_dynCon/promdata.csv',promdata)
%get promiscuity scores for each ROI for each run for whole-brain search
%can prob do this more effeciently but this is easier to see, harder to botch
prom_allrois=[];
for i=1:size(prom_cat,3)
prom_allrois=[prom_allrois;prom_cat(:,:,i)'];
end
dlmwrite('/danl/Harmon_dynCon/prom_allrois.csv',prom_allrois) forRest
%Adolescents
%load data and concatenate flexibility statistics
[a,b]=system('ls -d /danl/Harmon_dynCon/7*/Rest/Rest1.feat/prom.mat');
%do separately above and below
c=sort(strread(b,'%s'));
%get flexibility scores for each ROI for each run for whole-brain search
%can prob do this more effeciently but this is easier to see, harder to botch
prom_allrois=[];
prom_cat=[];
for i=1:size(prom_cat,3)
prom_allrois=[prom_allrois;prom_cat(:,:,i)’];
end
dlmwrite('/danl/Harmon_dynCon/prom_allroisRest1.03.04.2019.csv',prom_allrois)
%Rest2
%load data and concatenate flexibility statistics
[a,b]=system('ls -d /danl/Harmon_dynCon/7*/Rest/Rest2.feat/prom.mat');
%do separately above and below
c=sort(strread(b,'%s'));
%get flexibility scores for each ROI for each run for whole-brain search
%can prob do this more effeciently but this is easier to see, harder to botch
prom_allrois=[];
prom_cat=[];
for i=1:size(prom_cat,3)
prom_allrois=[prom_allrois;prom_cat(:,:,i)’];
end
dlmwrite('/danl/Harmon_dynCon/prom_allroisRest2.03.04.2019.csv',prom_allrois)
%Adults
%load data and concatenate flexibility statistics
[a,b]=system('ls -d /danl/Harmon_dynCon/4*/Rest/Rest1.feat/prom.mat');
%do separately above and below
c=sort(strread(b,'%s'));
%get flexibility scores for each ROI for each run for whole-brain search
%can prob do this more effeciently but this is easier to see, harder to botch
prom_allrois=[];
prom_cat=[];
for i=1:size(prom_cat,3)
prom_allrois=[prom_allrois;prom_cat(:,:,i)’];
end
dlmwrite('/danl/Harmon_dynCon/prom_allroisRest1Adults.03.04.2019.csv',prom_allrois)
%Rest2
%load data and concatenate flexibility statistics
[a,b]=system('ls -d /danl/Harmon_dynCon/4*/Rest/Rest2.feat/prom.mat');
%do separately above and below
c=sort(strread(b,'%s'));
%get flexibility scores for each ROI for each run for whole-brain search
%can prob do this more effeciently but this is easier to see, harder to botch
prom_allrois=[];
prom_cat=[];
for i=1:size(prom_cat,3)
prom_allrois=[prom_allrois;prom_cat(:,:,i)’];
end
dlmwrite('/danl/Harmon_dynCon/prom_allroisRest2Adults.03.04.2019.csv',prom_allrois)
forTask
%load data and concatenate flexibility statistics
[a,b]=system('ls -d /danl/Harmon_dynCon/4*/prom.mat');
%do separately above and below
c=sort(strread(b,'%s'));
prom_cat=[];
for j=1:size(c,1)
load(char(c(j)))
prom_cat=cat(3,prom_cat,prom)
end
plot(squeeze(mean(prom_cat)))
saveas(gcf,'meanProm22.png')
saveas(gcf,'meanProm22.fig')
block=repmat([1:4]',22,1);
sub=repmat([1:22]',1,4)'
sub=sub(:);
%reshape whole-brain average flexibility
meanprom=squeeze(mean(prom_cat));
meanprom=meanprom(:);
%get striatal average flexibility
%check to make sure indices are correct
[trash,roi_names]=system('ls /danl/Harmon_dynCon/Harvard-Oxford_ROIs/*nii.gz | xargs -n1 basename'); %SHOULD USE ORDER FROM CONN_CELLS different from them in this list
roi_names=sort(strread(roi_names,'%s'));
str_ind=[49,51,54,104,106,109];
roi_names(str_ind)
strprom=squeeze(mean(prom_cat(str_ind,:,:)));
strprom=strprom(:);
plot(squeeze(mean(prom_cat(str_ind,:,:))))
saveas(gcf,'strProm22.png')
saveas(gcf,'strProm22.fig')
%Hippocampus
%hipp_ind=[9, 69]; When using sort:
hipp_ind=[52, 107]
roi_names(hipp_ind)
%roi_names(hipp_ind)
left_hipp=[52];
right_hipp=[107];
hipprom=squeeze(mean(prom_cat(hipp_ind,:,:)));
hipprom=hipprom(:);
hippromR=squeeze(mean(prom_cat(right_hipp,:,:)));
hippromR=hippromR(:);
hippromL=squeeze(mean(prom_cat(right_hipp,:,:)));
hippromL=hippromL(:);
plot(squeeze((prom_cat(right_hipp,:,:))))
plot(squeeze((prom_cat(left_hipp,:,:))))
plot(squeeze(mean(prom_cat(hipp_ind,:,:))))
saveas(gcf,'hippProm22.png')
saveas(gcf,'hippProm22.fig')
%write out csv for modeling in R
promdata=[sub block meanprom strprom]
dlmwrite('/danl/Harmon_dynCon/promdata_adults.csv',promdata)
%get promiscuity scores for each ROI for each run for whole-brain search
%can prob do this more effeciently but this is easier to see, harder to botch
prom_allrois=[];
for i=1:size(prom_cat,3)
prom_allrois=[prom_allrois;prom_cat(:,:,i)'];
end
dlmwrite('/danl/Harmon_dynCon/prom_allrois_adults.csv',prom_allrois) a=dir('/danl/Harmon_dynCon/behavior/7*_tb_l*_svlo*.mat');
a={a.name};
longform=[];
for k = 1:length(a)
k
clear subnum chose_right shown_corr trial
filea=a{(k)};
load(a{(k)});
subnum=str2double(num2str(a{k}(1:3)))';
chose_right=double(strcmp(resp,'b'))';
shown_corr=double(shown_corr)';
no_resp=shown_corr==2;
chose_right(no_resp)=NaN;
shown_corr(no_resp)=NaN;
ntrials=size(resp',1);
trial=1:ntrials';
longform=[longform;repmat(subnum,ntrials,1) stim_shown' trial' chose_right shown_corr];
end
header={'sub','stim','trial','choice','fb'};
header2=sprintf('%s,',header{:});header2(end)=[];
dlmwrite('/danl/Harmon_dynCon/behavior/choice_fb_long.csv',...
header2,'')
dlmwrite('/danl/Harmon_dynCon/behavior/choice_fb_long.csv',...
longform,'-append','delimiter',',')
#do not run
!cp /danl/Harmon_dynCon/behavior/choice_fb_long.csv ~/GitHub/rl_flexibility/RL_model/choice_fb_long.csv#This code is also found in ~/Columbia/learninglab/dynCon/ch_scripts/RL_model/heirarBeysModel.R
library(rstan)
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())
dat <- read.csv('~/Columbia/learninglab/dynCon/ch_scripts/RL_model/choice_fb_long.csv')
choices<-unique(na.omit(dat$choice))
dat$chosen<-dat$choice
dat$chosen[dat$choice==choices[1]]=1
dat$chosen[dat$choice==choices[2]]=2
dat$unchosen=abs(dat$chosen-3)
dat$chose_two<-dat$chosen==2
subs = unique(dat$sub);
NS = length(subs);
NStim=length(unique(dat$stim));
MT=max(dat$trial);
NT = array(0,NS);
choice = array(0,c(NS,MT));
unchoice=choice;
choice_two=choice;
rew = choice;
stim=choice;
for (i in 1:NS) {
NT[i] = nrow(subset(dat,sub==subs[i]));
stim[i,1:NT[i]] = subset(dat,sub==subs[i])$stim;
#choice and reward history
choice[i,1:NT[i]] = subset(dat,sub==subs[i])$chosen;
unchoice[i,1:NT[i]] = subset(dat,sub==subs[i])$unchosen;
rew[i,1:NT[i]] = subset(dat,sub==subs[i])$fb;
#based on choosing second option
choice_two[i,1:NT[i]] = subset(dat,sub==subs[i])$chose_two;
}
choice[is.na(choice)]<--1
unchoice[is.na(unchoice)]<--1
rew[is.na(rew)]<--1
choice_two[is.na(choice_two)]<--1
rl_standata = list(NS=NS, NC=2, NStim=NStim, MT=MT, NT= NT, choice=choice,
stim=stim,choice_two=choice_two,rew=rew );
rl_fit <- stan(file = '~/Columbia/learninglab/dynCon/ch_scripts/RL_model/standard_rl.stan',
data = rl_standata, iter = 1250, warmup = 250, chains = 4)
save(rl_fit,file='~/Columbia/learninglab/dynCon/ch_scripts/RL_model/rl_fit')#Code can be found in ~/Columbia/learninglab/dynCon/ch_scripts/RL_model/extraxtParamsFromModel.R
library(rstan)
#install.packages('ppcor')
library(ppcor)
fit_rl<-load('~/Columbia/learninglab/dynCon/ch_scripts/RL_model/rl_fit')
fit_extract<-extract(rl_fit,permute=T)
#fit_extract$beta_mean<-fit_extract$b1/fit_extract$b2 #<- dividing is wrong #Group distribution - beta estimate mean for each posterior sample (stan gamma distribution) **Should be multiplied not divided!! because inverse. STAN uses rate - invesre of scale - rather than scale. THe inverse mean formula matters based on inverse
#beta = inverse temperature
hist(fit_extract$beta_mean) #model uncertainty about the mean
#could also calculate SD
fit_extract$alpha_mean<-fit_extract$a1/(fit_extract$a1+fit_extract$a2)
#alpha learning rate
hist(fit_extract$alpha_mean)
#variance
hist(fit_extract$b1*fit_extract$b2^2)
#When parameters dont change in the paradigm- inverese temp is similar to percent correct - lower alpha better learning, higher inverse temp, better learning
#Check with Juliet's paper similar values
#with unstable estimes - first thing to try is to make hyperparameters - non informed priors - may be worth shrinking the dist. = making priors find more extreme values less plausble - and if the models is still fitting well then maybe thats meaningul -keep priors the same if comparing the two
#beta how much does value influence decision
betas<-apply(fit_extract$beta,2,mean)
alphas<-apply(fit_extract$alpha,2,mean)
for (i in 1:25){
plot(fit_extract$alpha[,i], fit_extract$beta[,i], xlab="alpha", ylab='beta', main=paste('Subj', i))
}
#Check out Gamma distribution calculations for fits -
#https://en.wikipedia.org/wiki/Gamma_distribution
fit_extract$beta_mean<-fit_extract$b1*fit_extract$b2
fit_extract$alpha_mean<-fit_extract$a1/(fit_extract$a1+fit_extract$a2)
str_acorrel<-matrix(0,4000,6)
str_bcorrel<-matrix(0,4000,6)
str_apbcorrel<-matrix(0,4000,6)
str_bpacorrel<-matrix(0,4000,6)
wb_acorrel<-NULL
wb_bcorrel<-NULL
wb_apbcorrel<-NULL
wb_bpacorrel<-NULL
flexdat<-read.csv("flex_allrois.csv",header=F)
flex_behav<-read.csv("flex_behav.csv",header=T)
View(flex_behav)
flexdat$Sub<-rep(seq(1,25,1),each=4)
flexdat$Block<-rep(seq(1,4,1),times=25)
flexdat$Corr<-flex_behav$correct
flexdat$weights<-flex_behav$weights
flexdat<-melt(flexdat,id=c("Sub","Block","Corr","weights"))
names(flexdat)[c(5,6)]<-c("ROI","flex")
flexdat$ROI<-as.factor(as.numeric(flexdat$ROI))
meanflex_rois<-tapply(flexdat$flex,list(flexdat$Sub,flexdat$ROI),mean)
meanflex<-rowMeans(meanflex_rois)
str_ind=c(49,51,54,104,106,109);
for( i in seq(1,dim(fit_extract$beta)[1],1)){
betas_tmp=fit_extract$beta[i,]
alphas_tmp=fit_extract$alpha[i,]
k<-1
wb_acorrel[i]<-cor(meanflex,alphas_tmp)
wb_bcorrel[i]<-cor(meanflex,betas_tmp)
wb_apbcorrel[i]<-pcor.test(meanflex,alphas_tmp,betas_tmp)[1] #partial correlations
wb_bpacorrel[i]<-pcor.test(meanflex,betas_tmp,alphas_tmp)[1] #partial correlations
for(j in str_ind){
str_acorrel[i,k]<-cor(meanflex_rois[,k],alphas_tmp)
str_bcorrel[i,k]<-cor(meanflex_rois[,k],betas_tmp)
str_bpacorrel[i,k]<-pcor.test(meanflex_rois[,k],betas_tmp,alphas_tmp)$estimate
str_apbcorrel[i,k]<-pcor.test(meanflex_rois[,k],alphas_tmp,betas_tmp)$estimate
k<-k+1
}
}
par(mfrow=c(1,2))
plot(rowMeans(str_bcorrel),rowMeans(str_acorrel),
col=rgb(0,0,0,alpha=0),pch=21,bg=rgb(0,0,0,alpha=.03),
xlab="Correlation with Beta",ylab="Correlation with Alpha",
main="Striatum Flexibility",xlim=c(-.1,.5))
abline(v=0,lty=2)
abline(h=0,lty=2)
plot(wb_bcorrel,wb_acorrel,
col=rgb(0,0,0,alpha=0),pch=21,bg=rgb(0,0,0,alpha=.03),
xlab="Correlation with Beta",ylab="Correlation with Alpha",
main="Whole-brain Flexibility",xlim=c(-.1,.5))
abline(v=0,lty=2)
abline(h=0,lty=2)
par(mfrow=c(1,1))
plot(fit_extract$alpha[,1], fit_extract$beta[,1], xlab="alpha", ylab='beta', main='Subj 1')
plot(fit_extract$alpha[,2], fit_extract$beta[,2], xlab="alpha", ylab='beta', main='Subj 2')
plot(fit_extract$alpha[,3], fit_extract$beta[,3], xlab="alpha", ylab='beta', main='Subj 3')
plot(fit_extract$alpha[,1], fit_extract$beta[,1], xlab="alpha", ylab='beta', main='Subj 1')
Estimate the effect of striatal and whole-brain flexibility on reinforcement learning. See models.Rmd and models.pdf for more details. R
#This code can also be found ~/Columbia/learninglab/dynCon/ch_scripts/RL_model/MLBayesianHeirarchicalModels.R
#prepare data for binomial logistic modelling
library(lme4)
library(reshape2)
library(brms)
flexdat<-read.csv("flex_allrois.csv",header=F)
#read in trial-by-by trial behavioral data
data<-read.delim('behav_data.tsv',header=1)
data<-read.csv('behav_data.csv',header=1)
data <- data[data$group==7,]
#Only have behavioral data : no brain data
data <- data[!data$subjectNum==716,]
data <- data[!data$subjectNum==718,]
data$block<-rep(rep(seq(1,4,1),each=30),25)
flexdat$subject<-rep(seq(1,25,1),each=4)
flexdat$Block<-rep(seq(1,4,1),times=25)
flexdat$Corr<-melt(tapply(data$optCor,list(data$block,data$subjectNum),mean,na.rm=1))$value
nacount<-is.na(data$optCor)
weights<-melt(tapply(nacount,list(data$block,data$subjectNum),sum))
flexdat$weights<-30-weights[,3]
str_ind=c(49,51,54,104,106,109);
flexdat<-melt(flexdat,id=c("subject","Block","Corr","weights"))
names(flexdat)[c(5,6)]<-c("ROI","flex")
flexdat$ROI<-as.factor(as.numeric(flexdat$ROI))
flexdat_str<-subset(flexdat,ROI %in% str_ind)
flexdat_str$ROI<-as.factor(flexdat_str$ROI)
flexdat_str_avg<-melt(tapply(flexdat_str$flex,list(flexdat_str$subject,flexdat_str$ROI),mean,na.rm=1))
names(flexdat_str_avg)<-c("subject","ROI","flex")
flexdat_str_avg<-subset(flexdat_str_avg,ROI %in% str_ind)
flexdat_str$meanflex<-rep(flexdat_str_avg$flex,each=4)
flexdat_str_melt<-melt(tapply(flexdat_str$flex,list(flexdat_str$subject,flexdat_str$Block),mean))
names(flexdat_str_melt)<-c("subject","Block","flex")
flexdat_str_melt$correct<-melt(tapply(flexdat_str$Corr,list(flexdat_str$subject,flexdat_str$Block),mean))$value
flexdat_str_melt$weights<-melt(tapply(flexdat_str$weights,list(flexdat_str$subject,flexdat_str$Block),mean))$value
m_str<-glmer(correct~flex+(flex||subject),
family=binomial,weights=weights,
data=flexdat_str_melt)
write.csv(flexdat_str_melt,'flex_behav.csv')
#for posterior inference, run bayesian model using brms wrapper for stan
options(mc.cores = parallel::detectCores())
flexdat_str_melt$numcorr<-as.integer(flexdat_str_melt$correct*flexdat_str_melt$weights)
mlearn_stan<-brm(numcorr~flex+(flex|subject),data=flexdat_str_melt,family=binomial)
##
Warning messages:
1: There were 15 divergent transitions after warmup. Increasing adapt_delta above 0.8 may help. See
http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
2: Examine the pairs() plot to diagnose sampling problems
R
#This code can also be found ~/Columbia/learninglab/dynCon/ch_scripts/RL_model/wholeBrainSearchforEffectsofFlexonLearning.R
#load in data
library(lme4)
roi_data<-read.csv('flex_allrois.csv',header=0)
behav<-read.csv('flex_behav.csv')
p<-NULL
b<-NULL
for(i in 1:dim(roi_data)[2]){
mtmp<-summary(glmer(behav$correct~roi_data[,i]+(roi_data[,i]||behav$subject),family=binomial,weights=behav$weights))
if(dim(summary(mtmp)[10]$coefficients)[2]>=4){
p[i]<-summary(mtmp)[10]$coefficients[2,4]
b[i]<-summary(mtmp)$coefficients[2,1]
} else{
p[i]<-NA
b[i]<-NA
}
}
write(p,'p_vals_learning_glm.csv')
write(b,'b_vals_learning_glm.csv')
##
Warning message:
In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
Model failed to converge with max|grad| = 0.00105927 (tol = 0.001, component 1)
Matlab
%FOR ADULTS
str_ind=[49,51,54,104,106,109];
[a,b]=system('ls -d /data/engine/rgerraty/learn_dyncon/4*/a_mat.mat');
%FOR ADOLESCENTS
str_ind=[6,8,11,66,68,71];
[a,b]=system('ls -d /danl/Harmon_dynCon/7*/a_mat_rerun1.mat');
c=sort(strread(b,'%s'));
h=1
for s=1:length(c)
clear a_mat clear a_mat_str
load(c{s})
a_mat(a_mat==1)=NaN;
a_mat_str=a_mat(:,str_ind,:);
for i=1:size(a_mat_str,1)
for j=1:size(a_mat_str,2)
for k=1:size(a_mat_str,3)
a_mat_long(h,:)=[a_mat_str(i,j,k),i,j,ceil(k/8),s];
h=h+1;
end
end
end
end
header={'allegiance','ROI','str_roi','block','sub'};
header2=sprintf('%s,',header{:});header2(end)=[];
%ADULTS
dlmwrite('/data/engine/rgerraty/learn_dyncon/alleg_long.csv',...
header2,'')
dlmwrite('/data/engine/rgerraty/learn_dyncon/alleg_long.csv',...
a_mat_long,'-append','delimiter',',')
%ADOLESCENTS
dlmwrite('/danl/Harmon_dynCon/alleg_long03.04.2019.csv',...
header2,'')
dlmwrite('/danl/Harmon_dynCon/alleg_long03.04.2019.csv',...
a_mat_long,'-append','delimiter',',')
Adolescents
%read in all subject/run mc motion files
[a,b]=system('ls -d /danl/Harmon_dynCon/7*/Learn?_PEpriorD.feat/mc/prefiltered_func_data_mcf_rel_mean.rms');
%do separately above and below
c=sort(strread(b,'%s'))
motion = (1:length(c))'
for i=1:length(c)
load(c{i})
motion(i) = prefiltered_func_data_mcf_rel_mean
end
dlmwrite('/danl/Harmon_dynCon/motion.csv', motion)Adults
%read in all subject/run mc motion files
[a,b]=system('ls -d /data/engine/rgerraty/learn_dyncon/4*/Learn?_PEprior.feat/mc/prefiltered_func_data_mcf_rel_mean.rms');
%do separately above and below
c=sort(strread(b,'%s'))
motion = (1:length(c))'
for i=1:length(c)
load(c{i})
motion(i) = prefiltered_func_data_mcf_rel_mean
end
dlmwrite('/data/engine/charmon/motion.csv', motion)