diff --git a/dbs_cal_cp.m b/dbs_cal_cp.m new file mode 100644 index 0000000..88ef010 --- /dev/null +++ b/dbs_cal_cp.m @@ -0,0 +1,60 @@ +function cp = dbs_cal_cp (numNode, numperm, thrClst, initthre_range, dbs, perm_cp, pwr) +% DBS_CAL_CP +% ================================================================================================================ +% [ INPUTS ] +% numperm = the number of permutations performed for family (default = 5000). +% thrClst = DBS-based FWE-corrected cluster-level threshold for significance (default = 0.05) +% numNode, initthre_range, dbs, perm_cp, pwr +% ---------------------------------------------------------------------------------------------------------------- +% [ OUTPUTS ] +% cp +% ---------------------------------------------------------------------------------------------------------------- +% Last update: Aug 30, 2016. +% +% Copyright 2016. Kwangsun Yoo (K Yoo), PhD +% E-mail: rayksyoo@gmail.com / raybeam@kaist.ac.kr +% Laboratory for Cognitive Neuroscience and NeuroImaging (CNI) +% Department of Bio and Brain Engineering +% Korea Advanced Instititue of Science and Technology (KAIST) +% Daejeon, Republic of Korea +% ================================================================================================================ + +for thre = 1 : length(initthre_range) + cp.bd_max(thre,1) = dbs.bd.thre{thre}; +end + +for thre_limit = 2 : length(initthre_range) + cpMat = dbs_estm_cp(dbs.wd.org, initthre_range, thre_limit, numNode, pwr); + cp.org{thre_limit,1} = cpMat; + cp.org_sort{thre_limit,1} = sort(cp.org{thre_limit}, 'descend'); + perm_cp.max_sort{thre_limit,1} = sort(perm_cp.max{thre_limit}, 'descend'); +end + +cp.max_sort = perm_cp.max_sort; % clear perm_cp + +thrClst_int = ceil(thrClst * numperm); +cp.org_sig = cp.org; + +for thre_limit = 2 : length(initthre_range) + cp.thre(thre_limit,1) = perm_cp.max_sort{thre_limit}(thrClst_int); + cp.norm{thre_limit,1} = cp.org{thre_limit,1} / cp.thre(thre_limit); + cp.norm_sort{thre_limit,1} = sort(cp.norm{thre_limit,1}, 'descend'); + + cp.org_sig{thre_limit,1} = cp.org{thre_limit}; + cp.org_sig{thre_limit}(cp.org_sig{thre_limit} <= cp.thre(thre_limit)) = 0; + + cp.norm_sig{thre_limit,1} = cp.norm{thre_limit}; + cp.norm_sig{thre_limit}(cp.norm_sig{thre_limit} <= 1) = 0; + + cp.clst{thre_limit,1}(:,1) = find(cp.org_sig{thre_limit} > 0); + cp.clst{thre_limit,1}(:,2) = cp.org{thre_limit}(cp.org_sig{thre_limit} > 0); + + for i_cp = 1: size(cp.clst{thre_limit}, 1) + i_find = length(initthre_range); + while i_find > 0 + if find(dbs.wd.clst{i_find}(:,1) == cp.clst{thre_limit}(i_cp,1)) + break; end; + i_find = i_find-1; end; + cp.clstnode{thre_limit,1}{i_cp,1}(:,1) = find(dbs.s{i_find}(cp.clst{thre_limit}(i_cp,1), :) > 0); + end; +end; diff --git a/dbs_cal_df.m b/dbs_cal_df.m new file mode 100644 index 0000000..af1b958 --- /dev/null +++ b/dbs_cal_df.m @@ -0,0 +1,33 @@ +function dof = dbs_cal_df (hypotest , numSamples) +% DBS_CAL_DF Provide the degree of freedom (dof) for the hypothesis testing and the number of samples +% ================================================================================================================ +% [ INPUTS ] +% hypotest = type of test (default = 0). +% 0: two-sample paired t-test (ttest) +% 1: two-sample unpaired t-test (ttest2) (assumping the same variance for the two groups) +% 2: correlation analysis (corr) +% +% numSaples = the number of samples (or subjects) give to the hypothesis testing +% ---------------------------------------------------------------------------------------------------------------- +% [ OUTPUTS ] +% dof = the degree of freedom (DOF) for the hypothesis testing and the number of samples +% ---------------------------------------------------------------------------------------------------------------- +% Last update: Aug 30, 2016. +% +% Copyright 2016. Kwangsun Yoo (K Yoo), PhD +% E-mail: rayksyoo@gmail.com / raybeam@kaist.ac.kr +% Laboratory for Cognitive Neuroscience and NeuroImaging (CNI) +% Department of Bio and Brain Engineering +% Korea Advanced Instititue of Science and Technology (KAIST) +% Daejeon, Republic of Korea +% ================================================================================================================ + +if hypotest == 0 + dof = numSamples/2 - 1 ; + +elseif hypotest == 1 + dof = numSamples - 2 ; % for equal variance + +elseif hypotest == 2 + dof = numSamples - 2 ; +end \ No newline at end of file diff --git a/dbs_correction.m b/dbs_correction.m new file mode 100644 index 0000000..ecc7f79 --- /dev/null +++ b/dbs_correction.m @@ -0,0 +1,44 @@ +function dbs = dbs_correction(dbs_thre, icft, s, p) +% DBS_CORRECTION +% ================================================================================================================ +% [ INPUTS ] +% dbs_thre, icft, s, p +% ---------------------------------------------------------------------------------------------------------------- +% [ OUTPUTS ] +% dbs +% ---------------------------------------------------------------------------------------------------------------- +% Last update: Aug 30, 2016. +% +% Copyright 2016. Kwangsun Yoo (K Yoo), PhD +% E-mail: rayksyoo@gmail.com / raybeam@kaist.ac.kr +% Laboratory for Cognitive Neuroscience and NeuroImaging (CNI) +% Department of Bio and Brain Engineering +% Korea Advanced Instititue of Science and Technology (KAIST) +% Daejeon, Republic of Korea +% ================================================================================================================ + +s2 = abs(s); +for thre = 1:length(icft.s.range) + clear a b s_thre + [a, s_thre] = dbs_estm_dgr(p, s, icft.p.range(thre), icft.s.range(thre)); + dbs.bd.org{thre,1} = a.bd.org; + dbs.wd.org{thre,1} = a.wd.org; + dbs.s{thre,1} = s_thre; + + b = dbs_estm_sig(dbs.bd.org{thre}, dbs_thre.bd.thre{thre}, dbs.wd.org{thre}, dbs_thre.wd.thre{thre}); + dbs.bd.clst{thre,1} = b.bd.clst; + dbs.wd.clst{thre,1} = b.wd.clst; + for i_clst = 1 : size(dbs.bd.clst{thre},1) + dbs.bd.clstnode{thre,1}{i_clst,1}(:,1) = find(dbs.s{thre}(dbs.bd.clst{thre}(i_clst,1), :) ); + dbs.bd.clstnode{thre,1}{i_clst,1}(:,2) = s2(dbs.bd.clst{thre}(i_clst,1), dbs.bd.clstnode{thre,1}{i_clst,1}(:,1)); + dbs.bd.clstnode{thre,1}{i_clst,1}(:,3) = dbs.s{thre}(dbs.bd.clst{thre}(i_clst,1), dbs.bd.clstnode{thre,1}{i_clst,1}(:,1)); + end + for i_clst = 1 : size(dbs.wd.clst{thre},1) + dbs.wd.clstnode{thre,1}{i_clst,1}(:,1) = find(dbs.s{thre}(dbs.wd.clst{thre}(i_clst,1), :) ) ; + dbs.wd.clstnode{thre,1}{i_clst,1}(:,2) = s2(dbs.wd.clst{thre}(i_clst,1), dbs.wd.clstnode{thre,1}{i_clst,1}(:,1)); + dbs.wd.clstnode{thre,1}{i_clst,1}(:,3) = dbs.s{thre}(dbs.wd.clst{thre}(i_clst,1), dbs.wd.clstnode{thre,1}{i_clst,1}(:,1)); + end +end + +dbs.org.s{1,1} = s; dbs.org.s{2,1} = abs(dbs.org.s{1}); dbs.org.p = p; +dbs.bd.thre = dbs_thre.bd.thre; dbs.wd.thre = dbs_thre.wd.thre; diff --git a/dbs_estm_cp.m b/dbs_estm_cp.m new file mode 100644 index 0000000..2c3b840 --- /dev/null +++ b/dbs_estm_cp.m @@ -0,0 +1,28 @@ +function cpMat = dbs_estm_cp(wd_org, initthre_range, thre_limit, numnode, pwr) +% DBS_ESTM_CP +% ================================================================================================================ +% [ INPUTS ] +% wd_org, initthre_range, thre_limit, numnode, pwer% +% ---------------------------------------------------------------------------------------------------------------- +% [ OUTPUTS ] +% cpMat +% ---------------------------------------------------------------------------------------------------------------- +% Last update: Aug 30, 2016. +% +% Copyright 2016. Kwangsun Yoo (K Yoo), PhD +% E-mail: rayksyoo@gmail.com / raybeam@kaist.ac.kr +% Laboratory for Cognitive Neuroscience and NeuroImaging (CNI) +% Department of Bio and Brain Engineering +% Korea Advanced Instititue of Science and Technology (KAIST) +% Daejeon, Republic of Korea +% ================================================================================================================ + +cpMat = zeros(numnode,1); + +for thre = 1 : thre_limit + if thre == thre_limit + cpMat = cpMat + (wd_org{thre} .^ pwr.dgr * initthre_range(thre) ^ pwr.thr * (initthre_range(thre) - initthre_range(thre-1)) ) ; + else + cpMat = cpMat + (wd_org{thre} .^ pwr.dgr * initthre_range(thre) ^ pwr.thr * (initthre_range(thre+1) - initthre_range(thre)) ) ; + end +end \ No newline at end of file diff --git a/dbs_estm_dgr.m b/dbs_estm_dgr.m new file mode 100644 index 0000000..a420765 --- /dev/null +++ b/dbs_estm_dgr.m @@ -0,0 +1,36 @@ +function [dbs_estm, s3] = dbs_estm_dgr (p, s, thr_p, thr_s) +% DBS_ESTM_DGR Estimation of the degree in connectivitiy matrix for a initial cluster-forming threshold. +% ================================================================================================================ +% [ INPUTS ] +% p, s, thre_p, thr_s +% ---------------------------------------------------------------------------------------------------------------- +% [ OUTPUTS ] +% dbs_estm, s3 +% ---------------------------------------------------------------------------------------------------------------- +% Last update: Aug 30, 2016. +% +% Copyright 2016. Kwangsun Yoo (K Yoo), PhD +% E-mail: rayksyoo@gmail.com / raybeam@kaist.ac.kr +% Laboratory for Cognitive Neuroscience and NeuroImaging (CNI) +% Department of Bio and Brain Engineering +% Korea Advanced Instititue of Science and Technology (KAIST) +% Daejeon, Republic of Korea +% ================================================================================================================ + +p2 = p; +p2(isnan(p2)) = 1; +p2(p2 > thr_p) = 1; +dbs_estm.h = 1- p2; +dbs_estm.h(dbs_estm.h > 0) = 1; + +dbs_estm.bd.org = sum(dbs_estm.h)'; +dbs_estm.bd.max = max(dbs_estm.bd.org); + +s2 = s; +s2(isnan(s2)) = 0; +s2(dbs_estm.h == 0) = 0; +s3 = abs(s2); +s3(s3 > 0) = s3(s3 > 0) - abs(thr_s); + +dbs_estm.wd.org = sum(s3)'; +dbs_estm.wd.max = max(dbs_estm.wd.org); diff --git a/dbs_estm_dgr_range.m b/dbs_estm_dgr_range.m new file mode 100644 index 0000000..6c1d1b9 --- /dev/null +++ b/dbs_estm_dgr_range.m @@ -0,0 +1,27 @@ +function dgr_range = dbs_estm_dgr_range (p, s, icft) +% DBS_ESTM_DGR_RANGE +% ================================================================================================================ +% [ INPUTS ] +% p, s, icft +% ---------------------------------------------------------------------------------------------------------------- +% [ OUTPUTS ] +% dgr_range +% ---------------------------------------------------------------------------------------------------------------- +% Last update: Aug 30, 2016. +% +% Copyright 2016. Kwangsun Yoo (K Yoo), PhD +% E-mail: rayksyoo@gmail.com / raybeam@kaist.ac.kr +% Laboratory for Cognitive Neuroscience and NeuroImaging (CNI) +% Department of Bio and Brain Engineering +% Korea Advanced Instititue of Science and Technology (KAIST) +% Daejeon, Republic of Korea +% ================================================================================================================ + +for thre = 1:length(icft.s.range) + clear a + a = dbs_estm_dgr(p, s, icft.p.range(thre), icft.s.range(thre)); + dgr_range.bd.max{thre}(i_perm,1) = a.bd.max; + dgr_range.bd.org{thre} = a.bd.org; + dgr_range.wd.max{thre}(i_perm,1) = a.wd.max; + dgr_range.wd.org{thre} = a.wd.org; +end \ No newline at end of file diff --git a/dbs_estm_sig.m b/dbs_estm_sig.m new file mode 100644 index 0000000..82c4f6e --- /dev/null +++ b/dbs_estm_sig.m @@ -0,0 +1,28 @@ +function dbs = dbs_estm_sig(bdgr, bdgr_thre, wdgr, wdgr_thre) +% DBS_ESTM_SIG Estimation of the FWE-corrected significance of degree and cluster persistency of every node +% ================================================================================================================ +% [ INPUTS ] +% bdgr, bdgr_thre, wdgr, wdgr_thr +% ---------------------------------------------------------------------------------------------------------------- +% [ OUTPUTS ] +% dbs +% ---------------------------------------------------------------------------------------------------------------- +% Last update: Aug 30, 2016. +% +% Copyright 2016. Kwangsun Yoo (K Yoo), PhD +% E-mail: rayksyoo@gmail.com / raybeam@kaist.ac.kr +% Laboratory for Cognitive Neuroscience and NeuroImaging (CNI) +% Department of Bio and Brain Engineering +% Korea Advanced Instititue of Science and Technology (KAIST) +% Daejeon, Republic of Korea +% ================================================================================================================ + +bdgr_sig = bdgr; +bdgr_sig (bdgr < bdgr_thre) = 0; +dbs.bd.clst(:,1) = find(bdgr_sig > 0); +dbs.bd.clst(:,2) = bdgr(bdgr_sig > 0); + +wdgr_sig = wdgr; +wdgr_sig (wdgr <= wdgr_thre) = 0; +dbs.wd.clst(:,1) = find(wdgr_sig > 0); +dbs.wd.clst(:,2) = wdgr(wdgr_sig > 0); diff --git a/dbs_hypo_test.m b/dbs_hypo_test.m new file mode 100644 index 0000000..fcce20b --- /dev/null +++ b/dbs_hypo_test.m @@ -0,0 +1,43 @@ +function [s, t, p] = dbs_hypo_test(setmat, label, hypotest, direction) +% DBS_HYPO_TEST Set the hypothesis testing +% ================================================================================================================ +% [ INPUTS ] +% setmat = 3-D matrix which consists of a set of 2-D matrices from multiple subjects. +% A size of setmat should be [N by N by M]. +% N: the number of nodes. +% M: the number of subjects +% +% label = 1-D vector containing a list of labels +% with a value 0 (group 1) or 1 (group 2), indicating in which group each subject is included (for hypoTest = 0 or 1) +% or with individual measures of behavioral performance from each subjects for correlation (for hypoTest = 2) +% +% hypoTest = type of test (default = 0). +% 0: two-sample paired t-test (ttest) +% 1: two-sample unpaired t-test (ttest2) (assumping the same variance for the two groups) +% 2: correlation analysis (corr) +% +% direction +% 0: g1 = g2 (two-tail) +% 1: g1 > g2 (one-tail) +% -1: g2 < g1 (one-tail) +% ---------------------------------------------------------------------------------------------------------------- +% [ OUTPUTS ] +% s, t, p +% ---------------------------------------------------------------------------------------------------------------- +% Last update: Aug 30, 2016. +% +% Copyright 2016. Kwangsun Yoo (K Yoo), PhD +% E-mail: rayksyoo@gmail.com / raybeam@kaist.ac.kr +% Laboratory for Cognitive Neuroscience and NeuroImaging (CNI) +% Department of Bio and Brain Engineering +% Korea Advanced Instititue of Science and Technology (KAIST) +% Daejeon, Republic of Korea +% ================================================================================================================ + +if hypotest == 0 + [t, p] = dbs_pairT_matrix (setmat(:,:,label==0) , setmat(:,:,label==1), direction); s = t; +elseif hypotest == 1 + [t, p] = dbs_unpairT_matrix (setmat(:,:,label==0) , setmat(:,:,label==1), direction); s = t; +elseif hypotest == 2 + [s, t, p] = dbs_corr_matrix (setmat, label, direction); +end diff --git a/dbs_main.m b/dbs_main.m new file mode 100644 index 0000000..e58986a --- /dev/null +++ b/dbs_main.m @@ -0,0 +1,102 @@ +function dbs = dbs_main(setmat, label, hypoTest, direction, numperm, thrClst) +% DBS_MAIN FWE correction using Degree-based statistics (DBS) in connectivity analysis +% ================================================================================================================ +% [ INPUTS ] +% setmat = 3-D matrix which consists of a set of 2-D matrices from multiple subjects. +% A size of setmat should be [N by N by M]. +% N: the number of nodes. +% M: the number of subjects +% +% label = 1-D vector containing a list of labels +% with a value 0 (group 1) or 1 (group 2), indicating in which group each subject is included (for hypoTest = 0 or 1) +% or with individual measures of behavioral performance from each subjects for correlation (for hypoTest = 2) +% +% hypoTest = type of test (default = 0). +% 0: two-sample paired t-test (ttest) +% 1: two-sample unpaired t-test (ttest2) (assumping the same variance for the two groups) +% 2: correlation analysis (corr) +% +% direction +% 0: g1 = g2 (two-tail) +% 1: g1 > g2 (one-tail) +% -1: g2 < g1 (one-tail) +% +% numperm = the number of permutations performed for family (default = 5000). +% +% thrClst = DBS-based FWE-corrected cluster-level threshold for significance (default = 0.05) +% ---------------------------------------------------------------------------------------------------------------- +% [ OUTPUTS ] +% dbs.bd = DBS corrected results from the binary degree +% dbs.wd = DBS corrected results from the wiehgted degree +% dbs.cp = DBS corrected results using 'Cluster Persistency (cp) score' with the weighted degree +% +% dbs.s = +% dbs.org = +% dbs.time = time taken to run DBS analysis +% dbs.icft = a range of initial cluster-forming thresholds in analysis +% ---------------------------------------------------------------------------------------------------------------- +% Last update: Aug 30, 2016. +% +% Copyright 2016. Kwangsun Yoo (K Yoo), PhD +% E-mail: rayksyoo@gmail.com / raybeam@kaist.ac.kr +% Laboratory for Cognitive Neuroscience and NeuroImaging (CNI) +% Department of Bio and Brain Engineering +% Korea Advanced Instititue of Science and Technology (KAIST) +% Daejeon, Republic of Korea +% +% Paper: Degree-based statistic and center persistency for brain connectivity analysis (2016) Human Brain Mapping. +% ================================================================================================================ + +temp = clock; fprintf('\n\t\t%d, %d/%d, %d:%2.0f:%2.0f\tProcess started\n', temp) + +%% Check the input argument +if nargin < 3; error('Must input at least three parameters.\n'); end + +if ~exist('direction'); direction = 0; +elseif isempty(direction); direction = 0; end; + +if ~exist('numperm'); numperm = 5000; +elseif isempty(numperm); numperm = 5000; +else if numperm < 1000; warning('You should run permutations more for the exact estimation (e.g. 1,000, 5,000, 10,000 times, or more.)'); end; end; + +if ~exist('thrClst'); thrClst = 0.05; +elseif isempty(thrClst); thrClst = 0.05; +else if thrClst >= 1; error('Cluster-wise threshold should be below one.'); end; end; + +%% Set the default parameters & get the information from input +pwr.dgr = 1; pwr.thr = 0; % two paramters described in the Discussion section of the original paper. This is for further development. +numSubj = size(setmat,3); numNode = size(setmat, 1); + +df = dbs_cal_df (hypoTest, numSubj); % calculate degree of freedom (DOF or DF) given a hypothesis and the number of subjects +icft = dbs_set_p_range (); % initial cluster-forming threshold (icft) +icft_temp = dbs_set_initrange (hypoTest, icft.p.range, df, direction); +icft.s.range =icft_temp; + +%% Test the hypothesis +[s, ~, p] = dbs_hypo_test(setmat, label, hypoTest, direction); + +%% Permutations for correction +tic; +[perm_result] = dbs_perm(setmat, label, numperm, hypoTest, direction, icft, pwr); + +thrClst1 = ceil(thrClst*numperm); thrClst2 = floor(thrClst*numperm); +for thre = 1 : length(icft.s.range) + dbs_thre.bd.thre{thre,1} = perm_result.bd.max_sort{thre}(thrClst2); + dbs_thre.wd.thre{thre,1} = perm_result.wd.max_sort{thre}(thrClst1); +end +temp_time = toc; +fprintf('\t( %.2f minutes elapsed. )\n', temp_time/60) + +%% Correction with DBS +dbs = dbs_correction(dbs_thre, icft, s, p); +dbs.bd.max_sort = perm_result.bd.max_sort; dbs.wd.max_sort = perm_result.wd.max_sort; + +%% Calculate Center Persistency (CP) score +cp_result = dbs_cal_cp (numNode, numperm, thrClst, icft.s.range, dbs, perm_result.cp, pwr); +dbs.cp = cp_result; + +%% Finish process +dbs.icft = icft; +dbs.time{1} = temp; dbs.time{2,1} = clock; +fprintf('\n\t\t%d, %d/%d, %d:%2.0f:%2.0f\tProcess done\n\n', dbs.time{2}) + diff --git a/dbs_perm.m b/dbs_perm.m new file mode 100644 index 0000000..479acfb --- /dev/null +++ b/dbs_perm.m @@ -0,0 +1,78 @@ +function [dbs_perm] = dbs_perm(setmat, label, numperm, hypotest, direction, icft, pwr) +% DBS_PERM Permutation to acquire an emprical null distributino of the maximum measures targeted. +% for FWE correction using DBS in connectivity analysis +% ================================================================================================================ +% [ INPUTS ] +% setmat = 3-D matrix. a set of matrices from multiple subjects. One for each. +% A size of setmat should be [N by N by M]. +% N: the number of nodes. +% M: the number of subjects +% +% label = 1-D vector. a list of labels with a value 0 or 1, indicating in which group each subject is included (for hypotest = 0 or 2) +% or with individual measures of all subjects for correlation (for hypotest = 1) +% +% numperm = # of permutations performed for family (default = 5000). +% +% hypotest = type of test (default = 0). +% 0: two-sample paired t-test (ttest) +% 1: correlation analysis (corr) +% 2: two-sample unpaired t-test (ttest2) +% +% icft = 1-D vector. a range of initial cluster-forming thresholds to test. +% Refer function 'DBS_main'. +% +% pwr = the power of the threshold and the degree similar to the one in TFCE. +% Refer function 'DBS_main'. +% ---------------------------------------------------------------------------------------------------------------- +% [ OUTPUTS ] +% dbs_perm +% ---------------------------------------------------------------------------------------------------------------- +% Last update: Aug 30, 2016. +% +% Copyright 2016. Kwangsun Yoo (K Yoo), PhD +% E-mail: rayksyoo@gmail.com / raybeam@kaist.ac.kr +% Laboratory for Cognitive Neuroscience and NeuroImaging (CNI) +% Department of Bio and Brain Engineering +% Korea Advanced Instititue of Science and Technology (KAIST) +% Daejeon, Republic of Korea +% ================================================================================================================ + +numsub = size(setmat,3); +numnode = size(setmat, 1); % or numnode = size(setmat,2); + +for i_perm = 1 : numperm + permlist(i_perm, :) = randperm(numsub); % permlist_uniq = unique(permlist,'rows'); +end + +t1 = clock; reverseStr = ''; h = waitbar(0, 'DBS is running ...'); + +for i_perm = 1 : numperm + temp_setmat = zeros(numnode, numnode, numsub ); + for j = 1:numsub + temp_setmat(:,:,j) = setmat(:,:,permlist(i_perm,j)); + end + [s, ~, p] = dbs_hypo_test(temp_setmat, label, hypotest, direction); + + for thre = 1:length(icft.s.range); clear a; + a = dbs_estm_dgr(p, s, icft.p.range(thre), icft.s.range(thre)); + dbs_perm.bd.org{thre,1} = a.bd.org; dbs_perm.bd.max{thre,1}(i_perm,1) = a.bd.max; + dbs_perm.wd.org{thre,1} = a.wd.org; dbs_perm.wd.max{thre,1}(i_perm,1) = a.wd.max; + end + + for thre_limit = 2 : length(icft.s.range) + cpMat = dbs_estm_cp(dbs_perm.wd.org, icft.s.range, thre_limit, numnode, pwr); + dbs_perm.cp.org{thre_limit,1} = cpMat; + dbs_perm.cp.max{thre_limit,1}(i_perm,1) = max(dbs_perm.cp.org{thre_limit}); + end + + [reverseStr, msg] = dbs_progress_box(i_perm, numperm, t1, reverseStr); + waitbar(i_perm/numperm, h, sprintf('DBS is running ... %0.1f %%', 100*i_perm/numperm)) +end +reverseStr = repmat(sprintf('\b'), 1, length(msg)); fprintf(reverseStr); close(h); + +for thre = 1:length(icft.s.range) + dbs_perm.bd.max_sort{thre,1} = sort(dbs_perm.bd.max{thre,1}, 'descend'); + dbs_perm.wd.max_sort{thre,1} = sort(dbs_perm.wd.max{thre,1}, 'descend'); +end + +fprintf('\t\t[ %d permutations done ]', i_perm) diff --git a/dbs_progress_box.m b/dbs_progress_box.m new file mode 100644 index 0000000..ed409b7 --- /dev/null +++ b/dbs_progress_box.m @@ -0,0 +1,25 @@ +function [reverseStr, msg] = dbs_progress_box(i, max, t1, reverseStr) +% DBS_PROGRESS_BOX +% ================================================================================================================ +% [ INPUTS ] +% i, max, t1, reverseStr +% ---------------------------------------------------------------------------------------------------------------- +% [ OUTPUTS ] +% reverseStr, msg +% ---------------------------------------------------------------------------------------------------------------- +% Last update: Aug 30, 2016. +% +% Copyright 2016. Kwangsun Yoo (K Yoo), PhD +% E-mail: rayksyoo@gmail.com / raybeam2kaist.ac.kr +% Laboratory for Cognitive Neuroscience and NeuroImaging (CNI) +% Department of Bio and Brain Engineering +% Korea Advanced Instititue of Science and Technology (KAIST) +% Daejeon, Republic of Korea +% ================================================================================================================ + +t2 = clock; t3 = (etime(t2,t1)/i) * (max-i); +perc = 100*i/max; + +msg = sprintf('\t\t[ Running %d permutations >>> %.1f percent done, %.0f seconds left. ]', max, perc, t3); +fprintf([reverseStr, msg]) +reverseStr = repmat(sprintf('\b'), 1, length(msg)); \ No newline at end of file diff --git a/dbs_set_initrange.m b/dbs_set_initrange.m new file mode 100644 index 0000000..a06c413 --- /dev/null +++ b/dbs_set_initrange.m @@ -0,0 +1,59 @@ +function stats = dbs_set_initrange (hypotest, range, df, direction) +% DBS_SET_INITRANGE Set a range of the initial cluster-forming threshold statistical value +% (corresponding to the p-value) to be tested +% ================================================================================================================ +% [ INPUTS ] +% hypotest = type of test (default = 0). +% 0: two-sample paired t-test (ttest) +% 1: two-sample unpaired t-test (ttest2) (assumping the same variance for the two groups) +% 2: correlation analysis (corr) +% +% range +% +% df = the degree of freedom (DOF) for the hypothesis testing and the number of samples + +% direction +% 0: g1 = g2 (two-tail) +% 1: g1 > g2 (one-tail) +% -1: g2 < g1 (one-tail) +% ---------------------------------------------------------------------------------------------------------------- +% [ OUTPUTS ] +% stats +% ---------------------------------------------------------------------------------------------------------------- +% Last update: Aug 30, 2016. +% +% Copyright 2016. Kwangsun Yoo (K Yoo), PhD +% E-mail: rayksyoo@gmail.com / raybeam@kaist.ac.kr +% Laboratory for Cognitive Neuroscience and NeuroImaging (CNI) +% Department of Bio and Brain Engineering +% Korea Advanced Instititue of Science and Technology (KAIST) +% Daejeon, Republic of Korea +% ================================================================================================================ + +if hypotest == 0 + if direction == 0 + t = tinv( range/2 , df ); % for two-tailed (subdivision by 2 represents the two-tail) paired t-test. If <0, 1st grp has higher. + elseif direction == 1 %% testing whether [g1 > g2] (the same as ttest right) + t = tinv( 1 - range , df ); + elseif direction == - 1 %% testing whether [g1 < g2] (the same as ttest left) + t = tinv( range , df ); + end; stats = t; +elseif hypotest == 1 + if direction == 0 + t = tinv( range/2 , df ); % for unpaired t-test df = (numfirst-1) + (numsub-numfirst -1) = n_total - 2 for the equal variances + elseif direction == 1 %% testing whether [g1 > g2] (the same as ttest right) + t = tinv( 1- range , df ); + elseif direction == - 1 %% testing whether [g1 < g2] (the same as ttest left) + t = tinv( range , df ); + end; stats = t; +elseif hypotest == 2 + if direction == 0 + t = tinv( range/2 , df ); % for two-tailed (subdivision by 2 represents the two-tail) paired t-test. + elseif direction == 1 %% testing whether positive corrlation (the same as corr right) (with ttest right) + t = tinv( 1- range , df ); + elseif direction == - 1 %% testing whether negative corrlation (the same as corr left) (with ttest left) + t = tinv( range , df ); + end; r = sqrt( t.^2 ./ (df + t.^2)); % t to r transformation (from r to t transform) + stats = r; +end; +stats = abs(stats); diff --git a/dbs_set_p_range.m b/dbs_set_p_range.m new file mode 100644 index 0000000..ae308ce --- /dev/null +++ b/dbs_set_p_range.m @@ -0,0 +1,28 @@ +function a = dbs_set_p_range +% DBS_SET_P_RANGE Set a range of the initial cluster-forming threshold p-value to be tested +% ================================================================================================================ +% [ OUTPUTS ] +% a.p.lb = lower boundary of the range, p-value = 0.05. +% a.p.hb = higher boundary of the range, p-value = 0.0001. +% a.p.range = the range between a.p.lb and a.p.hb +% ---------------------------------------------------------------------------------------------------------------- +% Last update: Aug 30, 2016. +% +% Copyright 2016. Kwangsun Yoo (K Yoo), PhD +% E-mail: rayksyoo@gmail.com / raybeam@kaist.ac.kr +% Laboratory for Cognitive Neuroscience and NeuroImaging (CNI) +% Department of Bio and Brain Engineering +% Korea Advanced Instititue of Science and Technology (KAIST) +% Daejeon, Republic of Korea +% ================================================================================================================ + +a.p.lb = 0.05; % lower boundary +a.p.hb = 0.0001; % higher boundary +log10p.lb = log10(a.p.lb); +log10p.hb = log10(a.p.hb); +log10p.range = [ log10(a.p.lb) : -0.1 : log10(a.p.hb)-0.1 ]'; + +temp1 = 10000 * (10.^(log10p.range)); +temp2 = 10000 * [ a.p.lb : -0.005 : a.p.hb ]'; + +a.p.range = flipud(union(round(temp1), round(temp2)))/10000; diff --git a/hypo_test/dbs_corr_matrix.m b/hypo_test/dbs_corr_matrix.m new file mode 100644 index 0000000..e19abc0 --- /dev/null +++ b/hypo_test/dbs_corr_matrix.m @@ -0,0 +1,53 @@ +function [r, t, p, df] = dbs_corr_matrix (gr , label , direction) +% DBS_CORR_MATRIX +% ================================================================================================================ +% [ INPUTS ] +% gr = 3-D matrix which consists of a set of 2-D matrices from multiple subjects. +% A size of setmat should be [N by N by M]. +% N: the number of nodes. +% M: the number of subjects. +% +% label = 1-D vector containing a list of labels +% with individual measures of behavioral performance from each subjects for correlation +% +% direction +% 0: g1 = g2 (two-tail) +% 1: g1 > g2 (one-tail) +% -1: g2 < g1 (one-tail) +% ---------------------------------------------------------------------------------------------------------------- +% [ OUTPUTS ] +% r, t, p, df +% ---------------------------------------------------------------------------------------------------------------- +% Last update: Aug 30, 2016. +% +% Copyright 2016. Kwangsun Yoo (K Yoo), PhD +% E-mail: rayksyoo@gmail.com / raybeam@kaist.ac.kr +% Laboratory for Cognitive Neuroscience and NeuroImaging (CNI) +% Department of Bio and Brain Engineering +% Korea Advanced Instititue of Science and Technology (KAIST) +% Daejeon, Republic of Korea +% ================================================================================================================ + +setmat.mean = mean(gr, 3); +numsub = size(gr,3); +df = numsub - 2; + +lab.B = zeros(1,1, numsub); +lab.B(1,1,:) = label; +lab.mean = mean(lab.B,3); +lab.full = repmat(lab.B, [size(gr,1), size(gr,2)]); + +lab.mean2 = mean(lab.full, 3); +setmat.mean2 = repmat(setmat.mean, [1, 1, numsub]); +lab.B2 = repmat(lab.B, [1, 1, numsub]); +lab.mean3 = repmat(lab.mean2, [1, 1, numsub]); + +r = sum ( (gr - setmat.mean2) .* (lab.full - lab.mean3), 3) ./ sqrt ( sum( (gr - setmat.mean2).^2 ,3) .* sum( (lab.full - lab.mean3).^2 ,3) ); +t = r ./ sqrt( (1-r.^2) / df) ; +if direction == 0 + p = 2 * tcdf(abs(t), df, 'upper' ) ; +elseif direction == 1 %% testing whether positive corrlation (the same as corr right) (with ttest right) + p = tcdf(t, df, 'upper' ) ; +elseif direction == -1 %% testing whether negative corrlation (the same as corr left) (with ttest left) + p = tcdf(t, df ) ; +end \ No newline at end of file diff --git a/hypo_test/dbs_pairT_matrix.m b/hypo_test/dbs_pairT_matrix.m new file mode 100644 index 0000000..f413f2e --- /dev/null +++ b/hypo_test/dbs_pairT_matrix.m @@ -0,0 +1,51 @@ +function [t, p, df] = dbs_pairT_matrix (gr1 , gr2 , direction) +% DBS_PAIRT_MATRIX +% ================================================================================================================ +% [ INPUTS ] +% gr1, gr2 = a pair of 3-D matrix which consists of a set of 2-D matrices from multiple subjects. +% A size of setmat should be [N by N by M]. +% N: the number of nodes. +% M: the number of subjects. the same across two groups +% +% direction +% 0: g1 = g2 (two-tail) +% 1: g1 > g2 (one-tail) +% -1: g2 < g1 (one-tail) +% ---------------------------------------------------------------------------------------------------------------- +% [ OUTPUTS ] +% t, p, df +% ---------------------------------------------------------------------------------------------------------------- +% Last update: Aug 30, 2016. +% +% Copyright 2016. Kwangsun Yoo (K Yoo), PhD +% E-mail: rayksyoo@gmail.com / raybeam@kaist.ac.kr +% Laboratory for Cognitive Neuroscience and NeuroImaging (CNI) +% Department of Bio and Brain Engineering +% Korea Advanced Instititue of Science and Technology (KAIST) +% Daejeon, Republic of Korea +% ================================================================================================================ + +mu_zero = 0; % target difference to test + +setmat.diff = gr1 - gr2; +setmat.mean = mean(setmat.diff,3); + +numfirst = size(gr1,3); +df = numfirst - 1; + +setmat.mean2 = repmat(setmat.mean, [1, 1, numfirst]); + +setmat.csd = sqrt ( sum((setmat.diff - setmat.mean2).^2, 3) ./ df ); % csd: corrected sample standard deviation +t = (setmat.mean - mu_zero) ./ (setmat.csd / sqrt(numfirst)); + +if direction == 0 + p = 2 * tcdf(abs(t), df, 'upper') ; +elseif direction == 1 %% testing whether [g1 > g2] (the same as ttest right) + p = tcdf(t, df, 'upper') ; +elseif direction == -1 %% testing whether [g1 < g2] (the same as ttest left) + p = tcdf(t, df) ; +end + + + + \ No newline at end of file diff --git a/hypo_test/dbs_unpairT_matrix.m b/hypo_test/dbs_unpairT_matrix.m new file mode 100644 index 0000000..9c770c6 --- /dev/null +++ b/hypo_test/dbs_unpairT_matrix.m @@ -0,0 +1,48 @@ +function [t, p, df] = dbs_unpairT_matrix (gr1 , gr2 , direction) +% DBS_UNPAIRT_MATRIX +% ================================================================================================================ +% [ INPUTS ] +% gr1, gr2 = a pair of 3-D matrix which consists of a set of 2-D matrices from multiple subjects. +% A size of setmat should be [N by N by M]. +% N: the number of nodes. +% M: the number of subjects. Two groups could have different M. +% +% direction +% 0: g1 = g2 (two-tail) +% 1: g1 > g2 (one-tail) +% -1: g2 < g1 (one-tail) +% ---------------------------------------------------------------------------------------------------------------- +% [ OUTPUTS ] +% t, p, df +% ---------------------------------------------------------------------------------------------------------------- +% Last update: Aug 30, 2016. +% +% Copyright 2016. Kwangsun Yoo (K Yoo), PhD +% E-mail: rayksyoo@gmail.com / raybeam@kaist.ac.kr +% Laboratory for Cognitive Neuroscience and NeuroImaging (CNI) +% Department of Bio and Brain Engineering +% Korea Advanced Instititue of Science and Technology (KAIST) +% Daejeon, Republic of Korea +% ================================================================================================================ + +setmat.g1.mean = mean(gr1, 3); +setmat.g2.mean = mean(gr2, 3); + +numfirst = size(gr1,3); +numsecond = size(gr2,3); +df = numfirst+numsecond-2; + +setmat.g1.mean2 = repmat(setmat.g1.mean, [1, 1, numfirst]); +setmat.g2.mean2 = repmat(setmat.g2.mean, [1, 1, numsecond]); + +setmat.g1.usv = sum( (gr1 - setmat.g1.mean2) .^2 ,3) ./ (numfirst-1); % usv: unbiased sample variation +setmat.g2.usv = sum( (gr2 - setmat.g2.mean2) .^2 ,3) ./ (numsecond-1); + +t = ( setmat.g1.mean - setmat.g2.mean) ./ sqrt ( ( ( (numfirst-1) * setmat.g1.usv + (numsecond-1) * setmat.g2.usv ) / df) * ( (1/numfirst) + (1/numsecond) ) ); +if direction == 0 + p = 2 * tcdf(abs(t), df, 'upper' ) ; +elseif direction == 1 %% testing whether [g1 > g2] (the same as ttest right) + p = tcdf(t, df, 'upper' ) ; +elseif direction == -1 %% testing whether [g1 < g2] (the same as ttest left) + p = tcdf(t, df) ; +end