% Example 1: Run GLMdenoise on an example dataset
0001 %% Example 1: Run GLMdenoise on an example dataset 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 %% Inspect the data 0017 0018 % Look at the design matrix for the first run. Notice that there are 0019 % 35 conditions and each condition is presented once during the run. 0020 figure; 0021 set(gcf,'Units','points','Position',[100 100 350 500]); 0022 imagesc(design{1}); 0023 colormap(gray); 0024 colorbar; 0025 xlabel('Conditions'); 0026 ylabel('Time points'); 0027 title(sprintf('Design matrix (%d time points, %d conditions)',size(design{1},1),size(design{1},2))); 0028 %% 0029 0030 % Look at the first volume of the first run. We split the volume along 0031 % the third dimension (the slice dimension) and show slices starting at the 0032 % upper-left, proceeding downwards, and then from left to right. 0033 figure; 0034 imagesc(makeimagestack(data{1}(:,:,:,1))); 0035 colormap(gray); 0036 axis equal tight off; 0037 colorbar; 0038 title('fMRI data (first volume)'); 0039 %% 0040 0041 % Check various constants 0042 fprintf('There are %d runs in total.\n',length(design)); 0043 fprintf('The dimensions of the data for the first run are %s.\n',mat2str(size(data{1}))); 0044 fprintf('The stimulus duration is %.6f seconds.\n',stimdur); 0045 fprintf('The sampling rate (TR) is %.6f seconds.\n',tr); 0046 %% 0047 0048 %% Run GLMdenoise 0049 0050 % Call GLMdenoise using default parameters. Outputs will be stored in 0051 % 'results' and 'denoiseddata'. Figures will be written to 'example1figures' 0052 % in the current directory. 0053 [results,denoiseddata] = GLMdenoisedata(design,data,stimdur,tr,[],[],[],'example1figures'); 0054 %% 0055 0056 % For comparison purposes, we call GLMdenoise again but turn off the global 0057 % noise regressors. This is done by simply setting the maximum number of PCs 0058 % to 0. Outputs from this call will be stored in 'resultsALT' and 'denoiseddataALT', 0059 % and a separate set of figures will be written to 'example1figuresALT'. 0060 % We will inspect the results of this GLMdenoise call at a later point. 0061 [resultsALT,denoiseddataALT] = GLMdenoisedata(design,data,stimdur,tr,[],[],struct('numpcstotry',0),'example1figuresALT'); 0062 %% 0063 0064 %% Inspect figures 0065 0066 % GLMdenoisedata writes out a number of figures illustrating the various 0067 % computations that are performed. It is useful to browse through these 0068 % figures to check sanity. Here we highlight the most useful figures, 0069 % but there are additional figures (all figures are described in the 0070 % documentation in GLMdenoisedata.m.). 0071 0072 % 'MeanVolume.png' shows the mean across all volumes. 0073 % This is useful as a frame of reference. 0074 figure; 0075 imagesc(imread('example1figures/MeanVolume.png'),[0 255]); 0076 colormap(gray); 0077 axis equal tight off; 0078 title('Mean volume'); 0079 %% 0080 0081 % 'FinalModel.png' shows the R^2 of the final model fit. The color range is 0082 % 0% to 100% (but is nonlinear; see the color bar). These R^2 values are not 0083 % cross-validated (thus, even voxels with no signal have positive R^2 values). 0084 figure; 0085 imagesc(imread('example1figures/FinalModel.png'),[0 255]); 0086 colormap(hot); 0087 axis equal tight off; 0088 cb = colorbar; 0089 set(cb,'YTick',linspace(0,255,11),'YTickLabel',round((0:.1:1).^2 * 100)); 0090 title('Final model R^2 (not cross-validated)'); 0091 %% 0092 0093 % 'NoisePool.png' shows which voxels were used for the noise pool. Assuming 0094 % the default parameters, the noise pool consists of voxels whose (1) mean 0095 % intensity is greater than one-half of the 99th percentile of mean intensity 0096 % values and (2) cross-validated R^2 value is less than 0%. Notice that 0097 % the noise pool omits voxels near the top of the slice because the raw 0098 % signal level is relatively low there (due to the distance from the RF coil). 0099 figure; 0100 imagesc(imread('example1figures/NoisePool.png'),[0 255]); 0101 colormap(gray); 0102 axis equal tight off; 0103 title('Noise pool'); 0104 %% 0105 0106 % 'HRF.png' shows the initial HRF guess (red) and the final HRF estimate (blue). 0107 % Full flexibility is given to the HRF estimate, so the noisiness of the HRF 0108 % estimate gives an indication of the SNR level in the data. In this case, 0109 % the HRF estimate looks very nice, suggesting that the data have good SNR. 0110 % The first point in the HRF coincides with the onset of a condition (which is 0111 % denoted by 1s in the design matrix). The granularity of the HRF reflects 0112 % the sampling rate (TR) of the data. 0113 figure; 0114 imageactual('example1figures/HRF.png'); 0115 %% 0116 0117 % 'PCcrossvalidationXX.png' shows cross-validated R^2 values corresponding 0118 % to different numbers of PCs. The color range is 0% to 100%. Here we show 0119 % the initial model (no PCs) and the final model (number of PCs selected by 0120 % the code). Notice that because of cross-validation, most voxels are black, 0121 % indicating cross-validated R^2 values that are 0% or less. Also, notice 0122 % that there is some gain in cross-validation performance in the middle of the 0123 % slices. The most effective way to view and interpret these figures 0124 % is to flip between them in quick succession, and this is easiest to do 0125 % in your operating system (not in MATLAB). 0126 figure; 0127 set(gcf,'Units','points','Position',[100 100 900 325]); 0128 subplot(1,2,1); 0129 imagesc(imread('example1figures/PCcrossvalidation00.png'),[0 255]); 0130 colormap(hot); 0131 axis equal tight off; 0132 cb = colorbar; 0133 set(cb,'YTick',linspace(0,255,11),'YTickLabel',round((0:.1:1).^2 * 100)); 0134 title('Initial model (PC = 0) cross-validated R^2'); 0135 subplot(1,2,2); 0136 imagesc(imread(sprintf('example1figures/PCcrossvalidation%02d.png',results.pcnum))); 0137 colormap(hot); 0138 axis equal tight off; 0139 cb = colorbar; 0140 set(cb,'YTick',linspace(0,255,11),'YTickLabel',round((0:.1:1).^2 * 100)); 0141 title(sprintf('Final model (PC = %d) cross-validated R^2',results.pcnum)); 0142 %% 0143 0144 % 'PCscatterXX.png' shows scatter plots of the cross-validation performance of the 0145 % initial model (no PCs) against the performance of individual models with 0146 % different numbers of PCs. Here we show the scatter plot corresponding to 0147 % the final model (number of PCs selected by the code). Voxels that were 0148 % used to determine the optimal number of PCs are shown in red; other voxels 0149 % are shown in green. This figure shows that adding PCs produces a clear and 0150 % consistent improvement in cross-validation performance. 0151 figure; 0152 imageactual(sprintf('example1figures/PCscatter%02d.png',results.pcnum)); 0153 %% 0154 0155 % 'PCselection.png' illustrates how the number of PCs was selected. The median 0156 % cross-validated R^2 value across a subset of the voxels (marked as red 0157 % in the scatter plots) is plotted as a line, and the selected number 0158 % of PCs is marked with a circle. Assuming the default parameters, 0159 % the code chooses the minimum number of PCs such that the improvement 0160 % relative to the initial model is within 5% of the maximum improvement. 0161 % This is a slightly conservative selection strategy, designed to avoid 0162 % overfitting and to be robust across different datasets. 0163 figure; 0164 imageactual('example1figures/PCselection.png'); 0165 %% 0166 0167 %% Inspect outputs 0168 0169 % The outputs of GLMdenoisedata are contained in 0170 % the variables 'results' and 'denoiseddata'. 0171 % Here we do some basic inspections of the outputs. 0172 0173 % Select a voxel to inspect. This is done by finding voxels that have 0174 % cross-validated R^2 values between 0% and 5% under the initial model (no PCs), 0175 % and then selecting the voxel that shows the largest improvement when 0176 % using the final model. 0177 ix = find(results.pcR2(:,:,:,1) > 0 & results.pcR2(:,:,:,1) < 5); 0178 improvement = results.pcR2(:,:,:,1+results.pcnum) - results.pcR2(:,:,:,1); 0179 [mm,ii] = max(improvement(ix)); 0180 ix = ix(ii); 0181 [xx,yy,zz] = ind2sub(results.inputs.datasize{1}(1:3),ix); 0182 0183 % Inspect amplitude estimates for the example voxel. The first row shows results 0184 % obtained from the GLMdenoisedata call that does not involve global noise regressors; 0185 % the second row shows results obtained from the regular GLMdenoisedata call. 0186 % The left panels show amplitude estimates from individual bootstraps, whereas 0187 % the right panels show the median and 68% interval across bootstraps (this 0188 % corresponds to the final model estimate and the estimated error on the 0189 % estimate, respectively). 0190 figure; 0191 set(gcf,'Units','points','Position',[100 100 700 500]); 0192 for p=1:2 0193 if p==1 0194 ampboots = squeeze(resultsALT.models{2}(xx,yy,zz,:,:)); % conditions x boots 0195 amp = flatten(resultsALT.modelmd{2}(xx,yy,zz,:)); % 1 x conditions 0196 ampse = flatten(resultsALT.modelse{2}(xx,yy,zz,:)); % 1 x conditions 0197 else 0198 ampboots = squeeze(results.models{2}(xx,yy,zz,:,:)); % conditions x boots 0199 amp = flatten(results.modelmd{2}(xx,yy,zz,:)); % 1 x conditions 0200 ampse = flatten(results.modelse{2}(xx,yy,zz,:)); % 1 x conditions 0201 end 0202 n = length(amp); 0203 subplot(2,2,(p-1)*2 + 1); hold on; 0204 plot(ampboots); 0205 straightline(0,'h','k-'); 0206 xlabel('Condition number'); 0207 ylabel('BOLD signal (% change)'); 0208 title('Amplitude estimates (individual bootstraps)'); 0209 subplot(2,2,(p-1)*2 + 2); hold on; 0210 bar(1:length(amp),amp,1); 0211 errorbar2(1:length(amp),amp,ampse,'v','r-'); 0212 if p==1 0213 ax = axis; axis([0 n+1 ax(3:4)]); ax = axis; 0214 end 0215 xlabel('Condition number'); 0216 ylabel('BOLD signal (% change)'); 0217 title('Amplitude estimates (median and 68% interval)'); 0218 end 0219 for p=1:2 0220 subplot(2,2,(p-1)*2 + 1); axis(ax); 0221 subplot(2,2,(p-1)*2 + 2); axis(ax); 0222 end 0223 %% 0224 0225 % Compare SNR before and after the use of global noise regressors. 0226 % To focus on voxels that related to the experiment, we select voxels with 0227 % cross-validated R^2 values that are greater than 0% under the initial model. 0228 % To ensure that the SNR values reflect only changes in the noise level, 0229 % we ignore the SNR computed in each individual GLMdenoisedata call and 0230 % re-compute SNR, holding the numerator (the signal) constant across the 0231 % two calls. (Note: GLMdenoisedata automatically writes out a figure that 0232 % shows a before-and-after SNR comparison (SNRcomparebeforeandafter.png). 0233 % What is shown in this script does the comparison manually just for sake 0234 % of example.) 0235 ok = results.pcR2(:,:,:,1) > 0; 0236 signal = mean([results.signal(ok) resultsALT.signal(ok)],2); 0237 snr1 = signal ./ resultsALT.noise(ok); 0238 snr2 = signal ./ results.noise(ok); 0239 figure; hold on; 0240 scatter(snr1,snr2,'r.'); 0241 ax = axis; 0242 mx = max(ax(3:4)); 0243 axissquarify; 0244 axis([0 mx 0 mx]); 0245 xlabel('SNR before'); 0246 ylabel('SNR after'); 0247 %% 0248 0249 % An alternative to using the GLM estimates provided by GLMdenoisedata is 0250 % to use 'denoiseddata', which contains the original time-series data but 0251 % with the component of the data that is estimated to be due to the global 0252 % noise regressors subtracted off. Here we inspect the denoised data for 0253 % the same example voxel examined earlier. Note that the data components 0254 % that are present in 'denoiseddata' can be customized (see opt.denoisespec 0255 % in GLMdenoisedata.m). The default parameters leave in the estimated baseline 0256 % signal drift, which explains the drift in the plotted time-series. 0257 figure; hold on; 0258 set(gcf,'Units','points','Position',[100 100 700 250]); 0259 data1 = flatten(data{1}(xx,yy,zz,:)); 0260 data2 = flatten(denoiseddata{1}(xx,yy,zz,:)); 0261 n = length(data1); 0262 h1 = plot(data1,'r-'); 0263 h2 = plot(data2,'b-'); 0264 ax = axis; axis([0 n+1 ax(3:4)]); 0265 legend([h1 h2],{'Original' 'Denoised'}); 0266 xlabel('Time point'); 0267 ylabel('MR signal'); 0268 %%