% Example 2: Analyze different splits of runs in a session
0001 %% Example 2: Analyze different splits of runs in a session 0002 0003 %% Download dataset (if necessary) and add GLMdenoise to the MATLAB path 0004 0005 setup; 0006 0007 %% Load in the data 0008 0009 % Load in the data 0010 load('exampledataset.mat'); 0011 0012 % Check the workspace 0013 whos 0014 %% 0015 0016 %% Split session into two halves and analyze each half separately 0017 0018 % In some circumstances, you may want to split the runs acquired in a given 0019 % session into groups and analyze each group separately. For example, you might 0020 % want to train a classifier on data from some runs and then test the 0021 % classifer on data from other runs. Here we go through an example where 0022 % we split a session of data (10 runs) into halves and analyze each half separately. 0023 0024 % First, we define the way the runs are to be split. In the first 0025 % iteration, we will analyze the odd runs. In the second iteration, 0026 % we will analyze the even runs. 0027 splits = {[1:2:10] [2:2:10]}; 0028 0029 % For each split, we will analyze the data two ways. First, we analyze the data 0030 % using no denoising (REGULAR). Then, we analyze the data using denoising (DENOISE). 0031 % (To achieve no denoising, we set the number of PCs to try to zero.) 0032 clear resultsREGULAR resultsDENOISE; 0033 for p=1:length(splits) 0034 ix = splits{p}; 0035 0036 % analyze using no denoising 0037 resultsREGULAR(p) = GLMdenoisedata(design(ix),data(ix),stimdur,tr, ... 0038 [],[],struct('numpcstotry',0), ... 0039 sprintf('example2figures_REGULAR_split%d',p)); 0040 0041 % analyze using denoising 0042 resultsDENOISE(p) = GLMdenoisedata(design(ix),data(ix),stimdur,tr, ... 0043 [],[],[], ... 0044 sprintf('example2figures_DENOISE_split%d',p)); 0045 0046 end 0047 %% 0048 0049 % The default mode of operation involves estimating an HRF from the data (in each 0050 % call to GLMdenoisedata). Let's inspect the HRF results across the two splits of 0051 % the data. The reason this may be important is that the beta weights that are 0052 % estimated in the model are scale factors that are applied to the HRF estimate, 0053 % and if the HRF estimates are highly different across the splits, that may be some 0054 % reason for concern. In this case, we find that the HRF estimate is fairly 0055 % similar across the splits. (An alternative strategy is to just fix the HRF to 0056 % a canonical HRF and use that for each split.) 0057 figure; hold on; 0058 hrf1 = resultsREGULAR(1).modelmd{1}; 0059 hrf2 = resultsREGULAR(2).modelmd{1}; 0060 plot(0:tr:tr*(length(hrf1)-1),hrf1,'r-'); 0061 plot(0:tr:tr*(length(hrf2)-1),hrf2,'b-'); 0062 xlabel('Time from condition onset (s)'); 0063 ylabel('Response'); 0064 title('HRF estimate across splits (red = odd runs, blue = even runs)'); 0065 %% 0066 0067 % Let's now examine the stability of beta weights across the data splits. 0068 % First, we choose a subset of the voxels. 0069 ix = find(resultsREGULAR(1).R2 > 10 & resultsREGULAR(1).meanvol > 500); 0070 0071 % Then, we extract the beta weights corresponding to one of the 0072 % conditions (condition 10). Notice that we have four sets of beta weights. 0073 % This is because for each of the two splits of the data, we have two 0074 % beta weight estimates (one using no denoising, one using denoising). 0075 beta1REGULAR = subscript(squish(resultsREGULAR(1).modelmd{2},3),{ix 10}); 0076 beta2REGULAR = subscript(squish(resultsREGULAR(2).modelmd{2},3),{ix 10}); 0077 beta1DENOISE = subscript(squish(resultsDENOISE(1).modelmd{2},3),{ix 10}); 0078 beta2DENOISE = subscript(squish(resultsDENOISE(2).modelmd{2},3),{ix 10}); 0079 0080 % Visualize the results. 0081 figure; 0082 set(gcf,'Units','points','Position',[100 100 800 400]); 0083 subplot(1,2,1); hold on; 0084 scatter(beta1REGULAR,beta2REGULAR,'r.'); 0085 axissquarify; 0086 xlabel('BOLD signal (% change), odd runs'); 0087 ylabel('BOLD signal (% change), even runs'); 0088 title(sprintf('No denoising (r=%.3f)',corr(beta1REGULAR,beta2REGULAR))); 0089 subplot(1,2,2); hold on; 0090 scatter(beta1DENOISE,beta2DENOISE,'r.'); 0091 axissquarify; 0092 xlabel('BOLD signal (% change), odd runs'); 0093 ylabel('BOLD signal (% change), even runs'); 0094 title(sprintf('Denoising (r=%.3f)',corr(beta1DENOISE,beta2DENOISE))); 0095 %% 0096 0097 %% Split session into individual runs and analyze each run separately 0098 0099 % In some circumstances, you may wish to derive beta weights from individual 0100 % runs in a session. Here we go through an example where we analyze each 0101 % of the runs in a session separately. Note that denoising cannot be applied 0102 % if there is only a single run to work with (since cross-validation is performed 0103 % across runs). So, instead of calling GLMdenoisedata.m, we will call 0104 % GLMestimatemodel.m. Calling GLMestimatemodel.m also has the benefit that 0105 % we will avoid some of the overhead involved in GLMdenoisedata.m. 0106 0107 % First, we define the ways the runs are to be split. Every run will be 0108 % analyzed separately. 0109 splits = num2cell(1:10); 0110 0111 % For demonstration purposes, we will adopt the strategy of assuming a fixed HRF 0112 % in the analysis of each run. Here we compute a canonical HRF. 0113 hrf = getcanonicalhrf(stimdur,tr)'; 0114 0115 % Now we loop over each run, using GLMestimatemodel.m to fit a GLM to the data. 0116 % Notice that we specify that we want to use an assumed HRF. Also, the "0" 0117 % indicates that we just want to fit the data (no bootstrapping nor cross-validation). 0118 clear resultsIND; 0119 for p=1:length(splits) 0120 ix = splits{p}; 0121 resultsIND(p) = GLMestimatemodel(design(ix),data(ix),stimdur,tr,'assume',hrf,0); 0122 end 0123 %% 0124 0125 % Let's examine the beta weights from each run for an example voxel. 0126 xx = 46; yy = 24; zz = 1; 0127 figure; hold on; 0128 set(gcf,'Units','points','Position',[100 100 700 250]); 0129 betas = []; 0130 for p=1:length(resultsIND) 0131 betas(p,:) = flatten(resultsIND(p).modelmd{2}(xx,yy,zz,:)); 0132 plot(betas(p,:),'b-'); 0133 end 0134 plot(mean(betas,1),'k-','LineWidth',3); 0135 straightline(0,'h','k-'); 0136 xlabel('Condition number'); 0137 ylabel('BOLD signal (% change)'); 0138 title('Beta weights (blue = estimates from individual runs, black = mean across estimates)'); 0139 %%