%% This script shows how to compute the framework's coefficients, and the % ensuing results, for an example simulation. The paper's content is % essentially re-using similar code across simulation sub-types, across % methods that retrieve a coefficient set or another, and also for % experimental data investigations %% 1. Loading of the example data and addition of the needed paths to the % framework functions and plotting utilities load('ExampleData.mat'); addpath('./LogisticRegressionFramework/Framework'); addpath('./LogisticRegressionFramework/PostProcessing'); addpath(genpath('./LogisticRegressionFramework/Plotting')); % The loaded data contains: % - The ground truth parameters for co-activations and causal couplings, % respectively in Gamma_GT_regions (R x R matrix with network % assignments)/NetworkAssignment_GT (R x 1 vector with same information), % and in Beta_GT_networks (N x N causal modulations across % networks)/Beta_GT_regions (same in R x R format, across regions) % - The example simulated data, for the training set of S = 50 subjects (in % TC) and the cross-validation set of S = 30 subjects (in CV) % - The number of regions R = 35 and number of networks N = 7 %% 2. Binarizes the time courses TC_bin = RSCHMM_Binarize(TC); CV_bin = RSCHMM_Binarize(CV); clear TC clear CV %% 3. Selects the time points % This specifies the size of each individual structure containing data % points; it does not impact the outputs from the framework, but it plays a % role in computational speed. 5000 was empirically selected as a % satisfying value limit = 5000; % Selection of the time points: X_coact and X_causal are the activity % levels of the "non r" regions at times t+1 and t, respectively, while % Y denotes whether region r has transited or not at time t [X_coact,X_causal,Y] = RSCHMM_Select_WorkingSet_Optimized(TC_bin,limit); [X_coact_CV,X_causal_CV,Y_CV] = RSCHMM_Select_WorkingSet_Optimized(CV_bin,limit); clear TC_bin clear CV_bin %% 4. Computes the coefficients % Tolerance threshold used for convergence epsilon = 1e-2; % Number of iterations n_iter = 5; % Currently unused, but EN enables to perform elastic net regularisation % instead of purely L1 regularisation-based solutions as now. Keep EN = 1 % not to use elastic net EN = 1; % Range of probed lambda values for simulations lambda_range = [10000,9000,8000,7000,6000,5000,4600,4200,3800,... 3200,2800,2600,2400,2200,2000,1800,1600,1400,1200,1000,... 900,800,700,650,600,550,500,450,400,350,300,250,200,180,160,140,120,100,90,80,70,60,... 50,40,30,25,20,18,16,14,12,10,9.6,9.2,8.8,8.4,8,7.6,7.2,6.8,6.4,6,... 5.6,5.2,4.8,4.4,4,3.6,3.2,2.8,2.4,2,1.8,1.6,1.2,0.8,0.4,... 0.2,0.1,0]; % The commented range below is that used in experimental data applications % in the paper. It started at larger values because otherwise, % max(lambda_range) did not yield only null coefficients % lambda_range = [50000,48000,46000,44000,42000,40000,38000,36000,34000,... % 32000,30000,28000,26000,24000,22000,20000,18000,16000,14000,13000,... % 12000,11000,10000,9000,8600,8200,7800,7400,7000,6800,6600,6400,6200,... % 6000,5800,5600,5400,5200,5000,4800,4600,4400,4200,4000,3800,3600,... % 3400,3200,3000,2800,2600,2400,2200,2000,1800,1600,1400,1200,1000,... % 900,800,700,600,500,400,300,200,180,160,140,120,100,90,80,70,60,... % 50,40,30,20,18,16,14,12,10,9.6,9.2,8.8,8.4,8,7.6,7.2,6.8,6.4,6,... % 5.6,5.2,4.8,4.4,4,3.6,3.2,2.8,2.4,2,1.8,1.6,1.4,1.2,1,0.8,0.6,0.4,... % 0.2,0.1,0.06,0.02,0]; % lambda_range = downsample(lambda_range,2); % lambda_range = [lambda_range,0]; n_lambda = length(lambda_range); % Declaration of the coefficient matrices and times for each framework step Alpha = NaN(n_lambda,2,R,5); Beta = NaN(n_lambda,R,2,R,5); Gamma = NaN(n_lambda,R,2,R,5); t1 = NaN(n_lambda,R,5); t2 = NaN(n_lambda,R,5); t3 = NaN(n_lambda,R,5); % Solving per se to retrieve the coefficients. Note that parallel % computations are run across regions with the "parfor" syntax (remove if % you do not have the toolbox, and replace with a simple for loop) idx_rho = 1; % Five values for rho are probed only in the paper; it could be interesting % to examine a larger range, but it will make computations longer for rho = [0,0.25,0.5,0.75,1] parfor r = 1:R [Alpha(:,:,r,idx_rho),Beta(:,:,:,r,idx_rho),Gamma(:,:,:,r,idx_rho),... ~,t1(:,r,idx_rho),t2(:,r,idx_rho),t3(:,r,idx_rho)] = ... RSCHMM_Optimized_ROI_Deluxe(lambda_range,... epsilon,rho,EN,X_coact{r},X_causal{r},Y{r},r,n_iter); end idx_rho = idx_rho + 1; end clear idx_rho clear X_coact clear X_causal clear Y %% 5. Evaluates the log-likelihood on the cross-validation dataset % Initialises the counter for rho values iter_rho = 1; % Declares the log-likelihood matrix LL_RSCHMM = NaN(5,n_lambda,2,R); % Computes the log-likelihood across rho and lambda values on the % cross-validated data for rho = [0,0.25,0.5,0.75,1] iter_regul = 1; for lambda = lambda_range for r = 1:R for i = 1:2 LL_RSCHMM(iter_rho,iter_regul,i,r) = RSCHMM_ComputeLq_Optimized_Deluxe(Alpha(iter_regul,i,r,iter_rho),... Beta(iter_regul,(1:R)~=r,i,r,iter_rho)',Gamma(iter_regul,(1:R)~=r,i,r,iter_rho)',... X_causal_CV{r}(:,i),X_coact_CV{r}(:,i),Y_CV{r}(:,i)); end end iter_regul = iter_regul + 1; end iter_rho = iter_rho + 1; end clear iter_rho clear iter_regul clear X_causal_CV clear X_coact_CV clear Y_CV % Determines the optimal indices (where the log-likelihood is maximal) for % rho and for lambda idx_RSCHMM = NaN(2,R); idx1_RSCHMM = NaN(2,R); idx2_RSCHMM = NaN(2,R); for r = 1:R for i = 1:2 tmp_LL = squeeze(LL_RSCHMM(:,:,i,r)); [~,idx_RSCHMM(i,r)] = max(tmp_LL(:)); [idx1_RSCHMM(i,r),idx2_RSCHMM(i,r)] = ind2sub(size(squeeze(LL_RSCHMM(:,:,i,r))),idx_RSCHMM(i,r)); end end clear tmp_LL %% 6. Puts the coefficients into matrices as a function of coefficient % type (co-activation vs causal) and start state of activity for probed % region (baseline vs active) % Declares the coefficient matrices Beta_final_baseline = NaN(R,R); Beta_final_active = NaN(R,R); Gamma_final_baseline = NaN(R,R); Gamma_final_active = NaN(R,R); % Fills them for r = 1:R Gamma_final_baseline(:,r) = Gamma(idx2_RSCHMM(1,r),:,1,r,idx1_RSCHMM(1,r)); Beta_final_baseline(:,r) = Beta(idx2_RSCHMM(1,r),:,1,r,idx1_RSCHMM(1,r)); Gamma_final_active(:,r) = Gamma(idx2_RSCHMM(2,r),:,2,r,idx1_RSCHMM(2,r)); Beta_final_active(:,r) = Beta(idx2_RSCHMM(2,r),:,2,r,idx1_RSCHMM(2,r)); end % Also creates the same matrix for alpha: Alpha_final has baseline % coefficients as first column, and active ones as second column Alpha_final = RSCHMM_PickAlpha(Alpha,idx1_RSCHMM,idx2_RSCHMM); %% 7. Computes the resulting probabilistic couplings % Separate computations for the baseline and the active state; DP_gamma and % DP_beta outputs are R x R matrices [DP_gammaB,DP_betaB] = RSCHMM_ComputeDeltaP_Optimized(Alpha_final(:,1),Beta_final_baseline,Gamma_final_baseline); [DP_gammaA,DP_betaA] = RSCHMM_ComputeDeltaP_Optimized(Alpha_final(:,2),Beta_final_active,Gamma_final_active); % Subtraction of both state coefficients (used everywhere except in Figure % 5C) DP_gamma = DP_gammaB - DP_gammaA; DP_beta = DP_betaB - DP_betaA; %% 8. Computes the metrics of interest % From the co-activation matrix, determines network assignment assuming N % networks (hierarchical clustering through Ward's linkage analysis) [NetworkAssignment_SLR,Z_SLR] = RSCHMM_DetermineCommunities(DP_gamma,N); % Accuracy of the network assignments, knowing the ground truth number of % networks N Purity_SLR = compute_clustering_accuracy(NetworkAssignment_SLR,NetworkAssignment_GT,N); % Similarity of the retrieved co-activation coefficients with the ground % truth: we consider the upper diagonal and the lower diagonal together, % because the set of retrieved coefficients is actually not enforced to be % symmetric by the framework Sim_Gamma_SLR = corr([jUpperTriMatToVec(Gamma_GT_regions);jUpperTriMatToVec(Gamma_GT_regions')],... [jUpperTriMatToVec(DP_gamma);jUpperTriMatToVec((DP_gamma)')]); % Similarity of the retrieved causal coefficients with the ground truth, % computed in a similar manner Sim_Beta_SLR = corr([jUpperTriMatToVec(Beta_GT_regions);jUpperTriMatToVec(Beta_GT_regions')],... [jUpperTriMatToVec(DP_beta);jUpperTriMatToVec((DP_beta)')]); % Construction of the graph representation associated to the retrieved % causal outputs, knowing the ground truth networks [Graph_SLR] = RSCHMM_ConstructPopulationGraph(Beta_GT_regions,NetworkAssignment_GT,N); % Computation of sensitivity and specificity [Sensitivity_SLR,Specificity_SLR] = RSCHMM_ComputeSenSpec(Graph_SLR,Beta_GT_networks); %% 9. Plots example outputs from the framework % Red/Blue colormap CM = cbrewer('div','RdBu',1001); % Plot of the regional time courses at hand for 200 samples (on one % subject) figure; set(gcf,'color','w'); imagesc(TC{1}(:,1:200)); colormap(flipud(CM)); caxis([-6 6]); set(gca,'Box','off'); % The ground truth matrices figure; set(gcf,'color','w'); imagesc(Gamma_GT_regions); axis off axis square colormap(flipud(CM)); caxis([-1,1]); figure; set(gcf,'color','w'); imagesc(Beta_GT_regions); axis off axis square colormap(flipud(CM)); caxis([-1,1]); % The network-level ground truth graph, and the SLR output one RSCHMM_PlotTP(Beta_GT_networks); % The final solutions figure; set(gcf,'color','w'); imagesc(DP_gamma); axis off axis square colormap(flipud(CM)); caxis([-0.15,0.15]); figure; set(gcf,'color','w'); imagesc(DP_beta); axis off axis square colormap(flipud(CM)); caxis([-0.06,0.06]); % The SLR-derived graph RSCHMM_PlotTP(Graph_SLR);