?? fmri_taskpls_norotate.m
字號:
function [brainlv, s, designlv, b_scores, d_scores, lvintercorrs, design, ...
perm_result, boot_result, dev_evt_list] = ...
fmri_taskpls_norotate(st_datamat, design, num_conditions,...
evt_list, num_boot, num_perm, subj_group, ...
min_subj_per_group,is_boot_samples,boot_samples,new_num_boot)
progress_hdl = ShowProgress('initialize','Working on PLS ... ');
msg = sprintf('Working on PLS ... ');
ShowProgress(progress_hdl,msg)
if ~isempty(progress_hdl)
setappdata(progress_hdl,'ProgressScale',1/(num_perm+1)*0.8);
setappdata(progress_hdl,'ProgressStart',0.2);
ShowProgress(progress_hdl,0)
end;
if isempty(subj_group),
GROUP_ANALYSIS = 0; % for nongroup analysis
else
GROUP_ANALYSIS = 1;
end;
if (GROUP_ANALYSIS == 0),
% dev_data = st_datamat - ones(size(st_datamat,1),1)*mean(st_datamat,1);
% dev_evt_list = evt_list;
[dev_data dev_design] = ...
gen_dev_data(st_datamat,num_conditions,evt_list,design);
dev_evt_list = evt_list;
else
[dev_data dev_design] = ...
gen_grp_dev_data(st_datamat,num_conditions,evt_list,subj_group,design);
% dev_evt_list = repmat([1:num_conditions],1,length(subj_group));
dev_evt_list = evt_list;
end;
% instead of SVD
crossblock = normalize(dev_design)'*dev_data;
brainlv = crossblock';
s = sqrt(sum(crossblock.^2, 2));
designlv = dev_design;
normalized_brainlv = normalize(brainlv);
lvintercorrs = normalized_brainlv'*normalized_brainlv;
if (GROUP_ANALYSIS == 0),
d_scores = designlv(evt_list,:);
else
d_scores = designlv;
% need to expand d_score and reorder with evt_list grp by grp
%
d_scores = expand_d_scores(d_scores, num_conditions, evt_list, subj_group);
end;
b_scores = st_datamat * brainlv;
perm_result = [];
boot_result = [];
if num_perm > 0
% Begin permutation loop
%
ShowProgress(progress_hdl,1);
msg = sprintf('Computing permutation orders ...');
ShowProgress(progress_hdl,msg);
% generate the permutation orders
%
if (GROUP_ANALYSIS == 0),
num_repetitions = length(evt_list) / num_conditions;
perm_order = rri_perm_order(num_repetitions,num_conditions,num_perm);
else
perm_order = rri_perm_order(subj_group,num_conditions,num_perm);
end;
permcount = zeros(size(s));
for k=1:num_perm,
msg = sprintf('Working on permutation ... %d out of %d',k,num_perm);
ShowProgress(progress_hdl,msg);
new_order = perm_order(:,k);
if (GROUP_ANALYSIS == 0), % for nongroup analysis
[data_p design_p] = gen_dev_data(st_datamat(new_order,:),num_conditions, ...
evt_list,design);
else
[data_p design_p] = gen_grp_dev_data(st_datamat(new_order,:),num_conditions, ...
evt_list,subj_group,design);
end;
crossblock = normalize(design_p)'*data_p;
s_perm = sqrt(sum(crossblock.^2, 2));
permcount = permcount + (s_perm >= s);;
ShowProgress(progress_hdl,k+1);
end;
%
%-- permutation loop -------------------------------------------------
perm_result.s_prob = permcount / num_perm;
perm_result.num_perm = num_perm;
perm_result.permsamp = perm_order;
perm_result.sp = permcount;
perm_result.dp = [];
perm_result.designlv_prob = [];
end % end permutation loop
if num_boot > 0
% Begin bootstrap loop
%
max_subj_per_group = 8;
if isempty(subj_group),
num_groups = 0;
GROUP_ANALYSIS = 0; % for nongroup analysis
% moved from below
%
num_rep = length(evt_list) / num_conditions;
% make sure the datamat are blocked by subjects
[sort_evt,sort_idx] = sort(evt_list);
new_row_list = reshape(sort_idx,num_rep,num_conditions)';
new_evt_list = evt_list(new_row_list);
% boot_progress = rri_progress_ui('initialize');
% [boot_order, new_num_boot] = rri_boot_order(num_rep,1,num_boot,0,boot_progress, ...
[boot_order, new_num_boot] = rri_boot_order(num_rep,1,num_boot,0, ...
min_subj_per_group,is_boot_samples,boot_samples,new_num_boot);
if num_rep <= max_subj_per_group
is_boot_samples = 1;
else
is_boot_samples = 0;
end
else
num_groups = length(subj_group);
GROUP_ANALYSIS = 1;
% moved from below
%
% boot_progress = rri_progress_ui('initialize');
% [boot_order, new_num_boot] = rri_boot_order(subj_group,num_conditions,num_boot,0,boot_progress, ...
[boot_order, new_num_boot] = rri_boot_order(subj_group,num_conditions,num_boot,0, ...
min_subj_per_group,is_boot_samples,boot_samples,new_num_boot);
if (sum(subj_group <= length(subj_group)) == num_groups)
is_boot_samples = 1;
else
is_boot_samples = 0;
end
end;
if isempty(boot_order)
boot_result = [];
return;
end
ShowProgress(progress_hdl,1);
msg = sprintf('Compute bootstrap resampling orders ...');
ShowProgress(progress_hdl,msg);
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 ~isempty(progress_hdl)
setappdata(progress_hdl,'ProgressScale',1/(num_boot+1)*0.8);
end
end
if isempty(boot_order)
boot_result = [];
return;
end
brainlv_sq = brainlv.^2;
brainlv_sum = brainlv;
for k=1:num_boot,
msg = sprintf('Working on bootstrap ... %d out of %d',k,num_boot);
ShowProgress(progress_hdl,msg);
new_order = boot_order(:,k);
if (GROUP_ANALYSIS == 0), % for nongroup analysis
new_row_order = new_row_list(:,new_order);
[data_p design_p] = gen_dev_data(st_datamat(new_row_order,:), ...
num_conditions, new_evt_list, design);
else
[data_p design_p]= gen_grp_dev_data(st_datamat(new_order,:), ...
num_conditions, evt_list, subj_group, design);
end;
crossblock = normalize(design_p)'*data_p;
brainlv_sq = brainlv_sq + (crossblock.^2)';
brainlv_sum = brainlv_sum + crossblock';
ShowProgress(progress_hdl,k+1);
end;
%
%-- bootstrap loop -------------------------------------------------
boot_result.num_boot = num_boot;
boot_result.bootsamp = boot_order;
brainlv_sum2 = (brainlv_sum.^2) / (num_boot + 1);
boot_result.brainlv_se=sqrt((brainlv_sq - brainlv_sum2)/num_boot);
% check for zero standard errors - replace with ones
%
test_zeros=find(boot_result.brainlv_se<=0);
boot_result.zero_brain_se = test_zeros;
if ~isempty(test_zeros);
boot_result.brainlv_se(test_zeros)=1;
end
boot_result.compare = brainlv ./ boot_result.brainlv_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
end % end bootstrap loop
return; % fmri_taskls_norotate
%-------------------------------------------------------------------------
function [data, dev_design] = ...
gen_grp_dev_data(st_datamat, num_conditions, evt_list, subj_group, design)
%
% compute the average of st_datamat within the same group
%
num_group = length(subj_group);
g_end_idx = 0;
data = [];
dev_design = [];
for g = 1:num_group,
g_start_idx = g_end_idx + 1;
g_end_idx = g_start_idx+subj_group(g)*num_conditions-1;
g_range = [g_start_idx:g_end_idx];
% store the row indices for each task
%
g_evt_list = zeros(1,length(evt_list));
g_evt_list(g_range) = evt_list(g_range);
task_idx = cell(1,num_conditions);
for i=1:num_conditions,
task_idx{i} = find(g_evt_list==i);
end;
% compute the mean data of each condition for the group
%
if strcmpi(class(st_datamat),'single')
mean_datamat = single(zeros(num_conditions,size(st_datamat,2)));
else
mean_datamat = zeros(num_conditions,size(st_datamat,2));
end
for i=1:num_conditions,
mean_datamat(i,:) = mean(st_datamat(task_idx{i},:),1);
end;
grp_datamat = mean_datamat; % - ones(num_conditions,1)*mean(mean_datamat,1);
data = [data; grp_datamat];
span = [((g-1)*num_conditions+1) : g*num_conditions];
dev_design = [dev_design; design(span,:)];
end;
return; % gen_grp_dev_data
%-------------------------------------------------------------------------
function [data, design] = ...
gen_dev_data(st_datamat, num_conditions, evt_list, design)
task_idx = cell(1,num_conditions);
for i=1:num_conditions,
task_idx{i} = find(evt_list==i);
end;
% compute the mean data of each condition for the group
%
if strcmpi(class(st_datamat),'single')
mean_datamat = single(zeros(num_conditions,size(st_datamat,2)));
else
mean_datamat = zeros(num_conditions,size(st_datamat,2));
end
for i=1:num_conditions,
mean_datamat(i,:) = mean(st_datamat(task_idx{i},:),1);
end;
data = mean_datamat; % - ones(num_conditions,1)*mean(st_datamat,1);
return; % gen_dev_data
%-------------------------------------------------------------------------
function hdl = ShowProgress(progress_hdl,info)
% 'initialize' - return progress handle if any
%
if ischar(progress_hdl) & strcmp(lower(progress_hdl),'initialize'),
if ~isempty(gcf) & isequal(get(gcf,'Tag'),'ProgressFigure'),
hdl = gcf;
if ~isempty(info)
set(hdl,'Name',info);
end;
else
hdl = [];
end;
return;
end;
if ~isempty(progress_hdl)
if ischar(info)
rri_progress_status(progress_hdl,'Show_message',info);
else
rri_progress_status(progress_hdl,'Update_bar',info);
end;
return;
end;
if ischar(info),
disp(info)
end;
return; % ShowProgress
%-------------------------------------------------------------------------
function new_d_scores = expand_d_scores(d_scores,num_conditions,evt_list,subj_group)
new_d_scores = [];
num_in_grp = [0 num_conditions*subj_group];
for grp_idx = 1:length(subj_group)
num_subjs = subj_group(grp_idx);
first = sum(num_in_grp(1:grp_idx)) + 1;
last = sum(num_in_grp(1:(grp_idx+1)));
tmp_evt_list = evt_list(first:last);
tmp_d_scores = d_scores(((grp_idx-1)*num_conditions+1):grp_idx*num_conditions,:);
tmp_d_scores = tmp_d_scores(tmp_evt_list,:); % expand and reorder
new_d_scores = [new_d_scores;tmp_d_scores];
end
return;
?? 快捷鍵說明
復制代碼
Ctrl + C
搜索代碼
Ctrl + F
全屏模式
F11
切換主題
Ctrl + Shift + D
顯示快捷鍵
?
增大字號
Ctrl + =
減小字號
Ctrl + -