%PET_ANALYSIS_BOOT Apply Behavioral or Task PLS test and bootstrap test % on PET scan. % % Usage: [brainlv,s,behavlv,brainscores,behavscores,lvcorrs,boot_result, ... % datamatcorrs_lst] = ... % pet_analysis_boot(behavdata_lst,newdata_lst,num_cond_lst,... % num_subj_lst,num_boot,isbehav,Clim) % % See also PLS_BOOT_TEST, PLS_DEVIATION_BOOT_TEST, BEHAVPLS_BOOT, TASKPLS_BOOT % % Called by pet_analysis % % I (behavdata_lst) - A group list of behav data with selected columns; % I (newdata_lst) - A group list of datamat files; % I (num_cond_lst) - A group list of condition numbers; % I (num_subj_lst) - A group list of subject numbers; % I (num_boot) - Number of Bootstrap; % I (isbehav) - 1 if applying Behavioral PLS test; 0 if applying Task PLS test. % I (Clim) - upper limit of confidence interval estimated. % O (brainlv) - Left singular value vector. It is LV for the brain dimension. % O (s) - Singular values vector. % O (behavlv) - Right singular vector. It is LV for the behavior OR contrast % dimension for each scan. % O (brainscores) - The brain score. % O (behavscores) - If use GRAND_MEAN method, behavscores = behavlv, then % expand each condition for all the subjects. If not, it is % calculated for each set of condition, then stack together. % O (lvcorrs) - Correlates brain scores with behavior data for multiple scans, % return [] if for Task PLS or if use GRAND_MEAN method. % O (boot_result) - A Structure array containing the bootstrap result data. % % Created on 03-OCT-2002 by Jimmy Shen for PLS test % Modified on 27-OCT-2002 by Jimmy Shen to add bootstrap test % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function [brainlv,s,behavlv,brainscores,behavscores,lvcorrs, ... boot_result, datamatcorrs_lst] = ... pet_analysis_boot(behavdata_lst,newdata_lst,num_cond_lst,... num_subj_lst,num_boot,isbehav,Clim) % Init % brainlv = []; s = []; behavlv = []; brainscores = []; behavscores = []; lvcorrs = []; boot_result = []; datamatcorrs_lst = {}; stacked_behavdata = []; stacked_datamatcorrs = []; stacked_datamat = []; stacked_mask = []; num_groups = length(newdata_lst); % keeps track of number of times a new bootstrap had to be generated % countnewboot=0; badbeh = []; progress_hdl = rri_progress_ui('initialize'); msg = 'Working on PLS ...'; rri_progress_ui(progress_hdl, '', msg); boot_progress = rri_progress_ui('initialize'); if isbehav [reorder, new_num_boot] = rri_boot_order(num_subj_lst, num_cond_lst(1), num_boot, 0, boot_progress); else [reorder, new_num_boot] = rri_boot_order(num_subj_lst, num_cond_lst(1), num_boot, 1, boot_progress); end if isempty(reorder) return; end; progress_hdl = rri_progress_ui('initialize'); msg = 'Working on PLS ...'; rri_progress_ui(progress_hdl, '', msg); % loop accross the groups, and % calculate datamatcorrs for each group % for i = 1:num_groups k = num_cond_lst(i); n = num_subj_lst(i); datamat = newdata_lst{i}; rri_progress_ui(progress_hdl,'',2/10+5/10*(i-1)/(num_groups)+1/(10*num_groups)); % compute correlation & covariance % if ~isbehav % is Task PLS datamatcorrs = rri_task_mean(datamat,n)-ones(k,1)*mean(datamat); else % is Bahavior PLS datamatcorrs = rri_corr_maps(behavdata_lst{i}, datamat, n, k); end rri_progress_ui(progress_hdl,'',2/10+5/10*(i-1)/(num_groups)+3/(10*num_groups)); % if more than one group, stack data together % if isbehav stacked_behavdata = [stacked_behavdata; behavdata_lst{i}]; end stacked_datamat = [stacked_datamat; datamat]; stacked_datamatcorrs = [stacked_datamatcorrs; datamatcorrs]; datamatcorrs_lst = [datamatcorrs_lst, {datamatcorrs}]; rri_progress_ui(progress_hdl,'',2/10+5/10*(i-1)/(num_groups)+5/(10*num_groups)); end % for % actually, all the groups must have the same condition number % num_cond = num_cond_lst(1); % Singular Value Decomposition % [brainlv,s,behavlv] = svd(stacked_datamatcorrs',0); original_sal = brainlv * s; orig_behavlv = behavlv * s; s = diag(s); rri_progress_ui(progress_hdl,'',9/10); % calculate behav scores % if ~isbehav % is Task PLS brainscores = stacked_datamat * brainlv; % Here, behavlv is actually designlv % according to taskpls.m: fscores=design=testvec(designlv) % so, designscores = designlv num_col = size(behavlv, 2); % expand the num_subj for each row (cond) % did the samething as testvec % for i = 1:num_groups k = num_cond_lst(i); n = num_subj_lst(i); tmp = reshape(behavlv((i-1)*k+1:(i-1)*k+k,:), ... [1, num_col*k]); tmp = repmat(tmp, [n, 1]); % expand to num_subj tmp = reshape(tmp, [n*k, num_col]); behavscores = [behavscores; tmp]; % stack by groups end else % is Behavior PLS [brainscores, behavscores, lvcorrs] = ... rri_get_behavscores(stacked_datamat, stacked_behavdata, ... brainlv, behavlv, num_cond_lst(1), num_subj_lst); orig_corr = lvcorrs; [r1 c1] = size(orig_corr); distrib = zeros(r1, c1, num_boot+1); distrib(1:r1, 1:c1, 1) = orig_corr; end rri_progress_ui(progress_hdl,'',1); % Begin bootstrap loop % % reorder = rri_mkboot_order(num_cond_lst(1), num_subj_lst, num_boot); % move the following code up, so we don't need to run obserbation % test if bootstrap can not run % % [reorder, new_num_boot] = ... % rri_boot_order(num_subj_lst, num_cond_lst(1), num_boot, 1); if new_num_boot ~= num_boot num_boot = new_num_boot; h0 = findobj(0,'tag','PermutationOptionsFigure'); h = findobj(h0,'tag','NumBootstrapEdit'); set(h,'string',num2str(num_boot)); if isbehav distrib = zeros(r1, c1, num_boot+1); distrib(1:r1, 1:c1, 1) = orig_corr; end end max_subj_per_group = 8; if (sum(num_subj_lst <= max_subj_per_group) == num_groups) is_boot_samples = 1; else is_boot_samples = 0; end if isempty(reorder) return; end if ~isbehav sal_sq = zeros(size(brainlv)); sal_sum = zeros(size(brainlv)); dsal_sq = zeros(size(behavlv)); dsal_sum = zeros(size(behavlv)); else sal_sq = original_sal.^2; sal_sum = original_sal; dsal_sq = orig_behavlv.^2; dsal_sum = orig_behavlv; end for p = 1:num_boot msg = ['Working on Bootstrap: ',num2str(p),' out of ',num2str(num_boot)]; rri_progress_ui(progress_hdl, '', msg); rri_progress_ui(progress_hdl,'',p/num_boot); data_p = stacked_datamat(reorder(:,p),:); if isbehav behav_p = stacked_behavdata(reorder(:,p),:); end stacked_data = []; for g=1:num_groups k = num_cond_lst(g); n = num_subj_lst(g); span = sum(num_subj_lst(1:g-1)) * num_cond; if ~isbehav if num_groups == 1 meanmat = rri_task_mean(data_p, n); data = meanmat - (ones(k,1)*mean(meanmat)); else meanmat = rri_task_mean(data_p(1+span:n*k+span,:), n); data = meanmat - (ones(k,1)*mean(meanmat)); end else if ~is_boot_samples % the code below is mainly trying to find a proper % reorder matrix % init badbehav cell array to 0 % which is used to record the bad behav data caused by % bad re-order. This var. is for disp only. % badbehav = zeros(num_cond, size(behavdata_lst{g},2)); % Check for upcoming NaN and re-sample if necessary. % this only happened on behavior analysis, because the % 'xcor' inside of 'corr_maps' contains a 'stdev', which % is a divident. If it is 0, it will cause divided by 0 % problem. % since this happend very rarely, so the speed will not % be affected that much. % % For behavpls_boot, also need to account for multiple % scans and behavs % for c=1:k % traverse all conditions in this group stdmat(c,:) = std(stacked_behavdata(reorder((1+ ... (n*(c-1))+span):(n*c+span), p), :)); end % scanloop % now, check to see if any are zero % while sum(stdmat(:)==0)>0 countnewboot = countnewboot + 1; % keep track of scan & behav that force a resample % badbehav(find(stdmat(:)==0)) = ... badbehav(find(stdmat(:)==0)) + 1; % num_boot is just something to be picked to prevent % infinite loop % if countnewboot > num_boot msgbox(['countnewboot exceeds num_boot: ', ... num2str(countnewboot)],'modal'); breakon=1; break; end % reorder(:,p) = rri_mkboot_order(num_cond_lst(1), ... % num_subj_lst, 1); boot_progress = rri_progress_ui('initialize'); reorder(:,p) = rri_boot_order(num_subj_lst, num_cond_lst(1), 1, 0, boot_progress); for c=1:k % recalc stdmat stdmat(c,:)= std(stacked_behavdata(reorder((1+ ... (n*(c-1))+span):(n*c+span), p), :)); end % scanloop end % while badbeh{g} = badbehav; % save instead of disp end % if ~is_boot_samples % now, we can use this proper reorder matrix % to generate behav_p & data_p, and then % to calculate datamatcoors % behav_p = stacked_behavdata(reorder(:,p),:); data_p = stacked_datamat(reorder(:,p),:); if num_groups == 1 data = rri_corr_maps(behav_p, data_p, n, k); else data = rri_corr_maps(behav_p(1+span:n*k+span,:), ... data_p(1+span:n*k+span,:), n, k); end end stacked_data = [stacked_data; data]; end % for num_groups [pbrainlv, sboot, pbehavlv] = svd(stacked_data',0); % rotate pbehavlv to align with the original behavlv % rotatemat = rri_bootprocrust(behavlv,pbehavlv); % rescale the bootstrap vectors % pbrainlv = pbrainlv * sboot * rotatemat; pbehavlv = pbehavlv * sboot * rotatemat; if isbehav [brainsctmp, behavsctmp, bcorr] = ... rri_get_behavscores(data_p, behav_p, rri_normalize(pbrainlv), ... rri_normalize(pbehavlv), num_cond_lst(1), num_subj_lst); distrib(1:r1, 1:c1, p+1) = bcorr; end sal_sum = sal_sum + pbrainlv; sal_sq = sal_sq + pbrainlv.^2; dsal_sum = dsal_sum + pbehavlv; dsal_sq = dsal_sq + pbehavlv.^2; end % for num_boot boot_result.num_boot = num_boot; boot_result.bootsamp = reorder; if ~isbehav sal_sum2 = (sal_sum.^2) / (num_boot); dsal_sum2 = (dsal_sum.^2) / (num_boot); boot_result.original_sal = original_sal; boot_result.brain_se = sqrt((sal_sq - sal_sum2)/(num_boot-1)); % check for zero standard errors - replace with ones % test_zeros=find(boot_result.brain_se<=0); if ~isempty(test_zeros); boot_result.brain_se(test_zeros)=1; end % now compare the original unstandardized saliances with the % bootstrap salances % boot_result.compare = original_sal ./ boot_result.brain_se; % for zero standard errors - replace bootstrap ratios with zero % since the ratio makes no sense anyway % if ~isempty(test_zeros); boot_result.compare(test_zeros)=0; end boot_result.orig_designlv = orig_behavlv; boot_result.design_se = sqrt((dsal_sq - dsal_sum2)/(num_boot-1)); else sal_sum2 = (sal_sum.^2) / (num_boot+1); dsal_sum2 = (dsal_sum.^2) / (num_boot+1); % compute standard errors - standard deviation of bootstrap sample % since original sample is part of bootstrap, divide by number of % bootstrap iterations rather than number of bootstraps minus 1 % % add ceiling to calculations to prevent the following operations % from producing negative/complex numbers % brain_se = sqrt((ceil(sal_sq)-ceil(sal_sum2))/(num_boot)); behav_se = sqrt((ceil(dsal_sq)-ceil(dsal_sum2))/(num_boot)); % check for zero standard errors - replace with ones % test_zeros=find(brain_se<=0); if ~isempty(test_zeros); brain_se(test_zeros)=1; end compare = original_sal ./ brain_se; compare_behavlv = orig_behavlv ./ behav_se; % for zero standard errors - replace bootstrap ratios with zero % since the ratio makes no sense anyway % if ~isempty(test_zeros); compare(test_zeros)=0; end ul=Clim; ll=100-Clim; % e.g. 0.05 >> 0.025 for upper & lower tails, two-tailed % ClimNi = 0.5*(1-(Clim*0.01)); % loop to calculate upper and lower CI limits % for r=1:r1 for c=1:c1 ulcorr(r,c)=percentile(distrib(r,c,2:num_boot+1),ul); llcorr(r,c)=percentile(distrib(r,c,2:num_boot+1),ll); prop(r,c)=length( find(distrib(r,c,2:num_boot+1) ... <= orig_corr(r,c)) ) / num_boot; if prop(r,c)==1 |prop(r,c)==0 % can't calculate the cumulative_gaussian_inv llcorr_adj(r,c)=NaN; ulcorr_adj(r,c)=NaN; else % adjusted confidence intervals - in case the % bootstrap samples are extremely skewed % norm inverse to start to adjust conf int % ni=cumulative_gaussian_inv(prop(r,c)); % 1st part of recalc the lower conf interval, % this evaluates to +1.96 for 95%CI % uli=(2*ni) + cumulative_gaussian_inv(1-ClimNi); % 1st part of recalc the upper conf interval % e.g -1.96 for 95%CI % lli=(2*ni) + cumulative_gaussian_inv(ClimNi); ncdf_lli=cumulative_gaussian(lli)*100; % percentile for lower bounds ncdf_uli=cumulative_gaussian(uli)*100; % percentile for upper bounds % new percentile % llcorr_adj(r,c)=(percentile(distrib(r,c,2:num_boot+1), ... ncdf_lli)); ulcorr_adj(r,c)=(percentile(distrib(r,c,2:num_boot+1), ... ncdf_uli)); end % if end % for c end % for r boot_result.orig_corr = orig_corr; boot_result.ulcorr = ulcorr; boot_result.llcorr = llcorr; boot_result.ulcorr_adj = ulcorr_adj; boot_result.llcorr_adj = llcorr_adj; boot_result.prop = prop; boot_result.distrib = distrib; boot_result.orig_brainlv = original_sal; boot_result.brain_se = brain_se; boot_result.compare = compare; boot_result.compare_behavlv = compare_behavlv; boot_result.badbeh = badbeh; boot_result.countnewboot = countnewboot; end % enf if ~isbehav return; % pet_analysis_boot %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % load_contrast_file: turn the contrast file into a contrast matrix, % which is used for behavdata or design data. % % I (contrast_file) - filespec of the contrast file % I (num_subj) - number of subjects in each condition % O (contrasts) - contrast matrix % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function contrasts = load_contrast_file(contrast_file, num_subj) load(contrast_file); num_contrasts = length(pls_contrasts); num_conditions = length(pls_contrasts(1).value); contrasts = zeros(num_conditions,num_contrasts); for i=1:num_contrasts contrasts(:,i) = pls_contrasts(i).value'; end % expand each row of condition for all the subjects % tmp = contrasts(:)'; tmp = repmat(tmp,num_subj,1); contrasts = reshape(tmp, num_subj*num_conditions, num_contrasts); return; % load_contrast_file