00001 %GroupAna.m
00002 %
00003 %Purpose:
00004 %
00005 %
00006 %
00007 %Input Parameters:
00008 %
00009 %
00010 %
00011 %Output Parameters:
00012 % err : 0 No Problem
00013 % : 1 Problems
00014 %
00015 %
00016 %
00017 %Key Terms:
00018 %
00019 %More Info :
00020 %
00021 %
00022 %
00023 %
00024 % Author : Gang Chen
00025 % Date : Tue Mar 23 14:06:44 EST 2004
00026 % SSCC/NIMH/ National Institutes of Health, Bethesda MD 20892
00027
00028
00029 %Define the function name for easy referencing
00030 FuncName = 'GroupAna.m';
00031
00032 %Debug Flag
00033 DBG = 1;
00034
00035 %initailize return variables
00036 err = 1;
00037
00038 clear all;
00039
00040 diary('diary'); % save all subsequent command window input and most of the resulting command window output to be appended
00041 % to file 'diary'
00042 fprintf('\nPlease read the following carefully about group analysis setup:\n');
00043 fprintf('\n\t1. If the resolution of your EPI data is near n millimeters, during Talairach conversion use\n');
00044 fprintf('\t "the command adwarp -dxyz n" to greatly reduce runtime without sacrificing accuracy.\n');
00045 fprintf('\n\t2. We strongly suggest that factor names be labeled with short (2 or 3 capital letters) names\n');
00046 fprintf('\t so that subbrik labels can be shown on AFNI viewer.\n');
00047 fprintf('\n\t3. With nesting, arrange your design in such a way that the last factor is nested within the 1st factor.\n');
00048 %fprintf('\n\t4. The program can only handle balanced design now. \n');
00049 fprintf('\n\t4. Each input file should include only one subbrik. We suggest files be \n');
00050 fprintf('\t named by reflecting the hierarchy of the experiment design.\n');
00051 fprintf('\n\t5. Currently all of the following terms are modeled: main effects and applicable interactions in various orders, .\n');
00052 fprintf('\n\t6. One covariate is currently allowed in the analysis, which should be in the format of one-column text file.\n');
00053 fprintf('\t The column length has to be the same as the total number of input files.\n');
00054
00055
00056 % Grouop analysis for Volume or Surface data?
00057 flg = 0;
00058 while flg == 0,
00059 data_type = input('\nGroups analysis for volume or surface data (0 - volume; 1 - surface)? ');
00060 if (data_type ~= 0 & data_type ~= 1),
00061 flg = 0; fprintf(2,'Error: wrong input! Please try it again.\n');
00062 else flg = 1;
00063 end
00064 end
00065
00066 % If for surface data, acquire number of nodes
00067 flg = 0;
00068 if (data_type == 1),
00069 fprintf(1, '\nProgram @Purify_1D can be used to extract the regressor coefficient column to a 1D file.');
00070 Frame_N = input('\nWhich column corresponds to regressor coefficient? (1, 2, 3, ...) ');
00071 while flg == 0,
00072 node_n = input('\nHow many number of nodes in the surface data? ');
00073 if (isnumeric(node_n) == 0 | isempty(node_n)),
00074 flg = 0; fprintf(2,'Error: the input is not a number. Please try it again.\n');
00075 else flg = 1; end
00076 end
00077 else Opt.Frames = 1; % In the case of volumetric data, it's supposed to have only ONE subbrik!
00078 end
00079
00080 flg = 0;
00081 while flg == 0,
00082 NF = input('\nHow many factors? ');
00083 if (isnumeric(NF) == 0 | isempty(NF)),
00084 flg = 0; fprintf(2,'Error: the input is not a number. Please try it again.\n');
00085 else flg = 1;
00086 end
00087 end
00088
00089 switch NF
00090 case 1,
00091 fprintf('\nAvailable design types: \n');
00092 fprintf('\n\tType 1: one factor with fixed effect.\n');
00093
00094 case 2,
00095 fprintf('\nAvailabe design types:\n');
00096 fprintf('\n\tType 1: Factorial (crossed) design AXB - both factors are fixed.\n');
00097 fprintf('\n\tType 2: Factorial (crossed) design AXB - factor A is fixed while B is random. If B is subject, it is');
00098 fprintf('\n\t usually called 1-way design with A varying within subject. Notice: It is inappropriate to');
00099 fprintf('\n\t run any constrasts including mean estimates and differences for factor B with this design type.\n');
00100 fprintf('\n\tType 3: Factorial (crossed) design AXB - both factors are random.\n');
00101
00102 case 3,
00103 fprintf('\nAvailabe design types:\n');
00104 fprintf('\n\tType 1: Factorial (crossed) design AXBXC - all factors are fixed.\n');
00105 fprintf('\n\tType 2: Factorial (crossed) design AXBXC - factors A and B are fixed while C is random. If C is subject, it is');
00106 fprintf('\n\t usually called 2-way design with A and B varying within subject. Notice: It is inappropriate to');
00107 fprintf('\n\t run any constrast tests including mean estimates and differences for factor C with this design type.\n');
00108 fprintf('\n\tType 3: Mixed design BXC(A) - A and B are fixed while C (usually subject) is random and nested within A.\n');
00109 fprintf('\n\tType 4: Mixed design BXC(A) - Fixed factor C is nested within fixed factor A while B (usually subject) is random.\n');
00110
00111 case 4,
00112 fprintf('\nAvailabe design types:\n');
00113 fprintf('\n\tType 1: Factorial (crossed) design AXBXCXD - all 4 factors are fixed.\n');
00114 fprintf('\n\tType 2: Factorial (crossed) design AXBXCXD - only factor D is random. If D is subject it is also');
00115 fprintf('\n\t called 3-way design with all 3 factors A, B, and C varying within subject.');
00116 fprintf('\n\tType 3: Mixed design BXCXD(A)- only the nested (4th) factor D (usually subject) is random.');
00117 fprintf('\n\t Also called 3-way design with factors B and C varying within subject and factor A between subjects.');
00118 fprintf('\n\tType 4: Mixed design BXCXD(A)- D is nested within A, but only the 3rd factor (usually subject)');
00119 fprintf('\n\t is random.');
00120 fprintf('\n\tType 5: Mixed design CXD(AXB) - only the nested (4th) factor D (usually subject) is random,');
00121 fprintf('\n\t but factor D is nested within both factors A and B. If D is subject it is also called 3-way');
00122 fprintf('\n\t design with factor D varying within-subject and factors A and B between-subjects.\n');
00123 % fprintf('\nNotice: This is NOT an exhaustive list of design types for 4-way ANOVA. Other types might be implemented upon request.\n');
00124
00125 case 5,
00126 fprintf('\nAvailabe design types:\n');
00127 fprintf('\n\tType 1: Factorial (crossed) design AXBXCXDXE - all 5 factors are fixed.\n');
00128 fprintf('\n\tType 2: Factorial (crossed) design AXBXCXDXE - only factor E is random. If E is subject it is also');
00129 fprintf('\n\t called 4-way design with all 4 factors A, B, C and D varying within subject.\n');
00130 fprintf('\n\tType 3: Mixed design BXCXDXE(A) - only the nested (5th) factor E (usually subject) is random.');
00131 fprintf('\n\t Also called 4-way design with factors B, C and D varying within subject and factor A between subjects.\n');
00132 fprintf('\n\tType 4: Mixed design BXCXDXE(A) - the 5th factor E is nested within factor A, but factor D (usually subject)');
00133 fprintf('\n\t is random.\n');
00134
00135 end
00136
00137 flg = 0;
00138 while flg == 0,
00139 dsgn = input('\nChoose design type (1, 2, 3, 4, 5): ');
00140 if (isnumeric(dsgn) == 0 | isempty(dsgn)),
00141 flg = 0; fprintf(2,'Error: input is not a number. Please try it again.\n');
00142 else flg = 1;
00143 end
00144 end
00145
00146 %flg = 0;
00147 %while flg == 0,
00148 % slices = input('How many slices along the Z axis (run 3dinfo on one of the input files to find out)? ');
00149 % if (isnumeric(slices) == 0 | isempty(slices)),
00150 % flg = 0; fprintf(2,'Error: the input is not a number. Please try it again.\n');
00151 % else flg = 1;
00152 % end
00153 %end
00154
00155 % Balanced design?
00156 flg = 0;
00157 while flg == 0,
00158
00159 unbalanced.yes = 1 - input('\nIs the design balanced? (1 - Yes; 0 - No) ');
00160 if (unbalanced.yes ~= 0 & unbalanced.yes ~= 1),
00161 flg = 0; fprintf(2,'Error: inapproriate input. Please try it again.\n');
00162 else flg = 1;
00163 end
00164
00165 if (unbalanced.yes == 1),
00166 fprintf('\nThe following two kinds of unbalanced designs are currently allowed:\n')
00167 fprintf('\n(1) All factors are fixed - ');
00168 fprintf('\n\t 1-way ANOVA; 2-way ANOVA type 1: AXB; and 3-way ANOVA type 1: AXBXC');
00169 fprintf('\n\n(2) When a random factor (subject) is nested within another factor A,');
00170 fprintf('\n\t each level of factor A contains a unique and unequal number of subjects - ');
00171 fprintf('\n\t 3-way ANOVA type 3: BXC(A); and 4-way ANOVA type 3: BXCXD(A)')
00172 if (input('\n\nDoes your unbalanced design belong to either of the above types? (1 - Yes; 0 - No) ') == 0);
00173 while (1); fprintf(2,'Halted: Ctrl+c to exit'); pause; end
00174 end
00175 end
00176
00177 % fprintf('\nTwo kinds of unbalanced design are considered:\n');
00178 % fprintf('\nType 1: When subject is a random factor and nested within another factor (3-way ANOVA type 3, 4-way ANOVA types 3 and 5),');
00179 % fprintf('\n\t there can be different number of subjects for each level of factor A;\n');
00180 % fprintf('\nType 2: Not type 1, but there are different sample sizes for some cells (factor level combinations).');
00181 % fprintf('\n\t For example, add somethint here later...\n');
00182 %
00183 % unbalanced.type = input('\nChoose design type (1, 2): ');
00184 % if (unbalanced.type ~= 1 & unbalanced.type ~= 2),
00185 % flg = 0; fprintf(2,'Error: inapproriate input. Please try it again.\n');
00186 % else flg = 1;
00187 % end
00188 end
00189
00190 for (ii = 1:1:5)
00191 FL(ii).N_level = 0; % initialization for ContrVec.m so that their values are available for lower-way ANOVA
00192 end
00193
00194
00195 %Get the number of levels for each factor and also the total number of input files
00196
00197 flg = 0;
00198 ntot = 1;
00199
00200 if (unbalanced.yes == 1),
00201
00202 if ((NF == 1 | NF == 2 | NF == 3 | NF == 4) & dsgn == 1), % Basically for 1,2,3,4-way ANCOVA: This loop is the same as balanced. Maybe also for 4-way ANCOVA?
00203 for (i=1:1:NF),
00204 fprintf('\nLabel for No. %i ', i);
00205 FL(i).expr = input('factor: ', 's'); % Factor Label i
00206 fprintf(2,'How many levels does factor %c (%s) ', 64+i, FL(i).expr);
00207 FL(i).N_level = input('have? ');
00208 if (isnumeric(FL(i).N_level) == 0 | isempty(FL(i).N_level)),
00209 flg = 0; fprintf(2,'Error: the input is not a number. Please try again.\n');
00210 else flg = 1; end
00211 for (j=1:1:FL(i).N_level),
00212 fprintf('Label for No. %i level of factor %c (%s)', j, 64+i, FL(i).expr);
00213 FL(i).level(j).expr = input(' is: ', 's');
00214 end
00215 sz(i) = FL(i).N_level; % number of levels of factor i
00216 ntot = ntot * FL(i).N_level; %total number of combinations
00217 end
00218 end
00219
00220 if (((NF == 3 | NF == 4) & dsgn == 3)),
00221 % if (unbalanced.type == 1),
00222 for (i=1:1:(NF-1)),
00223 fprintf('\nLabel for No. %i ', i);
00224 FL(i).expr = input('factor: ', 's'); % Factor Label i
00225 fprintf(2,'How many levels does factor %c (%s) ', 64+i, FL(i).expr);
00226 FL(i).N_level = input('have? ');
00227 if (isnumeric(FL(i).N_level) == 0 | isempty(FL(i).N_level)),
00228 flg = 0; fprintf(2,'Error: the input is not a number. Please try again.\n');
00229 else flg = 1;
00230 end
00231 for (j=1:1:FL(i).N_level),
00232 fprintf('Label for No. %i level of factor %c (%s)', j, 64+i, FL(i).expr);
00233 FL(i).level(j).expr = input(' is: ', 's');
00234 end
00235 sz(i) = FL(i).N_level; % number of levels of factor i
00236 ntot = ntot * FL(i).N_level; %total number of combinations
00237 end % for (i=1:1:(NF-1))
00238
00239 FL(NF).N_level = 0;
00240 fprintf(2, '\nLabel for No. %i ', NF);
00241 FL(NF).expr = input('factor: ', 's'); % Label this unbalanced factor
00242 % fprintf(2, '\nNote: Since this is a nested design the label for levels (usuall subject names) of factor No. %i (%c - %s)', NF, 64+NF, FL(NF).expr);
00243 % fprintf(2, '\nhas to be DIFFERENT for each level of factor %c (%s)!!!\n\n', 64+1, FL(1).expr);
00244
00245 flag = 0;
00246 while flag == 0,
00247 combine = [];
00248 for (i = 1:1:FL(1).N_level),
00249 fprintf(2,'How many levels does factor %c (%s) corresponding to level %i (%s) of factor %c (%s) ', ...
00250 64+NF, FL(NF).expr, i, FL(1).level(i).expr, 64+1, FL(1).expr);
00251 FL(NF).UL(i).N_level = input('have? '); % unbalanced levels for this factor
00252 for (j=1:1:FL(NF).UL(i).N_level),
00253 fprintf('Label for No. %i level of factor %c (%s) in group %i of factor %c (%s)', j, 64+NF, FL(NF).expr, i, 64+1, FL(1).expr);
00254 FL(NF).UL(i).n(j).expr = input(' is: ', 's');
00255 combine = [combine {FL(NF).UL(i).n(j).expr}]; % Concatenate them to make a cell array
00256 end
00257 % FL(NF).N_level = FL(NF).N_level + FL(NF).UL(i).N_level;
00258 % FL(NF).N_level = max([FL(NF).UL(:).N_level]); % This is for positioning those contrast columns in the design matrix in PreProc.m
00259 end
00260
00261 if (length(unique(combine)) == max([FL(NF).UL(:).N_level])), % if same labels are used across groups
00262 FL(NF).N_level = max([FL(NF).UL(:).N_level]); flag = 1; % design matrix is built based on the longest group length
00263 ntot = ntot * sum([FL(NF).UL(:).N_level]) / FL(1).N_level;
00264 elseif (length(unique(combine)) == sum([FL(NF).UL(:).N_level])), % if different labels are used across groups
00265 FL(NF).N_level = sum([FL(NF).UL(:).N_level]); flag = 1; % design matrix is built based on the total length of the groups
00266 ntot = ntot * FL(NF).N_level / FL(1).N_level;
00267 else
00268 flag = 0;
00269 fprintf(2, '\nError: There is some overlap among the labels of factor %c (%s) across the groups of factor %c (%s)!\n', ...
00270 64+NF, FL(NF).expr, 64+1, FL(1).expr);
00271 end
00272 end % while flag == 0
00273
00274 end % if (((NF == 3 | NF == 4) & dsgn == 3))
00275
00276 else % Balanced designs
00277
00278 for (i=1:1:NF),
00279 fprintf('\nLabel for No. %i ', i);
00280 FL(i).expr = input('factor: ', 's'); % Factor Label i
00281 fprintf(2,'How many levels does factor %c (%s) ', 64+i, FL(i).expr);
00282 FL(i).N_level = input('have? ');
00283 if (isnumeric(FL(i).N_level) == 0 | isempty(FL(i).N_level)),
00284 flg = 0; fprintf(2,'Error: the input is not a number. Please try again.\n');
00285 else flg = 1;
00286 end
00287 for (j=1:1:FL(i).N_level),
00288 fprintf('Label for No. %i level of factor %c (%s)', j, 64+i, FL(i).expr);
00289 FL(i).level(j).expr = input(' is: ', 's');
00290 end
00291 sz(i) = FL(i).N_level; % number of levels of factor i
00292 ntot = ntot * FL(i).N_level; %total number of combinations
00293 end
00294
00295 end % if (unbalanced)
00296
00297 if ~((NF == 1 | NF == 2 | NF == 3 | NF == 4) & dsgn == 1), % Mainly for ANCOVA?
00298 fprintf(2,'\nEnter the sample size (number of observations) per cell');
00299 FL(NF+1).N_level = input(': '); % Should it be changed to a different name instead of FL???
00300 ntot = ntot * FL(NF+1).N_level; % total number of factor combinations including repeats
00301 sz(NF+1) = FL(NF+1).N_level;
00302 end
00303
00304 % Running ANOCVA?
00305 flg = 0;
00306 %cov.label = [];
00307 while flg == 0,
00308 cov.do = input('\nAny covariate (concomitant variable)? (1 - Yes; 0 - No) ');
00309 if (cov.do ~= 0 & cov.do ~= 1),
00310 flg = 0; fprintf(2,'Error: inapproriate input. Please try it again.\n');
00311 else flg = 1;
00312 end
00313 end
00314
00315 %Only allow one covariate now!
00316 if (cov.do),
00317
00318 FL(NF+1).N_level = 1; % This is for PreProc.m to get the correct 'shift' position
00319
00320 if (NF == 2 & dsgn == 3), % 2-way ANOVA of A(R)XB(R) possible in FMRI?
00321 fprintf('\nThe 2-way ANCOVA for this design is currently NOT available.\n');
00322 while (1); fprintf(2,'Halted: Ctrl+c to exit'); pause; end
00323 end
00324
00325 flg = 0;
00326 cov.label = input('Label for the covariate is: ', 's');
00327 fprintf('The 1D file for the covariate has to be in the format of one-column text.\n');
00328 fprintf('And the column should have exactly the same number and order of those input files.\n');
00329 while flg == 0,
00330 cov.FN = input('\nConvariate file name: ', 's');
00331 fid0 = fopen(cov.FN,'r');
00332 if (fid0 == -1),
00333 flg = 0; fprintf(2,'Error: File %s does not exist. Please try it again. \n', cov.FN);
00334 else flg = 1;
00335 [cov.vec, count] = fscanf(fid0, '%f');
00336 if (count ~= ntot & ~((NF ==1 | NF == 2 | NF == 3) & dsgn == 1)), % Check length of the 1D file
00337 fprintf(2, '\nError: The column length of the covariate has to equal to the total number of input files!\n');
00338 fprintf(2,'Halted: Ctrl+c to exit'); pause;
00339 end
00340 end
00341 end % while
00342 end
00343
00344 %flg = 0;
00345 %while flg == 0,
00346 % fprintf('\nIs input in the format of one-subbrik-one-file or one-file-per-subject?\n');
00347 % file_format = input('(1 - Single brik; 0 - Multiple subbriks) ');
00348 % if (file_format ~= 0 & file_format ~= 1),
00349 % flg = 0; fprintf(2,'Error: inapproriate input. Please try it again.\n');
00350 % else flg = 1;
00351 % end
00352 %end
00353
00354 fprintf('\nAll input files are supposed to contain only one subbrik.\n');
00355 file_format = 1; % Ignore that file_format = 0 now since it leads to too much trouble
00356
00357 if (file_format == 0), % unused now!
00358 flg = 0;
00359 file_num = input('Number of input files: ');
00360 file_SB = input('Number of subbriks for analysis in each file: ');
00361 if (isnumeric(file_num) == 0 | isempty(file_num) | isnumeric(file_SB) == 0 | isempty(file_SB)),
00362 flg = 0; fprintf(2,'Error: the input is not a number. Please try it again.\n');
00363 else flg = 1; end
00364 if (file_num*file_SB ~= ntot),
00365 fprintf(2, '\nError: The number of files and subbriks are not consistent!\n');
00366 fprintf(2,'Halted: Ctrl+c to exit'); pause;
00367 end
00368 end % if (file_format == 0)
00369
00370 if (file_format == 1 & ~((NF == 1 | NF == 2 | NF == 3 | NF == 4) & dsgn == 1)),
00371 fprintf(2,'\nThere should be totally %i input files. \n', ntot);
00372 corr = input('Correct? (1 - Yes; 0 - No) ');
00373 if (corr == 0),
00374 fprintf(2,'Error: Check the inconsistency. \n');
00375 while (1); fprintf(2,'Halted: Ctrl+c to exit'); pause; end
00376 end
00377 end
00378
00379 %generate the subcripts and obtain all input files
00380 GP = cell(NF, ntot); %creat a cell array to reflect the structure of all the combinations
00381
00382 if (file_format == 1),
00383
00384 if (unbalanced.yes == 1), % Try 4-way design 3 first
00385
00386 % if ((NF == 1 | NF == 2 | NF == 3 | NF == 4) & dsgn == 1), % Meant for 1,2,3,4-way ANCOVA with unequal sample size
00387 if (dsgn == 1),
00388 FI = 0; % File index
00389 flg = 0;
00390 ntot = input('Total number of input files: ');
00391 if (isnumeric(ntot) == 0 | isempty(ntot)),
00392 flg = 0; fprintf(2,'Error: the input is not a number. Please try again.\n');
00393 else flg = 1;
00394 end
00395
00396 switch NF
00397 case 1,
00398 for (ii = 1:1:FL(1).N_level),
00399 fprintf (2,'\nFor factor %c (%s) at level %i (%s),', 64+1, FL(1).expr, ii, FL(1).level(ii).expr);
00400 sz = input('\nsample size is: ');
00401 fprintf (2,'\nProvide those %i input files:\n', sz);
00402 for (rr = 1:1:sz),
00403 FI = FI + 1;
00404 GP(1, FI) = {FL(1).level(ii).expr};
00405
00406 flg = 0;
00407 while flg == 0,
00408 fprintf (2,'No. %i file ', rr);
00409 file(FI).nm = input('is: ', 's');
00410 fid = fopen (file(FI).nm,'r');
00411 if (fid == -1), flg = 0; fprintf(2,'Error: File %s does not exist. Please try it again. \n', file(FI).nm);
00412 else flg = 1; fclose (fid); end
00413 if isempty(strfind(file(FI).nm, 'tlrc')) == 0
00414 format = 'tlrc';
00415 elseif isempty(strfind(file(FI).nm, 'orig')) == 0
00416 format = 'orig';
00417 else
00418 if isempty(strfind(file(FI).nm, '1D')) == 0 % if 1D file (surface data)
00419 format = '1D';
00420 else
00421 while (1); fprintf(2,'Error: format of file %s is incorrect!\n', file(FI).nm);
00422 fprintf(2,'Halted: Ctrl+c to exit'); pause; end
00423 end
00424 end
00425 end
00426
00427 end % for (rr = 1:1:sz)
00428 end % for (ii = 1:1:FL(1).N_level)
00429
00430 case 2,
00431 for (ii = 1:1:FL(1).N_level),
00432 for (jj = 1:1:FL(2).N_level),
00433 fprintf (2,'\nFor factor %c (%s) at level %i (%s),', 64+1, FL(1).expr, ii, FL(1).level(ii).expr);
00434 fprintf (2,'\n factor %c (%s) at level %i (%s),', 64+2, FL(2).expr, jj, FL(2).level(jj).expr);
00435 sz = input('\nsample size is: ');
00436 fprintf (2,'\nProvide those %i input files:\n', sz);
00437 for (rr = 1:1:sz),
00438 FI = FI + 1;
00439 GP(1, FI) = {FL(1).level(ii).expr};
00440 GP(2, FI) = {FL(2).level(jj).expr};
00441
00442 flg = 0;
00443 while flg == 0,
00444 fprintf (2,'No. %i file ', rr);
00445 file(FI).nm = input('is: ', 's');
00446 fid = fopen (file(FI).nm,'r');
00447 if (fid == -1), flg = 0; fprintf(2,'Error: File %s does not exist. Please try it again. \n', file(FI).nm);
00448 else flg = 1; fclose (fid); end
00449 if isempty(strfind(file(FI).nm, 'tlrc')) == 0
00450 format = 'tlrc';
00451 elseif isempty(strfind(file(FI).nm, 'orig')) == 0
00452 format = 'orig';
00453 else
00454 if isempty(strfind(file(FI).nm, '1D')) == 0 % if 1D file (surface data)
00455 format = '1D';
00456 else
00457 while (1); fprintf(2,'Error: format of file %s is incorrect!\n', file(FI).nm);
00458 fprintf(2,'Halted: Ctrl+c to exit'); pause; end
00459 end
00460 end
00461 end
00462
00463 end % for (rr = 1:1:sz)
00464 end % for (jj = 1:1:FL(2).N_level)
00465 end % for (ii = 1:1:FL(1).N_level)
00466
00467 case 3,
00468 for (ii = 1:1:FL(1).N_level),
00469 for (jj = 1:1:FL(2).N_level),
00470 for (kk = 1:1:FL(3).N_level),
00471 fprintf (2,'\nFor factor %c (%s) at level %i (%s),', 64+1, FL(1).expr, ii, FL(1).level(ii).expr);
00472 fprintf (2,'\n factor %c (%s) at level %i (%s),', 64+2, FL(2).expr, jj, FL(2).level(jj).expr);
00473 fprintf (2,'\n factor %c (%s) at level %i (%s),', 64+3, FL(3).expr, kk, FL(3).level(kk).expr);
00474 sz = input('\nsample size is: ');
00475 fprintf (2,'\nProvide those %i input files:\n', sz);
00476 for (rr = 1:1:sz),
00477 FI = FI + 1;
00478 GP(1, FI) = {FL(1).level(ii).expr};
00479 GP(2, FI) = {FL(2).level(jj).expr};
00480 GP(3, FI) = {FL(3).level(kk).expr};
00481
00482 flg = 0;
00483 while flg == 0,
00484 fprintf (2,'No. %i file ', rr);
00485 file(FI).nm = input('is: ', 's');
00486 fid = fopen (file(FI).nm,'r');
00487 if (fid == -1), flg = 0; fprintf(2,'Error: File %s does not exist. Please try it again. \n', file(FI).nm);
00488 else flg = 1; fclose (fid); end
00489 if isempty(strfind(file(FI).nm, 'tlrc')) == 0
00490 format = 'tlrc';
00491 elseif isempty(strfind(file(FI).nm, 'orig')) == 0
00492 format = 'orig';
00493 else
00494 if isempty(strfind(file(FI).nm, '1D')) == 0 % if 1D file (surface data)
00495 format = '1D';
00496 else
00497 while (1); fprintf(2,'Error: format of file %s is incorrect!\n', file(FI).nm);
00498 fprintf(2,'Halted: Ctrl+c to exit'); pause; end
00499 end
00500 end
00501 end
00502
00503 end % for (rr = 1:1:sz)
00504 end % for (kk = 1:1:FL(3).N_level)
00505 end % for (jj = 1:1:FL(2).N_level)
00506 end % for (ii = 1:1:FL(1).N_level)
00507
00508 case 4,
00509 for (ii = 1:1:FL(1).N_level),
00510 for (jj = 1:1:FL(2).N_level),
00511 for (kk = 1:1:FL(3).N_level),
00512 for (ll = 1:1:FL(4).N_level),
00513 fprintf (2,'\nFor factor %c (%s) at level %i (%s),', 64+1, FL(1).expr, ii, FL(1).level(ii).expr);
00514 fprintf (2,'\n factor %c (%s) at level %i (%s),', 64+2, FL(2).expr, jj, FL(2).level(jj).expr);
00515 fprintf (2,'\n factor %c (%s) at level %i (%s),', 64+3, FL(3).expr, kk, FL(3).level(kk).expr);
00516 fprintf (2,'\n factor %c (%s) at level %i (%s),', 64+3, FL(3).expr, kk, FL(4).level(ll).expr);
00517 sz = input('\nsample size is: ');
00518 fprintf (2,'\nProvide those %i input files:\n', sz);
00519 for (rr = 1:1:sz),
00520 FI = FI + 1;
00521 GP(1, FI) = {FL(1).level(ii).expr};
00522 GP(2, FI) = {FL(2).level(jj).expr};
00523 GP(3, FI) = {FL(3).level(kk).expr};
00524 GP(4, FI) = {FL(4).level(ll).expr};
00525 flg = 0;
00526 while flg == 0,
00527 fprintf (2,'No. %i file ', rr);
00528 file(FI).nm = input('is: ', 's');
00529 fid = fopen (file(FI).nm,'r');
00530 if (fid == -1), flg = 0; fprintf(2,'Error: File %s does not exist. Please try it again. \n', file(FI).nm);
00531 else flg = 1; fclose (fid); end
00532 if isempty(strfind(file(FI).nm, 'tlrc')) == 0
00533 format = 'tlrc';
00534 elseif isempty(strfind(file(FI).nm, 'orig')) == 0
00535 format = 'orig';
00536 else
00537 if isempty(strfind(file(FI).nm, '1D')) == 0 % if 1D file (surface data)
00538 format = '1D';
00539 else
00540 while (1); fprintf(2,'Error: format of file %s is incorrect!\n', file(FI).nm);
00541 fprintf(2,'Halted: Ctrl+c to exit'); pause; end
00542 end
00543 end
00544 end
00545
00546 end % for (rr = 1:1:sz)
00547 end % for (ll = 1:1:FL(4).N_level)
00548 end % for (kk = 1:1:FL(3).N_level)
00549 end % for (jj = 1:1:FL(2).N_level)
00550 end % for (ii = 1:1:FL(1).N_level)
00551
00552 end % switch NF
00553 if (ntot == FI), fprintf (2,'\n%i input files have been read in. \n', FI);
00554 else fprintf(2,'Error: Total number of files do not match up. \n');
00555 while (1); fprintf(2,'Halted: Ctrl+c to exit'); pause; end
00556 end
00557 end % if (dsgn == 1)
00558
00559 if ((NF == 3 & dsgn == 3)),
00560 % acc = 0;
00561 FI = 0; % File index
00562 % if (unbalanced.type == 1),
00563 for (i = 1:1:FL(1).N_level),
00564 for (j = 1:1:FL(2).N_level),
00565 for (k = 1:1:FL(3).UL(i).N_level),
00566
00567 % FI = acc + (j-1)*FL(3).UL(i).N_level + k; % file index
00568 % Create a matrix for group indices
00569 % GP(1, FI) = {FL(1).level(i).expr};
00570 % GP(2, FI) = {FL(2).level(j).expr};
00571 % GP(3, FI) = {FL(3).UL(i).n(k).expr};
00572
00573 for (r = 1:1:FL(NF+1).N_level), % if there is any repeated observations
00574 FI = FI + 1;
00575 GP(1, FI) = {FL(1).level(i).expr};
00576 GP(2, FI) = {FL(2).level(j).expr};
00577 GP(3, FI) = {FL(3).UL(i).n(k).expr};
00578
00579 flg = 0;
00580 while flg == 0,
00581 fprintf (2,'\n(%i) factor combination:\n', FI);
00582 for (m=1:1:NF),
00583 fprintf('\tfactor %c (%s) at level %s \n', 64+m, FL(m).expr, char(GP(m, FI)));
00584 end
00585 fprintf('\tat repeat %i \n', r);
00586 file(FI).nm = input('is: ', 's');
00587 fid = fopen (file(FI).nm,'r');
00588 if (fid == -1),
00589 flg = 0; fprintf(2,'Error: File %s does not exist. Please try it again. \n', file(FI).nm);
00590 else flg = 1; fclose (fid); end
00591 if isempty(strfind(file(FI).nm, 'tlrc')) == 0
00592 format = 'tlrc';
00593 elseif isempty(strfind(file(FI).nm, 'orig')) == 0
00594 format = 'orig';
00595 else
00596 if isempty(strfind(file(FI).nm, '1D')) == 0 % if 1D file (surface data)
00597 format = '1D';
00598 else
00599 while (1); fprintf(2,'Error: format of file %s is incorrect!\n', file(FI).nm);
00600 fprintf(2,'Halted: Ctrl+c to exit'); pause; end
00601 end
00602 end
00603 end
00604 end % for (r = 1:1:FL(NF+1).N_level),
00605
00606 end % k
00607 end % j
00608 % acc = FI;
00609 end % i
00610
00611 % end % if (unbalanced.type == 1),
00612 end % if ((NF == 3 & dsgn == 3))
00613
00614 if ((NF == 4 & dsgn == 3)),
00615 % acc = 0;
00616 % if (unbalanced.type == 1),
00617 FI = 0; % File index
00618 for (i = 1:1:FL(1).N_level),
00619 for (j = 1:1:FL(2).N_level),
00620 for (k = 1:1:FL(3).N_level),
00621 for (l = 1:1:FL(4).UL(i).N_level),
00622
00623 % FI = acc + (j-1)*FL(3).N_level*FL(4).UL(i).N_level+(k-1)*FL(4).UL(i).N_level + l; % file index
00624 % Create a matrix for group indices
00625 % GP(1, FI) = {FL(1).level(i).expr};
00626 % GP(2, FI) = {FL(2).level(j).expr};
00627 % GP(3, FI) = {FL(3).level(k).expr};
00628 % GP(4, FI) = {FL(4).UL(i).n(l).expr};
00629
00630 for (r = 1:1:FL(NF+1).N_level),
00631 FI = FI +1;
00632 GP(1, FI) = {FL(1).level(i).expr};
00633 GP(2, FI) = {FL(2).level(j).expr};
00634 GP(3, FI) = {FL(3).level(k).expr};
00635 GP(4, FI) = {FL(4).UL(i).n(l).expr};
00636
00637 flg = 0;
00638 while flg == 0,
00639 fprintf (2,'\n(%i) factor combination:\n', FI);
00640 for (m=1:1:NF),
00641 fprintf('\tfactor %c (%s) at level %s \n', 64+m, FL(m).expr, char(GP(m, FI)));
00642 end
00643 fprintf('\tat repeat %i \n', r);
00644 file(FI).nm = input('is: ', 's');
00645 fid = fopen (file(FI).nm,'r');
00646 if (fid == -1),
00647 flg = 0; fprintf(2,'Error: File %s does not exist. Please try it again. \n', file(FI).nm);
00648 else flg = 1; fclose (fid); end
00649
00650 if isempty(strfind(file(FI).nm, 'tlrc')) == 0
00651 format = 'tlrc';
00652 elseif isempty(strfind(file(FI).nm, 'orig')) == 0
00653 format = 'orig';
00654 else
00655 if isempty(strfind(file(FI).nm, '1D')) == 0 % if 1D file (surface data)
00656 format = '1D';
00657 else
00658 while (1); fprintf(2,'Error: format of file %s is incorrect!\n', file(FI).nm);
00659 fprintf(2,'Halted: Ctrl+c to exit'); pause; end
00660 end
00661 end
00662 end
00663 end % for (r = 1:1:FL(NF+1).N_level),
00664
00665 end
00666 end
00667 end
00668 % acc = FI;
00669 end
00670
00671 % end % if (unbalanced.type == 1),
00672 end % if ((NF == 4 & dsgn == 3))
00673 %end % (unbalanced.yes == 1),
00674
00675
00676 else % Balanced designs
00677
00678 fprintf('\nPlease provide input files.');
00679 for (i=1:1:ntot),
00680 %Ziad Saad modified Matlab function ind2sub for the purpose here
00681 %Converting the index into multiple subscripts
00682 %I want vary the levels starting the last factor in stead of the first.
00683 [err, file(i).v] = gind2sub (fliplr(sz), i); %flip before subscripting. Doing flip because
00684 %I want vary the levels starting the last factor in stead of the first.
00685 scpt = fliplr(file(i).v); %flip back to restore the original order
00686
00687 flg = 0;
00688 while flg == 0,
00689 fprintf (2,'\n(%i) factor combination:\n', i);
00690 for (k=1:1:NF),
00691 fprintf('\tfactor %c (%s) at level %i (%s) \n', 64+k, FL(k).expr, scpt(k), FL(k).level(scpt(k)).expr);
00692 GP(k, i) = {FL(k).level(scpt(k)).expr};
00693 end
00694 fprintf('\tat repeat %i \n', scpt(NF+1));
00695 file(i).nm = input('is: ', 's');
00696 fid = fopen (file(i).nm,'r');
00697 if (fid == -1),
00698 flg = 0; fprintf(2,'Error: File %s does not exist. Please try it again. \n', file(i).nm);
00699 else flg = 1; fclose (fid); end
00700
00701 if isempty(strfind(file(i).nm, 'tlrc')) == 0
00702 format = 'tlrc';
00703 elseif isempty(strfind(file(i).nm, 'orig')) == 0
00704 format = 'orig';
00705 else
00706 if isempty(strfind(file(i).nm, '1D')) == 0 % if 1D file (surface data)
00707 format = '1D';
00708 else
00709 while (1); fprintf(2,'Error: format of file %s is incorrect!\n', file(i).nm);
00710 fprintf(2,'Halted: Ctrl+c to exit'); pause; end
00711 end
00712 end
00713 end
00714 end
00715 end % if (unbalanced.yes == 1),
00716 end % if (file_format == 1),
00717
00718 if (file_format == 0), % unused now!
00719
00720 if (unbalanced.yes == 1), % Try 4-way design 3 first
00721 if ((NF == 4 & dsgn == 3)),
00722 acc = 0;
00723 if (unbalanced.type == 1),
00724 for (i = 1:1:FL(1).N_level),
00725 for (j = 1:1:FL(2).N_level),
00726 for (k = 1:1:FL(3).N_level),
00727 for (l = 1:1:FL(4).UL(i).N_level),
00728
00729 FI = acc + (j-1)*FL(3).N_level*FL(4).UL(i).N_level+(k-1)*FL(4).UL(i).N_level + l; % file index
00730
00731 GP(1, FI) = {FL(1).level(i).expr};
00732 GP(2, FI) = {FL(2).level(j).expr};
00733 GP(3, FI) = {FL(3).level(k).expr};
00734 GP(4, FI) = {FL(4).UL(i).n(l).expr};
00735
00736 end
00737 end
00738 end
00739 acc = FI;
00740 end
00741
00742 for (m = 1:1:FL(1).N_level),
00743 counter(m) = 0;
00744 end %Just initiation
00745
00746 for (i=1:1:file_num),
00747
00748 flg = 0;
00749 while flg == 0,
00750 fprintf('Input file for subject #%i ', i);
00751 file(i).nm = input('is: ', 's');
00752 fid = fopen (file(i).nm,'r');
00753 fprintf('\tThe corresponding level of factor 1 for subject %i ', i);
00754 file(i).F1L = input('is: '); % level for factor 1
00755 if (isnumeric(file(i).F1L) == 0 | isempty(file(i).F1L)),
00756 flg = 0; fprintf(2,'Error: the input is not a number. Please try it again.\n');
00757 end
00758
00759 for (m = 1:1:FL(1).N_level),
00760 if (file(i).F1L == m),
00761 counter(m) = counter(m) + 1;
00762 file(i).cntr(m) = counter(m); % Any better way to implement this without the temporary array of counter(j)?
00763 end % set a counter of 4th factor level marker for later use
00764 end
00765
00766 if (fid == -1),
00767 flg = 0; fprintf(2,'Error: File %s does not exist. Please try it again. \n', file(i).nm);
00768 else flg = 1; fclose (fid); end
00769 end
00770 end % for (i=1:1:file_num),
00771
00772 fprintf('\nIt is assumed that all files have the same subbrik order.');
00773 fprintf('\nFor each subbrik of an input file, specify the factor level of 2nd and 3rd factors. \n');
00774
00775 for (j = 1:1:file_SB),
00776 fprintf('\n\tFor subbrik %i: \n', j-1);
00777 for (k = 1:1:(NF-2)), % Here we only keep record of 2nd and 3rd factor levels
00778 fprintf('\t\tlevel of factor %c (%s) ', 65+k, FL(k+1).expr);
00779 SB(j).lv(k) = input('is: ');
00780 end %for (j = 1:1:file_SB),
00781 end % for (j = 1:1:file_SB),
00782
00783 end % if (unbalanced.type == 1),
00784 end % if ((NF == 4 & dsgn == 3))
00785 %end % (unbalanced.yes == 1),
00786 else
00787
00788
00789
00790 for (i=1:1:ntot),
00791 [err, file(i).v] = gind2sub (fliplr(sz), i); %flip before subscripting. Doing flip because
00792 %I want vary the levels starting the last factor in stead of the first.
00793 scpt = fliplr(file(i).v); %flip back to restore the original order
00794 for (k=1:1:NF),
00795 GP(k, i) = {FL(k).level(scpt(k)).expr};
00796 end
00797 end
00798
00799 for (j = 1:1:FL(1).N_level),
00800 counter(j) = 0;
00801 end %Just initiation
00802 for (i = 1:1:file_num),
00803 flg = 0;
00804 while flg == 0,
00805 fprintf('Input file for subject #%i ', i);
00806 file(i).nm = input('is: ', 's');
00807 fid = fopen (file(i).nm,'r');
00808 fprintf('\tThe corresponding level of factor 1 for subject %i ', i);
00809 file(i).F1L = input('is: '); % level for factor 1
00810 if (isnumeric(file(i).F1L) == 0 | isempty(file(i).F1L)),
00811 flg = 0; fprintf(2,'Error: the input is not a number. Please try it again.\n');
00812 end
00813
00814 for (j = 1:1:FL(1).N_level),
00815 if (file(i).F1L == j),
00816 counter(j) = counter(j) + 1;
00817 file(i).cntr(j) = counter(j); % Any better way to implement this without the temporary array of counter(j)?
00818 end % set a counter of 4th factor level marker for later use
00819 end
00820
00821 if (fid == -1),
00822 flg = 0; fprintf(2,'Error: File %s does not exist. Please try it again. \n', file(i).nm);
00823 else flg = 1; fclose (fid); end
00824 end
00825 end
00826 fprintf('\nIt is assumed that all files have the same subbrik order.');
00827 fprintf('\nFor each subbrik of an input file, specify the factor level of 2nd and 3rd factors. \n');
00828
00829 for (j = 1:1:file_SB),
00830 fprintf('\n\tFor subbrik %i: \n', j-1);
00831 for (k = 1:1:(NF-2)), %Here we only keep record of 2nd and 3rd factor levels
00832 fprintf('\t\tlevel of factor %c (%s) ', 65+k, FL(k+1).expr);
00833 SB(j).lv(k) = input('is: ');
00834 end
00835 end
00836 end % if (unbalanced.yes == 1)
00837 end % if (file_format == 0),
00838
00839
00840
00841 %Obtain output file name
00842 flg = 0;
00843 while flg == 0,
00844 OutFN = input('\nOutput file name (in bucket format): ', 's');
00845 OutFull = sprintf('%s+%s.HEAD', OutFN, format);
00846 fid2 = fopen(OutFull,'r');
00847 if (fid2 ~= -1),
00848 flg = 0; fprintf(2,'Error: File %s exists already. Please give another name. \n', OutFN);
00849 else flg = 1; end
00850 end
00851
00852 % Gather contrast information: Everything stored in structure Contr
00853
00854 flg = 0;
00855 while flg == 0,
00856 Contr.do = input('\nAny contrast test (1 - Yes, 0 - No)? ');
00857 if (Contr.do ~= 0 & Contr.do ~= 1),
00858 flg = 0; fprintf(2,'Error: invalid answer. Try it again. \n', OutFN);
00859 else flg = 1; end
00860 end
00861
00862 if (Contr.do == 0),
00863 Contr.ord1.tot = 0; % assign these for output later
00864 Contr.ord2.tot = 0;
00865 Contr.ord3.tot = 0;
00866 Contr.ord4.tot = 0;
00867 else
00868
00869 fprintf('\nNow coding contrasts:\n');
00870 fprintf('\n\t1. Each contrast should contain at least 2 terms;');
00871 fprintf('\n\t2. Each term in a contrast should have %i character(s), \n', NF);
00872 for (i = 1:1:NF), fprintf('\tNo. %i character corresponds to the level of factor %c (%s),\n', i, 'A'+i-1, FL(i).expr); end
00873 fprintf('\n\tUse 0 if a factor is collapsed.\n');
00874 fprintf('\n\tIf a factor level is smaller than 9, use its ordinal number;');
00875 fprintf('\n\tIf a factor level is bigger than 9, use a, b, c, ... (no capitals) for 10, 11, 12, ... \n');
00876 fprintf('\n\t3. All weights/coeficients in a contrast have to add up to zero.\n');
00877
00878 if (NF == 1),
00879 % 1st order contrasts
00880 flg = 0;
00881
00882 while flg == 0,
00883 fprintf('\n1st order contrasts have %i factor(s) collapsed.\n', NF-1);
00884 Contr.ord1.tot = input('\nHow many 1st-order contrasts? (0 if none) ');
00885 if (isnumeric(Contr.ord1.tot) == 0 | Contr.ord1.tot < 0),
00886 flg = 0; fprintf(2,'Error: inapproriate input. Please try again.\n');
00887 else flg = 1;
00888 end
00889 end
00890 % fprintf('\nUse factor level to code each term in a contrast. For example, 1 means the factor is at first level.\n');
00891 for (i = 1:1:Contr.ord1.tot),
00892 fprintf('\nLabel for 1st order contrast No. %i ', i);
00893 Contr.ord1.label(i).nm = input('is: ', 's');
00894 flg = 0;
00895 while flg == 0,
00896 Contr.ord1.cnt(i).NT = input('How many terms are involved in this contrast? '); % NT = number of terms involved in this contrast
00897 if (isnumeric(Contr.ord1.cnt(i).NT) == 0 | Contr.ord1.cnt(i).NT < 2),
00898 flg = 0; fprintf(2,'\nError: inapproriate input. Please try again.\n\n');
00899 else flg = 1; end
00900 end
00901
00902 flg0 = 0;
00903 while flg0 == 0,
00904 for (j = 1:1:Contr.ord1.cnt(i).NT),
00905 flg = 0;
00906 while flg == 0,
00907 fprintf('Factor index for No. %i term ', j);
00908 Contr.ord1.cnt(i).code(j).str = input('is (e.g., 2): ', 's');
00909 if (length(Contr.ord1.cnt(i).code(j).str) ~= NF),
00910 flg = 0; fprintf(2,'Error: invalid input. Try it again. \n', OutFN);
00911 else flg = 1; end
00912 end
00913 Contr.ord1.cnt(i).coef(j) = input('Corresponding coefficient (e.g., 1 or -1): ');
00914 end
00915 if (sum(Contr.ord1.cnt(i).coef) ~= 0),
00916 flg0 = 0; fprintf(2,'\nError: All the coeficients of a contrast have to sum up to 0! Try again...\n\n');
00917 else flg0 = 1; end
00918 end
00919 end
00920 Contr.ord2.tot = 0;
00921 Contr.ord3.tot = 0;
00922 Contr.ord4.tot = 0;
00923 fprintf(1,'Done with 1st order contrast information.\n');
00924 end
00925
00926 if (NF == 2),
00927 % 1st order contrasts ONLY at this point
00928 flg = 0;
00929 while flg == 0,
00930 fprintf('\n1st order contrasts have %i factor(s) collapsed.\n', NF-1);
00931 Contr.ord1.tot = input('\nHow many 1st-order contrasts? (0 if none) ');
00932 if (isnumeric(Contr.ord1.tot) == 0 | Contr.ord1.tot < 0),
00933 flg = 0; fprintf(2,'\nError: inapproriate input. Please try again.\n');
00934 else flg = 1;
00935 end
00936 end
00937 % fprintf('\nUse factor level to code each term in a contrast. For example, 1 means the factor is at first level.\n');
00938 for (i = 1:1:Contr.ord1.tot),
00939 fprintf('\nLabel for 1st order contrast No. %i ', i);
00940 Contr.ord1.label(i).nm = input('is: ', 's');
00941 flg = 0;
00942 while flg == 0,
00943 Contr.ord1.cnt(i).NT = input('How many terms are involved in this contrast? '); % NT = number of terms involved in this contrast
00944 if (isnumeric(Contr.ord1.cnt(i).NT) == 0 | Contr.ord1.cnt(i).NT < 2),
00945 flg = 0; fprintf(2,'\nError: inapproriate input. Please try again.\n\n');
00946 else flg = 1; end
00947 end
00948
00949 flg0 = 0;
00950 while flg0 == 0,
00951 for (j = 1:1:Contr.ord1.cnt(i).NT),
00952 flg = 0;
00953 while flg == 0,
00954 fprintf('Factor index for No. %i term ', j);
00955 Contr.ord1.cnt(i).code(j).str = input('is (e.g., 10): ', 's');
00956 if (length(Contr.ord1.cnt(i).code(j).str) ~= NF),
00957 flg = 0; fprintf(2,'\nError: invalid input. Try it again. \n', OutFN);
00958 else flg = 1; end
00959 end
00960 Contr.ord1.cnt(i).coef(j) = input('Corresponding coefficient (e.g., 1 or -1): ');
00961 end
00962 if (sum(Contr.ord1.cnt(i).coef) ~= 0),
00963 flg0 = 0; fprintf(2,'\nError: All the coeficients of a contrast have to sum up to 0! Try again...\n\n');
00964 else flg0 = 1; end
00965 end
00966 end
00967 fprintf(1,'Done with 1st order contrast information.\n');
00968 % 2nd order contrasts
00969 flg = 0;
00970 while flg == 0,
00971 fprintf('\n2nd order contrasts have %i factor(s) collapsed.\n', NF-2);
00972 fprintf('\nNotice: Contrasts for random factor are NOT feasible.\n');
00973 Contr.ord2.tot = input('\nHow many 2nd-order contrasts? (0 if none) ');
00974 if (isnumeric(Contr.ord2.tot) == 0 | Contr.ord2.tot < 0),
00975 flg = 0; fprintf(2,'\nError: inapproriate input. Please try again.\n');
00976 else flg = 1;
00977 end
00978 end
00979 % fprintf('\nUse factor level to code each term in a contrast. For example, 12 means ');
00980 % fprintf('\nthe 1st and 2nd factors are at their first and second level respectively.\n');
00981 for (i = 1:1:Contr.ord2.tot),
00982 fprintf('\nLabel for 2nd order contrast No. %i: ', i);
00983 Contr.ord2.label(i).nm = input('is: ', 's');
00984 flg = 0;
00985 while flg == 0,
00986 Contr.ord2.cnt(i).NT = input('How many terms are involved in this contrast? '); % NT = number of terms involved in this contrast
00987 if (isnumeric(Contr.ord2.cnt(i).NT) == 0 | Contr.ord2.cnt(i).NT < 2),
00988 flg = 0; fprintf(2,'\nError: inapproriate input. Please try again.\n\n');
00989 else flg = 1; end
00990 end
00991
00992 flg0 = 0;
00993 while flg0 == 0,
00994 for (j = 1:1:Contr.ord2.cnt(i).NT),
00995 flg = 0;
00996 while flg == 0,
00997 fprintf('Factor index for No. %i term ', j);
00998 Contr.ord2.cnt(i).code(j).str = input('is (e.g., 12): ', 's');
00999 if (length(Contr.ord2.cnt(i).code(j).str) ~= NF),
01000 flg = 0; fprintf(2,'\nError: invalid input. Try it again. \n', OutFN);
01001 else flg = 1; end
01002 end
01003 Contr.ord2.cnt(i).coef(j) = input('Corresponding coefficient (e.g., 1 or -1): ');
01004 end
01005 if (sum(Contr.ord2.cnt(i).coef) ~= 0),
01006 flg0 = 0; fprintf(2,'\nError: All the coeficients of a contrast have to sum up to 0! Try again...\n\n');
01007 else flg0 = 1; end
01008 end
01009 end
01010
01011 Contr.ord3.tot = 0;
01012 Contr.ord4.tot = 0;
01013 fprintf(1,'Done with 2nd order contrast information.\n');
01014 end
01015
01016 if (NF == 3 | NF == 4),
01017
01018 % 1st order contrasts
01019 flg = 0;
01020 while flg == 0,
01021 fprintf('\n1st order contrasts have %i factor(s) collapsed.\n', NF-1);
01022 %fprintf('\t1. Simple mean for a factor level, such as [1 0 0];\n');
01023 %fprintf('\t2. Difference between two factor levels, such as [1 -1 0];\n');
01024 %fprintf('\t3. Linear combination of factor levels, such as [0.5 0.5 -1];\n');
01025 fprintf('\nNotice: Contrasts for random factor are NOT feasible.\n');
01026 Contr.ord1.tot = input('\nHow many 1st-order contrasts? (0 if none) ');
01027 if (isnumeric(Contr.ord1.tot) == 0 | Contr.ord1.tot < 0),
01028 flg = 0; fprintf(2,'\nError: inapproriate input. Please try again.\n');
01029 else flg = 1;
01030 end
01031 end
01032 % fprintf('\nUse factor level to code each term in a contrast. For example, 0100 means the first, third ');
01033 % fprintf('\nand fourth factors are collapsed while 2nd factor is at first level.\n');
01034 for (i = 1:1:Contr.ord1.tot),
01035 fprintf('\nLabel for 1st order contrast No. %i ', i);
01036 Contr.ord1.label(i).nm = input('is: ', 's');
01037 flg = 0;
01038 while flg == 0,
01039 Contr.ord1.cnt(i).NT = input('How many terms are involved in this contrast? '); % NT = number of terms involved in this contrast
01040 if (isnumeric(Contr.ord1.cnt(i).NT) == 0 | Contr.ord1.cnt(i).NT < 2),
01041 flg = 0; fprintf(2,'\nError: inapproriate input. Please try again.\n\n');
01042 else flg = 1; end
01043 end
01044
01045 flg0 = 0;
01046 while flg0 == 0,
01047 for (j = 1:1:Contr.ord1.cnt(i).NT),
01048 flg = 0;
01049 while flg == 0,
01050 fprintf('Factor index for No. %i term ', j);
01051 Contr.ord1.cnt(i).code(j).str = input('is (e.g., 0020): ', 's');
01052 if (length(Contr.ord1.cnt(i).code(j).str) ~= NF),
01053 flg = 0; fprintf(2,'\nError: invalid input. Try it again. \n', OutFN);
01054 else flg = 1; end
01055 end
01056 Contr.ord1.cnt(i).coef(j) = input('Corresponding coefficient (e.g., 1 or -1): ');
01057 end
01058 if (sum(Contr.ord1.cnt(i).coef) ~= 0),
01059 flg0 = 0; fprintf(2,'\nError: All the coeficients of a contrast have to sum up to 0! Try again...\n\n');
01060 else flg0 = 1; end
01061 end
01062 end
01063 fprintf(1,'Done with 1st order contrast information.\n');
01064
01065 % 2nd order contrasts
01066 flg = 0;
01067 while flg == 0,
01068 fprintf('\n2nd order contrasts have %i factor(s) collapsed.\n', NF-2);
01069 fprintf('\nNotice: Contrasts for random factor are NOT feasible.\n');
01070 Contr.ord2.tot = input('\nHow many 2nd-order contrasts? (0 if none) ');
01071 if (isnumeric(Contr.ord2.tot) == 0 | Contr.ord2.tot < 0),
01072 flg = 0; fprintf(2,'\nError: inapproriate input. Please try again.\n');
01073 else flg = 1;
01074 end
01075 end
01076 % fprintf('\nUse factor level to code each term in a contrast. For example, 0120 means both the first ');
01077 % fprintf('\nand fourth factors are collapsed while 2nd and 3rd factors are at first and second level respectively.\n');
01078 for (i = 1:1:Contr.ord2.tot),
01079 fprintf('\nLabel for 2nd order contrast No. %i ', i);
01080 Contr.ord2.label(i).nm = input('is: ', 's');
01081 flg = 0;
01082 while flg == 0,
01083 Contr.ord2.cnt(i).NT = input('How many terms are involved in this contrast? '); % NT = number of terms involved in this contrast
01084 if (isnumeric(Contr.ord2.cnt(i).NT) == 0 | Contr.ord2.cnt(i).NT < 2),
01085 flg = 0; fprintf(2,'\nError: inapproriate input. Please try again.\n\n');
01086 else flg = 1; end
01087 end
01088
01089 flg0 = 0;
01090 while flg0 == 0,
01091 for (j = 1:1:Contr.ord2.cnt(i).NT),
01092 flg = 0;
01093 while flg == 0,
01094 fprintf('Factor index for No. %i term ', j);
01095 Contr.ord2.cnt(i).code(j).str = input('is (e.g., 0100): ', 's');
01096 if (length(Contr.ord2.cnt(i).code(j).str) ~= NF),
01097 flg = 0; fprintf(2,'Error: invalid input. Try it again. \n', OutFN);
01098 else flg = 1; end
01099 end
01100 Contr.ord2.cnt(i).coef(j) = input('Corresponding coefficient (e.g., 1 or -1): ');
01101 end
01102 if (sum(Contr.ord2.cnt(i).coef) ~= 0),
01103 flg0 = 0; fprintf(2,'\nError: All the coeficients of a contrast have to sum up to 0! Try again...\n\n');
01104 else flg0 = 1; end
01105 end %flg0 = 0;
01106 end
01107
01108 fprintf(1,'Done with 2nd order contrast information.\n');
01109
01110
01111 % 3rd-order contrasts
01112 flg = 0;
01113 while flg == 0,
01114 fprintf('\n3rd order contrasts have %i factor(s) collapsed.\n', NF-3);
01115 fprintf('\nNotice: Contrasts for random factor are NOT feasible.\n');
01116 Contr.ord3.tot = input('\nHow many 3rd-order contrasts? (0 if none) ');
01117 if (isnumeric(Contr.ord3.tot) == 0 | Contr.ord3.tot < 0),
01118 flg = 0; fprintf(2,'\nError: inapproriate input. Please try again.\n');
01119 else flg = 1;
01120 end
01121 end
01122 % fprintf('\nUse factor level to code each term in a contrast. For example, 1230 means the fourth');
01123 % fprintf('\nfactor is collapsed while all the other 3 factors are at first, second and third level respectively.\n');
01124 for (i = 1:1:Contr.ord3.tot),
01125 fprintf('\nLabel for 3rd order contrast No. %i ', i);
01126 Contr.ord3.label(i).nm = input('is: ', 's');
01127 flg = 0;
01128 while flg == 0,
01129 Contr.ord3.cnt(i).NT = input('How many terms are involved in this contrast? '); % NT = number of terms involved in this contrast
01130 if (isnumeric(Contr.ord3.cnt(i).NT) == 0 | Contr.ord3.cnt(i).NT < 2),
01131 flg = 0; fprintf(2,'\nError: inapproriate input. Please try again.\n\n');
01132 else flg = 1; end
01133 end
01134
01135 flg0 = 0;
01136 while flg0 == 0,
01137 for (j = 1:1:Contr.ord3.cnt(i).NT),
01138 flg = 0;
01139 while flg == 0,
01140 fprintf('Factor index for No. %i term ', j);
01141 Contr.ord3.cnt(i).code(j).str = input('is (e.g., 1230): ', 's');
01142 if (length(Contr.ord3.cnt(i).code(j).str) ~= NF),
01143 flg = 0; fprintf(2,'\nError: invalid input. Try it again. \n', OutFN);
01144 else flg = 1; end
01145 end
01146 Contr.ord3.cnt(i).coef(j) = input('Corresponding coefficient (e.g., 1 or -1): ');
01147 end
01148 if (sum(Contr.ord3.cnt(i).coef) ~= 0),
01149 flg0 = 0; fprintf(2,'\nError: All the coeficients of a contrast have to sum up to 0! Try again...\n\n');
01150 else flg0 = 1; end
01151 end %flg0 = 0;
01152 end
01153 Contr.ord4.tot = 0;
01154 fprintf(1,'Done with 3rd order contrast information.\n');
01155
01156 end % if (NF == 4)
01157
01158 if (NF == 5),
01159
01160 % 1st order contrasts
01161 flg = 0;
01162 while flg == 0,
01163 fprintf('\n1st order contrasts have %i factor(s) collapsed.\n', NF-1);
01164 %fprintf('\t1. Simple mean for a factor level, such as [1 0 0];\n');
01165 %fprintf('\t2. Difference between two factor levels, such as [1 -1 0];\n');
01166 %fprintf('\t3. Linear combination of factor levels, such as [0.5 0.5 -1];\n');
01167 fprintf('\nNotice: Contrasts for random factor are NOT feasible.\n');
01168 Contr.ord1.tot = input('\nHow many 1st-order contrasts? (0 if none) ');
01169 if (isnumeric(Contr.ord1.tot) == 0 | Contr.ord1.tot < 0),
01170 flg = 0; fprintf(2,'\nError: inapproriate input. Please try again.\n');
01171 else flg = 1;
01172 end
01173 end
01174 % fprintf('\nUse factor level to code each term in a contrast. For example, 0100 means the first, third ');
01175 % fprintf('\nand fourth factors are collapsed while 2nd factor is at first level.\n');
01176 for (i = 1:1:Contr.ord1.tot),
01177 fprintf('\nLabel for 1st order contrast No. %i ', i);
01178 Contr.ord1.label(i).nm = input('is: ', 's');
01179 flg = 0;
01180 while flg == 0,
01181 Contr.ord1.cnt(i).NT = input('How many terms are involved in this contrast? '); % NT = number of terms involved in this contrast
01182 if (isnumeric(Contr.ord1.cnt(i).NT) == 0 | Contr.ord1.cnt(i).NT < 2),
01183 flg = 0; fprintf(2,'\nError: inapproriate input. Please try again.\n\n');
01184 else flg = 1; end
01185 end
01186
01187 flg0 = 0;
01188 while flg0 == 0,
01189 for (j = 1:1:Contr.ord1.cnt(i).NT),
01190 flg = 0;
01191 while flg == 0,
01192 fprintf('Factor index for No. %i term ', j);
01193 Contr.ord1.cnt(i).code(j).str = input('is (e.g., 00200): ', 's');
01194 if (length(Contr.ord1.cnt(i).code(j).str) ~= NF),
01195 flg = 0; fprintf(2,'Error: invalid input. Try it again. \n', OutFN);
01196 else flg = 1; end
01197 end
01198 Contr.ord1.cnt(i).coef(j) = input('Corresponding coefficient (e.g., 1 or -1): ');
01199 end % for (j = 1:1:Contr.ord1.cnt(i).NT)
01200 if (sum(Contr.ord1.cnt(i).coef) ~= 0),
01201 flg0 = 0; fprintf(2,'\nError: All the coeficients of a contrast have to sum up to 0! Try again...\n\n');
01202 else flg0 = 1; end
01203 end %flg0 = 0;
01204 end
01205 fprintf(1,'Done with 1st order contrast information.\n');
01206
01207 % 2nd order contrasts
01208 flg = 0;
01209 while flg == 0,
01210 fprintf('\n2nd order contrasts have %i factor(s) collapsed.\n', NF-2);
01211 fprintf('\nNotice: Contrasts for random factor are NOT feasible.\n');
01212 Contr.ord2.tot = input('\nHow many 2nd-order contrasts? (0 if none) ');
01213 if (isnumeric(Contr.ord2.tot) == 0 | Contr.ord2.tot < 0),
01214 flg = 0; fprintf(2,'Error: inapproriate input. Please try again.\n');
01215 else flg = 1;
01216 end
01217 end % while flg == 0
01218 % fprintf('\nUse factor level to code each term in a contrast. For example, 0120 means both the first ');
01219 % fprintf('\nand fourth factors are collapsed while 2nd and 3rd factors are at first and second level respectively.\n');
01220 for (i = 1:1:Contr.ord2.tot),
01221 fprintf('\nLabel for 2nd order contrast No. %i ', i);
01222 Contr.ord2.label(i).nm = input('is: ', 's');
01223 flg = 0;
01224 while flg == 0,
01225 Contr.ord2.cnt(i).NT = input('How many terms are involved in this contrast? '); % NT = number of terms involved in this contrast
01226 if (isnumeric(Contr.ord2.cnt(i).NT) == 0 | Contr.ord2.cnt(i).NT < 2),
01227 flg = 0; fprintf(2,'\nError: inapproriate input. Please try again.\n\n');
01228 else flg = 1; end
01229 end
01230
01231 flg0 = 0;
01232 while flg0 == 0,
01233 for (j = 1:1:Contr.ord2.cnt(i).NT),
01234 flg = 0;
01235 while flg == 0,
01236 fprintf('Factor index for No. %i term ', j);
01237 Contr.ord2.cnt(i).code(j).str = input('is (e.g., 01200): ', 's');
01238 if (length(Contr.ord2.cnt(i).code(j).str) ~= NF),
01239 flg = 0; fprintf(2,'\nError: invalid input. Try it again. \n', OutFN);
01240 else flg = 1; end
01241 end % while flg == 0
01242 Contr.ord2.cnt(i).coef(j) = input('Corresponding coefficient (e.g., 1 or -1): ');
01243 end % for (j = 1:1:Contr.ord2.cnt(i).NT)
01244 if (sum(Contr.ord2.cnt(i).coef) ~= 0),
01245 flg0 = 0; fprintf(2,'\nError: All the coeficients of a contrast have to sum up to 0! Try again...\n\n');
01246 else flg0 = 1; end
01247 end %flg0 = 0;
01248 end % for (i = 1:1:Contr.ord2.tot)
01249 fprintf(1,'Done with 2nd order contrast information.\n');
01250
01251
01252 % 3rd-order contrasts
01253 flg = 0;
01254 while flg == 0,
01255 fprintf('\n3rd order contrasts have %i factor(s) collapsed.\n', NF-3);
01256 fprintf('\nNotice: Contrasts for random factor are NOT feasible.\n');
01257 Contr.ord3.tot = input('\nHow many 3rd-order contrasts? (0 if none) ');
01258 if (isnumeric(Contr.ord3.tot) == 0 | Contr.ord3.tot < 0),
01259 flg = 0; fprintf(2,'\nError: inapproriate input. Please try again.\n');
01260 else flg = 1;
01261 end
01262 end
01263 % fprintf('\nUse factor level to code each term in a contrast. For example, 1230 means the fourth');
01264 % fprintf('\nfactor is collapsed while all the other 3 factors are at first, second and third level respectively.\n');
01265 for (i = 1:1:Contr.ord3.tot),
01266 fprintf('\nLabel for 3rd order contrast No. %i ', i);
01267 Contr.ord3.label(i).nm = input('is: ', 's');
01268 flg = 0;
01269 while flg == 0,
01270 Contr.ord3.cnt(i).NT = input('How many terms are involved in this contrast? '); % NT = number of terms involved in this contrast
01271 if (isnumeric(Contr.ord3.cnt(i).NT) == 0 | Contr.ord3.cnt(i).NT < 2),
01272 flg = 0; fprintf(2,'\nError: inapproriate input. Please try again.\n\n');
01273 else flg = 1; end
01274 end
01275
01276 flg0 = 0;
01277 while flg0 == 0,
01278 for (j = 1:1:Contr.ord3.cnt(i).NT),
01279 flg = 0;
01280 while flg == 0,
01281 fprintf('Factor index for No. %i term ', j);
01282 Contr.ord3.cnt(i).code(j).str = input('is (e.g., 1230): ', 's');
01283 if (length(Contr.ord3.cnt(i).code(j).str) ~= NF),
01284 flg = 0; fprintf(2,'\nError: invalid input. Try it again. \n', OutFN);
01285 else flg = 1; end
01286 end
01287 Contr.ord3.cnt(i).coef(j) = input('Corresponding coefficient (e.g., 1 or -1): ');
01288 end
01289 if (sum(Contr.ord3.cnt(i).coef) ~= 0),
01290 flg0 = 0; fprintf(2,'\nError: All the coeficients of a contrast have to sum up to 0! Try again...\n\n');
01291 else flg0 = 1; end
01292 end %flg0 = 0;
01293 end
01294 fprintf(1,'Done with 3rd order contrast information.\n');
01295
01296 % 4th-order contrasts
01297 flg = 0;
01298 while flg == 0,
01299 fprintf('\n4th order contrasts have %i factor(s) collapsed.\n', NF-4);
01300 fprintf('\nNotice: Contrasts for random factor are NOT feasible.\n');
01301 Contr.ord4.tot = input('\nHow many 4th-order contrasts? (0 if none) ');
01302 if (isnumeric(Contr.ord4.tot) == 0 | Contr.ord4.tot < 0),
01303 flg = 0; fprintf(2,'\nError: inapproriate input. Please try again.\n');
01304 else flg = 1;
01305 end
01306 end
01307 % fprintf('\nUse factor level to code each term in a contrast. For example, 1230 means the fourth');
01308 % fprintf('\nfactor is collapsed while all the other 3 factors are at first, second and third level respectively.\n');
01309 for (i = 1:1:Contr.ord4.tot),
01310 fprintf('\nLabel for 4th order contrast No. %i ', i);
01311 Contr.ord4.label(i).nm = input('is: ', 's');
01312 flg = 0;
01313 while flg == 0,
01314 Contr.ord4.cnt(i).NT = input('How many terms are involved in this contrast? '); % NT = number of terms involved in this contrast
01315 if (isnumeric(Contr.ord4.cnt(i).NT) == 0 | Contr.ord4.cnt(i).NT < 2),
01316 flg = 0; fprintf(2,'\nError: inapproriate input. Please try again.\n\n');
01317 else flg = 1; end
01318 end
01319
01320 flg0 = 0;
01321 while flg0 == 0,
01322 for (j = 1:1:Contr.ord4.cnt(i).NT),
01323 flg = 0;
01324 while flg == 0,
01325 fprintf('Factor index for No. %i term ', j);
01326 Contr.ord4.cnt(i).code(j).str = input('is (e.g., 01230): ', 's');
01327 if (length(Contr.ord4.cnt(i).code(j).str) ~= NF),
01328 flg = 0; fprintf(2,'\nError: invalid input. Try it again. \n', OutFN);
01329 else flg = 1; end
01330 end
01331 Contr.ord4.cnt(i).coef(j) = input('Corresponding coefficient (e.g., 1 or -1): ');
01332 end
01333 if (sum(Contr.ord4.cnt(i).coef) ~= 0),
01334 flg0 = 0; fprintf(2,'\nError: All the coeficients of a contrast have to sum up to 0! Try again...\n\n');
01335 else flg0 = 1; end
01336 end %flg0 = 0;
01337 end
01338 fprintf(1,'Done with 4th order contrast information.\n');
01339
01340 end % if (NF == 5)
01341
01342
01343 end % end of if (Contrast)
01344
01345 % Done for gathering info from the user.
01346
01347 %=========================================
01348
01349 cov.marker = 10000; % temporary solution when no covariate exists in PreProc.m
01350
01351 if (cov.do), % new design type after covariate is added: Seems this section not needed?
01352 switch NF
01353 case 1,
01354 NF = 2; % 1-ay ANCOVA
01355 dsgn = 1;
01356 case 2,
01357 NF = 3; % 2-ay ANCOVA
01358 if (dsgn == 1), dsgn = 1; end
01359 if (dsgn == 2), dsgn = 2; end
01360 case 3,
01361 NF = 4; % 3-ay ANCOVA
01362 if (dsgn == 1), dsgn = 1; end
01363 if (dsgn == 2), dsgn = 2; end
01364 if (dsgn == 3), dsgn = 3; end
01365 end
01366 end
01367
01368 for (i = 1:1:NF),
01369 if (i == NF & cov.do),
01370 group(NF) = {cov.vec'}; % Convert the cloumn into one row
01371 varnames(NF) = {cov.label};
01372 else
01373 group(i) = {GP(i,:)};
01374 varnames(i) = {FL(i).expr};
01375 end
01376 end
01377
01378 if (cov.do), % swap the covariate with another factor for some special cases: more elegant way to do this?
01379
01380 switch NF
01381 case 2, % 1-way ANCOVA
01382
01383 cov.marker = 2;
01384
01385 for (i = 1:1:Contr.ord1.tot), %1st order contrasts
01386 for (j = 1:1:Contr.ord1.cnt(i).NT),
01387 Contr.ord1.cnt(i).code(j).str = [Contr.ord1.cnt(i).code(j).str(1) '0'];
01388 end
01389 end
01390
01391 case 3, % 2-way ANCOVA
01392
01393 cov.marker = 3; %default
01394
01395 if (dsgn == 2), % for 2-way ANCOVA with design of A(F)XB(R), make covarite the 2nd factor so that 2-way ANCOVA would be 3-way ANOVA of A(F)XB(F)XC(R);
01396 temp1 = group{2};
01397 group{2} = group{3};
01398 group{3} = temp1;
01399
01400 temp2 = varnames(2);
01401 varnames(2) = varnames(3);
01402 varnames(3) = temp2;
01403
01404 temp3 = FL(2);
01405 FL(2) = FL(3);
01406 FL(3) = temp3;
01407
01408 cov.marker = 2; % now 2nd factor is covariate
01409
01410 for (i = 1:1:Contr.ord1.tot),
01411 for (j = 1:1:Contr.ord1.cnt(i).NT),
01412 Contr.ord1.cnt(i).code(j).str = [Contr.ord1.cnt(i).code(j).str(1) '0' Contr.ord1.cnt(i).code(j).str(2)];
01413 end
01414 end
01415
01416 for (i = 1:1:Contr.ord2.tot),
01417 for (j = 1:1:Contr.ord2.cnt(i).NT),
01418 Contr.ord2.cnt(i).code(j).str = [Contr.ord2.cnt(i).code(j).str(1) '0' Contr.ord2.cnt(i).code(j).str(2)];
01419 end
01420 end
01421
01422 else % other two designs
01423 for (i = 1:1:Contr.ord1.tot),
01424 for (j = 1:1:Contr.ord1.cnt(i).NT),
01425 Contr.ord1.cnt(i).code(j).str = [Contr.ord1.cnt(i).code(j).str(1) Contr.ord1.cnt(i).code(j).str(2) '0'];
01426 end
01427 end
01428
01429 for (i = 1:1:Contr.ord2.tot),
01430 for (j = 1:1:Contr.ord2.cnt(i).NT),
01431 Contr.ord2.cnt(i).code(j).str = [Contr.ord2.cnt(i).code(j).str(1) Contr.ord2.cnt(i).code(j).str(2) '0'];
01432 end
01433 end
01434
01435 end % case 3
01436
01437 case 4,
01438
01439 cov.marker = 4; %default
01440
01441 if (dsgn == 2 | dsgn == 3), % for 3-way ANCOVA, swap the 3rd factor and the covariate
01442 temp1 = group{3};
01443 group{3} = group{4};
01444 group{4} = temp1;
01445
01446 temp2 = varnames(3);
01447 varnames(3) = varnames(4);
01448 varnames(4) = temp2;
01449
01450 temp3 = FL(3);
01451 FL(3) = FL(4);
01452 FL(4) = temp3;
01453
01454 cov.marker = 3;
01455
01456 for (i = 1:1:Contr.ord1.tot),
01457 for (j = 1:1:Contr.ord1.cnt(i).NT),
01458 Contr.ord1.cnt(i).code(j).str = [Contr.ord1.cnt(i).code(j).str(1) Contr.ord1.cnt(i).code(j).str(2) '0' Contr.ord1.cnt(i).code(j).str(3)];
01459 end
01460 end
01461
01462 for (i = 1:1:Contr.ord2.tot),
01463 for (j = 1:1:Contr.ord2.cnt(i).NT),
01464 Contr.ord2.cnt(i).code(j).str = [Contr.ord2.cnt(i).code(j).str(1) Contr.ord2.cnt(i).code(j).str(2) '0' Contr.ord2.cnt(i).code(j).str(3)];
01465 end
01466 end
01467
01468 for (i = 1:1:Contr.ord3.tot),
01469 for (j = 1:1:Contr.ord3.cnt(i).NT),
01470 Contr.ord3.cnt(i).code(j).str = [Contr.ord3.cnt(i).code(j).str(1) Contr.ord3.cnt(i).code(j).str(2) '0' Contr.ord3.cnt(i).code(j).str(3)];
01471 end
01472 end
01473
01474 else
01475 for (i = 1:1:Contr.ord1.tot),
01476 for (j = 1:1:Contr.ord1.cnt(i).NT),
01477 Contr.ord1.cnt(i).code(j).str = [Contr.ord1.cnt(i).code(j).str(1) Contr.ord1.cnt(i).code(j).str(2) Contr.ord1.cnt(i).code(j).str(3) '0'];
01478 end
01479 end
01480
01481 for (i = 1:1:Contr.ord2.tot),
01482 for (j = 1:1:Contr.ord2.cnt(i).NT),
01483 Contr.ord2.cnt(i).code(j).str = [Contr.ord2.cnt(i).code(j).str(1) Contr.ord2.cnt(i).code(j).str(2) Contr.ord2.cnt(i).code(j).str(3) '0'];
01484 end
01485 end
01486
01487 for (i = 1:1:Contr.ord3.tot),
01488 for (j = 1:1:Contr.ord3.cnt(i).NT),
01489 Contr.ord3.cnt(i).code(j).str = [Contr.ord3.cnt(i).code(j).str(1) Contr.ord3.cnt(i).code(j).str(2) Contr.ord3.cnt(i).code(j).str(3) '0'];
01490 end
01491 end
01492 end % case 4
01493
01494 end % switch
01495 end
01496
01497 N_Brik = 2^NF - 1; % For 4way with crossed design, there are 15 effect terms: A, B, C, D, AB, AC, AD, BC, BD, ABC, ABD, BCD, and ABCD.
01498 switch NF
01499 case 1,
01500
01501 case 2,
01502
01503 case 3,
01504 if (dsgn == 3 | dsgn == 4), N_Brik = 5; end
01505
01506 case 4,
01507 if (dsgn == 3 | dsgn == 4), N_Brik = 11; end
01508 if (dsgn == 5), N_Brik = 9; end
01509
01510 case 5,
01511 if (dsgn == 3 | dsgn == 4), N_Brik = 23; end
01512
01513 end
01514
01515 %Voxel independent stuff
01516 [err, Qd, s, termname, nterms, sindices, dfbothSS, modw, modwo, tnames, dfterm, dfe, Contr] = PreProc(ntot, NF, group, varnames, FL, Contr, cov, unbalanced);
01517 %[err, Qd, s, termname, nterms, sindices, dfbothSS, modw, modwo, tnames, dfterm, dfe, Contr] = PreProc(ntot, NF, group, varnames, FL, Contr, dsgn, cov);
01518
01519 %Use Ziad's function BrikLoad to load all 'ntot' of datasets in a column format
01520 %since matlab function 'anovan' only accepts vectors.
01521
01522 Opt.Format = 'matrix';
01523
01524 t0 = clock;
01525
01526 [err, FileInfo] = BrikInfo(file(1).nm);
01527 if (data_type == 0),
01528 slices = FileInfo.DATASET_DIMENSIONS(3); % Get the number of slices along Z axis from file header
01529 elseif (data_type == 1),
01530 if (Frame_N == 1), Opt.method = 1; % if being purified
01531 else Opt.method = 2; end
01532 Opt.SliceSize_1D = 50000; % each time run 50000 nodes due to memory limit
01533 slices = ceil(node_n/Opt.SliceSize_1D);
01534 end
01535
01536 fprintf(1, '\nTotal slices along Z axis: %i - You can estimate the total runtime\n', slices);
01537 fprintf(1, '\tby multiplying the runtime for each slice from down below.\n');
01538 fprintf(1, '\nRunning analysis on slice:\n');
01539
01540 for (sn = 1:1:slices),
01541 tic,
01542
01543 fprintf('\t#%d... ', sn);
01544
01545 if (file_format == 1), % Each file contains only single subbrik
01546 % fprintf(1,'\nReading %d input files...', ntot);
01547 Opt.Slices = sn; % Right now just try once slice
01548 % Opt.Frames = 1;
01549 for (i=1:1:ntot), % Read in each file
01550 %V(i) is supposed to be a N X M X K matrix? The information regarding the dimensions
01551 %should be in Info: Info.DATASET_DIMENSIONS(1), Info.DATASET_DIMENSIONS(2),
01552 %Info.DATASET_DIMENSIONS(3), Info.DATASET_RANK(2)?
01553 % [err, X(:, :, :, i), Info, ErrMessage] = BrikLoad(file(i).nm, Opt);
01554
01555 if (data_type == 0), % Volume data
01556 [err, X(:, :, i), Info, ErrMessage] = BrikLoad(file(i).nm, Opt);
01557 elseif (data_type == 1), % Surface data
01558 Opt.Frames = Frame_N;
01559 [err, Z, Info, ErrMessage] = BrikLoad(file(i).nm, Opt);
01560 if (sn ~= 1), X(:,1,i) = zeros(size(X(:,1,i))); end
01561 % For the last chunk of nodes, which are not a whole set of 50,000.
01562 % Not an elegant solution: This would create some dangling 0's in the output file!
01563 X(1:size(Z,1), 1, i) = Z;
01564 clear Z;
01565 end
01566
01567 if (err == 1),
01568 fprintf(2,'Error: failure on loading file %s -- %s. \n', file(i).nm, ErrMessage);
01569 while (1); fprintf(2,'Halted: Ctrl+c to exit'); pause; end
01570 end
01571 end
01572 end
01573
01574 if (file_format == 0), % Each file contains mutliple subbriks from one subject
01575 % fprintf(1,'\nReading slice #%d from %d input files... ', sn, file_num);
01576 Opt.Slices = sn; % Right now just try once slice
01577 % Create a subbrik array. Here we assume the first file_SB subbriks are for ANOVA. Need to change for general situation?
01578 for (i = 1:1:file_SB), Opt.Frames = [Opt.Frames, i]; end
01579
01580 for (i=1:1:file_num), % Read in each file. For each subject (file), there should have FL(1).N_level*FL(4).N_level files
01581 % [err, tmp(:, :, :, (i-1)*file_SB+1:i*file_SB), Info, ErrMessage] = BrikLoad(file(i).nm,Opt);
01582 [err, tmp(:, :, (i-1)*file_SB+1:i*file_SB), Info, ErrMessage] = BrikLoad(file(i).nm,Opt);
01583 if (err == 1),
01584 fprintf(2,'Error: failure on loading file %s -- %s. \n', file(i).nm, ErrMessage);
01585 while (1); fprintf(2,'Halted: Ctrl+c to exit'); pause; end
01586 end
01587 % for (j = 1:1:file_SB), % For each subbrik, there should have FL(2).N_level*FL(3).N_level subbriks
01588 % X(:, :, :, (file(i).F1L-1)*file_SB*FL(4).N_level + (SB(j).lv(1)-1)*FL(3).N_level*FL(4).N_level + ...
01589 % (SB(j).lv(2)-1)*FL(4).N_level + file(i).cntr(file(i).F1L)) = tmp(:, :, :, (i-1)*file_SB+j);
01590
01591 if (NF == 4),
01592 for (j = 1:1:file_SB), % For each subbrik, there should have FL(2).N_level*FL(3).N_level subbriks
01593 X(:, :, (file(i).F1L-1)*file_SB*FL(4).N_level + (SB(j).lv(1)-1)*FL(3).N_level*FL(4).N_level + ...
01594 (SB(j).lv(2)-1)*FL(4).N_level + file(i).cntr(file(i).F1L)) = tmp(:, :, (i-1)*file_SB+j);
01595 end
01596 else fprintf(2,'Sorry, I cannot handle mutliple subbriks for this case. \n');
01597 while (1); fprintf(2,'Halted: Ctrl+c to exit'); pause; end
01598 end
01599
01600 end
01601
01602 % Here we run 4-way ANOVA one slice a time due to swap memory issue. In matrix tmp, the first two
01603 % dimensions are voxels along X and Y axes while the 4th dimension is subbrik tacked together for
01604 % each subject. We need to shuffle the 4th dimension in tmp so that the 4th dimention is arranged
01605 % in the following order: varying 4th factor first, then 3rd, then 2nd and 1st factors.
01606
01607 % 1st factor: (file(i).F1L-1)*file_SB*FL(4).N_level
01608 % 2nd factor: (SB(j).lv(1)-1)*FL(3).N_level*FL(4).N_level
01609 % 3rd factor: (SB(j).lv(2)-1)*FL(4).N_level
01610 % 4th factor: file(i).cntr(file(i).F1L)
01611
01612 % X(:, :, :, (i-1)*file_SB*FL(1).N_level + (j-1)*FL(3).N_level*FL(4).N_level + (k-1)*FL(4).N_level + l)
01613 % Here i, j, k, and l are indeces for 1st, 2nd, 3rd and 4th factor levels.
01614
01615 % file_SB = FL(3).N_level*FL(2).N_level Maybe I should set a checking condition to make sure this equality does hold?
01616
01617 % Is there a better way to load those subbriks directly into matrix X instead of relaying them into
01618 % matrix tmp???
01619
01620 end
01621
01622 D1 = Info.DATASET_DIMENSIONS(1); % X
01623 D2 = Info.DATASET_DIMENSIONS(2); % Y
01624 D3 = Info.DATASET_DIMENSIONS(3); % Z
01625 dim = size(X);
01626
01627 %fstat = zeros(D1, D2, D3, N_Brik);
01628 fstat = zeros(D1, D2, N_Brik);
01629 %intensity = zeros(D1, D2, D3, N_Brik);
01630 intensity = zeros(D1, D2, N_Brik);
01631
01632 %Y = reshape(X, ntot, D1, D2); % Have not figured out how to change the shape of the matrix yet!
01633 for (i = 1:1:D1),
01634 for (j = 1:1:D2),
01635 %for (k = 1:1:D3),
01636 % [err,fstat(i, j, k, :), intensity(i, j, k, :), dfterm_new, dfdenom, tnames_new] = SumsOfSquares(reshape(X(i, j, k, :), ntot, 1), ...
01637 % nterms, Qd, s, sindices, dfbothSS, modw, modwo, tnames, dfterm, dfe, Nest);
01638 [err, fstat(i, j, :), intensity(i, j, :), dfterm_new, dfdenom, tnames_new, LC(i, j)] = SumsOfSquares(reshape(X(i, j, :), ntot, 1), ...
01639 NF, FL, ntot, nterms, Qd, s, sindices, dfbothSS, modw, modwo, tnames, dfterm, dfe, dsgn, N_Brik, Contr);
01640 %end
01641 end
01642 end
01643
01644
01645 %====================================
01646
01647 %M = zeros(D1, D2, D3, 2*N_Brik); %initialization
01648
01649 % if (cov.do == 1),
01650 % M = zeros(D1, D2, 2*(N_Brik+Contr.ord1.tot+Contr.ord2.tot+Contr.ord3.tot+1)); %initialization
01651 % else M = zeros(D1, D2, 2*(N_Brik+Contr.ord1.tot+Contr.ord2.tot+Contr.ord3.tot)); %initialization
01652 % end
01653
01654 M = zeros(D1, D2, 2*(N_Brik+Contr.ord1.tot+Contr.ord2.tot+Contr.ord3.tot+Contr.ord4.tot)); %initialization
01655
01656 % Assemble a matrix with intensity and F value sandwiched with each other
01657 for (i = 1:1:D1),
01658 for (j = 1:1:D2),
01659 %for (k = 1:1:D3),
01660 for (l = 1:1:N_Brik),
01661 % M(i, j, k, 2*l-1) = intensity(i, j, k, l)*isfinite(intensity(i, j, k, l));
01662 M(i, j, 2*l-1) = intensity(i, j, l)*isfinite(intensity(i, j, l));
01663 % M(i, j, k, 2*l) = fstat(i, j, k, l)*isfinite(fstat(i, j, k, l));
01664 M(i, j, 2*l) = fstat(i, j, l)*isfinite(fstat(i, j, l));
01665 end
01666 if (Contr.ord1.tot>0),
01667 for (l = 1:1:Contr.ord1.tot),
01668 M(i, j, 2*(l+N_Brik)-1) = LC(i, j).t1(l).value*isfinite(LC(i, j).t1(l).value);
01669 M(i, j, 2*(l+N_Brik)) = LC(i, j).t1(l).t*isfinite(LC(i, j).t1(l).t);
01670 end
01671 end
01672 if (Contr.ord2.tot>0),
01673 for (l = 1:1:Contr.ord2.tot),
01674 M(i, j, 2*(l+N_Brik+Contr.ord1.tot)-1) = LC(i, j).t2(l).value*isfinite(LC(i, j).t2(l).value);
01675 M(i, j, 2*(l+N_Brik+Contr.ord1.tot)) = LC(i, j).t2(l).t*isfinite(LC(i, j).t2(l).t);
01676 end
01677 end
01678 if (Contr.ord3.tot>0),
01679 for (l = 1:1:Contr.ord3.tot),
01680 M(i, j, 2*(l+N_Brik+Contr.ord1.tot+Contr.ord2.tot)-1) = LC(i, j).t3(l).value*isfinite(LC(i, j).t3(l).value);
01681 M(i, j, 2*(l+N_Brik+Contr.ord1.tot+Contr.ord2.tot)) = LC(i, j).t3(l).t*isfinite(LC(i, j).t3(l).t);
01682 end
01683 end
01684 if (Contr.ord4.tot>0),
01685 for (l = 1:1:Contr.ord4.tot),
01686 M(i, j, 2*(l+N_Brik+Contr.ord1.tot+Contr.ord2.tot+Contr.ord3.tot)-1) = LC(i, j).t4(l).value*isfinite(LC(i, j).t4(l).value);
01687 M(i, j, 2*(l+N_Brik+Contr.ord1.tot+Contr.ord2.tot+Contr.ord3.tot)) = LC(i, j).t4(l).t*isfinite(LC(i, j).t4(l).t);
01688 end
01689 end
01690 end
01691 end
01692
01693 % Get df for 1st order contrasts
01694 % if ((NF == 3 | NF == 4) & Contr.ord1.tot > 0),
01695 if (Contr.ord1.tot > 0),
01696 for (i = 1:1:Contr.ord1.tot),
01697 Contr.ord1.df(i) = dfdenom(Contr.ord1.cnt(i).idx1);
01698 end
01699 end
01700
01701 % Get df for 2nd order contrasts
01702
01703 if ((NF == 3) & Contr.ord2.tot > 0),
01704 for (i = 1:1:Contr.ord2.tot),
01705 switch Contr.ord2.cnt(i).idx1 % 1st factor whose level is fixed
01706 case 1,
01707 switch Contr.ord2.cnt(i).idx2 % 2nd factor whose level is fixed
01708 case 2, Contr.ord2.df(i) = dfdenom(4); % MSAB
01709 case 3, Contr.ord2.df(i) = dfdenom(5) * (dsgn == 1 | dsgn == 2) + dfdenom(3) * (dsgn == 4); % MSAC
01710 end
01711 case 2,
01712 if (Contr.ord2.cnt(i).idx2 == 3), Contr.ord2.df(i) = dfdenom(6) * (dsgn == 1 | dsgn == 2) + dfdenom(5) * (dsgn == 3 | dsgn == 4); % MSBC
01713 else fprintf('\nSomething is wrong in the contrast coding!\n'); fprintf(2,'Halted: Ctrl+c to exit'); pause;
01714 end
01715 case 3,
01716 fprintf('\nSomething is wrong in the contrast coding!\n'); fprintf(2,'Halted: Ctrl+c to exit'); pause;
01717 end
01718 end
01719 end
01720
01721 if ((NF == 4) & Contr.ord2.tot > 0),
01722 for (i = 1:1:Contr.ord2.tot),
01723 % if (dsgn == 3),
01724 switch Contr.ord2.cnt(i).idx1
01725 case 1,
01726 switch Contr.ord2.cnt(i).idx2
01727 case 2, Contr.ord2.df(i) = dfdenom(5); % MSAB
01728 case 3, Contr.ord2.df(i) = dfdenom(6); % MSAC
01729 case 4, Contr.ord2.df(i) = dfdenom(7)*(dsgn == 1 | dsgn == 2) + dfdenom(4) * (dsgn == 4); % MSAD
01730 end
01731 case 2,
01732 switch Contr.ord2.cnt(i).idx2
01733 case 3, Contr.ord2.df(i) = dfdenom(8) * (dsgn == 1 | dsgn == 2) + dfdenom(7) * (dsgn == 3 | dsgn == 4 | dsgn == 5); % MSBC
01734 case 4, Contr.ord2.df(i) = dfdenom(9) * (dsgn == 1 | dsgn == 2) + dfdenom(8) * (dsgn == 3 | dsgn == 4 | dsgn == 5); % Less likely occur: MSBD
01735 end
01736 case 3, % Less likely occur
01737 switch Contr.ord2.cnt(i).idx2
01738 case 4, Contr.ord2.df(i) = dfdenom(10) * (dsgn == 1 | dsgn == 2) + dfdenom(9)*(dsgn == 3 | dsgn == 4 | dsgn == 5); % Less likely occur: MSCD
01739 end
01740 end
01741 % end
01742 end
01743 end
01744
01745 if ((NF == 5) & Contr.ord2.tot > 0), % Refer SumsOfSquares.m for details
01746 for (i = 1:1:Contr.ord2.tot),
01747 switch Contr.ord2.cnt(i).idx1
01748 case 1,
01749 switch Contr.ord2.cnt(i).idx2
01750 case 2, Contr.ord2.df(i) = dfdenom(6); % MSAB
01751 case 3, Contr.ord2.df(i) = dfdenom(7); % MSAC
01752 case 4, Contr.ord2.df(i) = dfdenom(8); % MSAD
01753 case 5, Contr.ord2.df(i) = dfdenom(9)*(dsgn == 1 | dsgn == 2); % MSAD
01754 end
01755 case 2,
01756 switch Contr.ord2.cnt(i).idx2
01757 case 3, Contr.ord2.df(i) = dfdenom(10) * (dsgn == 1 | dsgn == 2) + dfdenom(9) * (dsgn == 3 | dsgn == 4); % MSBC
01758 case 4, Contr.ord2.df(i) = dfdenom(11) * (dsgn == 1 | dsgn == 2) + dfdenom(10) * (dsgn == 3 | dsgn == 4); % Less likely occur: MSBD
01759 case 5, Contr.ord2.df(i) = dfdenom(12) * (dsgn == 1 | dsgn == 2) + dfdenom(11) * (dsgn == 3 | dsgn == 4); % Less likely occur: MSBE
01760 end
01761 case 3, % Less likely occur
01762 switch Contr.ord2.cnt(i).idx2
01763 case 4, Contr.ord2.df(i) = dfdenom(13) * (dsgn == 1 | dsgn == 2) + dfdenom(12)*(dsgn == 3 | dsgn == 4); % Less likely occur: MSCD
01764 case 5, Contr.ord2.df(i) = dfdenom(14) * (dsgn == 1 | dsgn == 2) + dfdenom(13)*(dsgn == 3 | dsgn == 4); % Less likely occur: MSCE
01765 end
01766 case 4, % Less likely occur
01767 switch Contr.ord2.cnt(i).idx2
01768 case 5, Contr.ord2.df(i) = dfdenom(15) * (dsgn == 1 | dsgn == 2) + dfdenom(14)*(dsgn == 3 | dsgn == 4); % Less likely occur: MSDE
01769 end
01770 end % switch Contr.ord2.cnt(i).idx1
01771 % end
01772 end
01773 end % if ((NF == 5) & Contr.ord2.tot > 0)
01774
01775
01776
01777 % Get df for 3rd order contrasts
01778
01779 if ((NF == 3) & Contr.ord3.tot > 0),
01780 for (i = 1:1:Contr.ord3.tot),
01781 switch Contr.ord3.cnt(i).idx1
01782 case 1,
01783 switch Contr.ord3.cnt(i).idx2
01784 case 2,
01785 if (Contr.ord3.cnt(i).idx3 == 3), Contr.ord3.df(i) = dfdenom(7)*(dsgn == 1 | dsgn == 2); % MSABC
01786 else fprintf('\nSomething is wrong in the contrast coding!\n'); fprintf(2,'Halted: Ctrl+c to exit'); pause;
01787 end
01788 case 3, fprintf('\nSomething is wrong in the contrast coding!\n'); fprintf(2,'Halted: Ctrl+c to exit'); pause;
01789 end
01790 case 2, fprintf('\nSomething is wrong in the contrast coding!\n'); fprintf(2,'Halted: Ctrl+c to exit'); pause;
01791 case 3, fprintf('\nSomething is wrong in the contrast coding!\n'); fprintf(2,'Halted: Ctrl+c to exit'); pause;
01792 end
01793 end
01794 end
01795
01796 if ((NF == 4) & Contr.ord3.tot > 0),
01797 for (i = 1:1:Contr.ord3.tot),
01798 % if (dsgn == 3),
01799 switch Contr.ord3.cnt(i).idx1
01800 case 1,
01801 switch Contr.ord3.cnt(i).idx2
01802 case 2,
01803 switch Contr.ord3.cnt(i).idx3
01804 case 3, Contr.ord3.df(i) = dfdenom(11)*(dsgn == 1 | dsgn == 2) + dfdenom(10)* (dsgn == 3 | dsgn == 4 | dsgn == 5); % MSABC
01805 case 4, Contr.ord3.df(i) = dfdenom(12)*(dsgn == 1 | dsgn == 2); % MSABD not exist for (dsgn == 3 | dsgn == 4 | dsgn == 5)
01806 end
01807 case 3,
01808 if (Contr.ord3.cnt(i).idx3 == 4), Contr.ord3.df(i) = dfdenom(13)*(dsgn == 1 | dsgn == 2); % MSACD
01809 else fprintf('\nSomething is wrong in the contrast coding!\n'); fprintf(2,'Halted: Ctrl+c to exit'); pause;
01810 end
01811 end
01812 case 2,
01813 switch Contr.ord3.cnt(i).idx2
01814 case 3,
01815 if (Contr.ord3.cnt(i).idx3 == 4), Contr.ord3.df(i) = dfdenom(14)*(dsgn == 1 | dsgn == 2) + dfdenom(11)* (dsgn == 3 | dsgn == 4 | dsgn == 5); % MSBCD
01816 else fprintf('\nSomething is wrong in the contrast coding!\n'); fprintf(2,'Halted: Ctrl+c to exit'); pause;end
01817 case 4,
01818 fprintf('\nSomething is wrong in the contrast coding!\n'); fprintf(2,'Halted: Ctrl+c to exit'); pause;
01819 end
01820 case 3,
01821 fprintf('\nSomething is wrong in the contrast coding!\n');
01822 fprintf(2,'Halted: Ctrl+c to exit'); pause;
01823 case 4,
01824 fprintf('\nSomething is wrong in the contrast coding!\n');
01825 fprintf(2,'Halted: Ctrl+c to exit'); pause;
01826 end
01827 % end
01828 end
01829 end
01830
01831 if ((NF == 5) & Contr.ord3.tot > 0),
01832 for (i = 1:1:Contr.ord3.tot),
01833 switch Contr.ord3.cnt(i).idx1
01834 case 1,
01835 switch Contr.ord3.cnt(i).idx2
01836 case 2,
01837 switch Contr.ord3.cnt(i).idx3
01838 case 3, Contr.ord3.df(i) = dfdenom(16)*(dsgn == 1 | dsgn == 2) + dfdenom(15)* (dsgn == 3 | dsgn == 4); % MSABC
01839 case 4, Contr.ord3.df(i) = dfdenom(17)*(dsgn == 1 | dsgn == 2) + dfdenom(16)* (dsgn == 3 | dsgn == 4); % MSABD
01840 case 5, Contr.ord3.df(i) = dfdenom(18)*(dsgn == 1 | dsgn == 2); % MSABE
01841 end
01842 case 3,
01843 switch Contr.ord3.cnt(i).idx3
01844 case 4, Contr.ord3.df(i) = dfdenom(19)*(dsgn == 1 | dsgn == 2) + dfdenom(17)* (dsgn == 3 | dsgn == 4); % MSACD
01845 case 5, Contr.ord3.df(i) = dfdenom(20)*(dsgn == 1 | dsgn == 2); % MSACE
01846 end
01847 case 4,
01848 if (Contr.ord3.cnt(i).idx3 == 5), Contr.ord3.df(i) = dfdenom(21)*(dsgn == 1 | dsgn == 2); % MSADE
01849 else fprintf('\nSomething is wrong in the contrast coding!\n'); fprintf(2,'Halted: Ctrl+c to exit'); pause;
01850 end
01851 end % switch Contr.ord3.cnt(i).idx2
01852 case 2,
01853 switch Contr.ord3.cnt(i).idx2
01854 case 3,
01855 switch Contr.ord3.cnt(i).idx3
01856 case 4, Contr.ord3.df(i) = dfdenom(22)*(dsgn == 1 | dsgn == 2) + dfdenom(18)* (dsgn == 3 | dsgn == 4); % MSBCD
01857 case 5, Contr.ord3.df(i) = dfdenom(23)*(dsgn == 1 | dsgn == 2) + dfdenom(19)* (dsgn == 3 | dsgn == 4); % MSBCE
01858 end
01859 case 4,
01860 if (Contr.ord3.cnt(i).idx3 == 5), Contr.ord3.df(i) = dfdenom(24)*(dsgn == 1 | dsgn == 2) + dfdenom(20)* (dsgn == 3 | dsgn == 4); % MSBDE
01861 else fprintf('\nSomething is wrong in the contrast coding!\n'); fprintf(2,'Halted: Ctrl+c to exit'); pause;end
01862 case 4,
01863 fprintf('\nSomething is wrong in the contrast coding!\n'); fprintf(2,'Halted: Ctrl+c to exit'); pause;
01864 end % switch Contr.ord3.cnt(i).idx2
01865 case 3,
01866 if (Contr.ord3.cnt(i).idx2 == 4 & Contr.ord3.cnt(i).idx3 == 5),
01867 if (dsgn == 1 | dsgn == 2),
01868 Contr.ord3.df(i) = dfdenom(25);
01869 elseif (dsgn == 3),
01870 dfdenom(21); % MSCDE
01871 end
01872 else fprintf('\nSomething is wrong in the contrast coding!\n'); fprintf(2,'Halted: Ctrl+c to exit'); pause;end
01873
01874 case 4,
01875 fprintf('\nSomething is wrong in the contrast coding!\n');
01876 fprintf(2,'Halted: Ctrl+c to exit'); pause;
01877 end % switch Contr.ord3.cnt(i).idx1
01878 end
01879 end
01880
01881 if ((NF == 5) & Contr.ord4.tot > 0),
01882 for (i = 1:1:Contr.ord4.tot),
01883 switch Contr.ord4.cnt(i).idx1
01884 case 1,
01885 switch Contr.ord4.cnt(i).idx2
01886 case 2,
01887 switch Contr.ord4.cnt(i).idx3
01888 case 3,
01889 switch Contr.ord4.cnt(i).idx4
01890 case 4,
01891 if (dsgn == 1 | dsgn == 2),
01892 Contr.ord4.df(i) = dfdenom(26);
01893 elseif (dsgn == 3 | dsgn == 4),
01894 Contr.ord4.df(i) = dfdenom(22); % MSABCD
01895 end
01896 case 5,
01897 if (dsgn == 1 | dsgn == 2),
01898 Contr.ord4.df(i) = dfdenom(27);
01899 end % MSABCE
01900 end % switch Contr.ord4.cnt(i).idx4
01901 case 4,
01902 if (Contr.ord4.cnt(i).idx4 == 5),
01903 if (dsgn == 1 | dsgn == 2),
01904 Contr.ord4.df(i) = dfdenom(28); % MSABDE
01905 end
01906 else fprintf('\nSomething is wrong in the contrast coding!\n'); fprintf(2,'Halted: Ctrl+c to exit'); pause;
01907 end % if (Contr.ord4.cnt(i).idx4 == 5)
01908 case 5, fprintf('\nSomething is wrong in the contrast coding!\n');
01909 fprintf(2,'Halted: Ctrl+c to exit'); pause;
01910 end % switch Contr.ord4.cnt(i).idx3
01911 case 3,
01912 switch Contr.ord4.cnt(i).idx3
01913 case 4,
01914 if (Contr.ord4.cnt(i).idx4 == 5),
01915 if (dsgn == 1 | dsgn == 2),
01916 Contr.ord4.df(i) = dfdenom(29); % MSACDE
01917 end
01918 else fprintf('\nSomething is wrong in the contrast coding!\n'); fprintf(2,'Halted: Ctrl+c to exit'); pause;
01919 end
01920 case 5, fprintf('\nSomething is wrong in the contrast coding!\n');
01921 fprintf(2,'Halted: Ctrl+c to exit'); pause;
01922 end % switch Contr.ord4.cnt(i).idx3
01923 case 4, fprintf('\nSomething is wrong in the contrast coding!\n');
01924 fprintf(2,'Halted: Ctrl+c to exit'); pause;
01925 end % switch Contr.ord4.cnt(i).idx2
01926
01927 case 2,
01928 switch Contr.ord4.cnt(i).idx2
01929 case 3,
01930 switch Contr.ord4.cnt(i).idx3
01931 case 4,
01932 if (Contr.ord4.cnt(i).idx4 == 5),
01933 if (dsgn == 1 | dsgn == 2),
01934 Contr.ord4.df(i) = dfdenom(30); % MSBCDE
01935 end
01936 else fprintf('\nSomething is wrong in the contrast coding!\n'); fprintf(2,'Halted: Ctrl+c to exit'); pause;
01937 end
01938 case 5, fprintf('\nSomething is wrong in the contrast coding!\n');
01939 fprintf(2,'Halted: Ctrl+c to exit'); pause;
01940 end % switch Contr.ord4.cnt(i).idx3
01941 case 4,
01942 fprintf('\nSomething is wrong in the contrast coding!\n'); fprintf(2,'Halted: Ctrl+c to exit'); pause;
01943 end % switch Contr.ord4.cnt(i).idx2
01944 case {3,4,5},
01945 fprintf('\nSomething is wrong in the contrast coding!\n');
01946 fprintf(2,'Halted: Ctrl+c to exit'); pause;
01947
01948 end % switch Contr.ord4.cnt(i).idx1
01949 end % for (i = 1:1:Contr.ord4.tot)
01950 end % if ((NF == 5) & Contr.ord4.tot > 0)
01951
01952 % if (cov.do == 1), covoutdf = dfe; end % df of MSE
01953
01954 %Set up options for writing out the results
01955
01956 Opt.Prefix = sprintf('%s', OutFN);
01957 Opt.NoCheck = 0;
01958 Opt.Scale = 1;
01959 if (data_type == 0), Opt.View = sprintf('+%s', format);
01960 elseif (data_type == 1) Opt.View = '.1D.dset'; end
01961
01962 Info.TYPESTRING = '3DIM_HEAD_FUNC';
01963 %Info.DATASET_RANK(2) = 2*N_Brik; % Number of subbriks
01964 % if (cov.do == 1),
01965 % Info.DATASET_RANK(2) = 2*(N_Brik + Contr.ord1.tot + Contr.ord2.tot + Contr.ord3.tot +1);
01966 % else Info.DATASET_RANK(2) = 2*(N_Brik + Contr.ord1.tot + Contr.ord2.tot + Contr.ord3.tot); % Number of subbriks
01967 % end
01968 Info.DATASET_RANK(2) = 2*(N_Brik + Contr.ord1.tot + Contr.ord2.tot + Contr.ord3.tot + Contr.ord4.tot); % Number of subbriks
01969
01970 Info.DATASET_DIMENSIONS = [D1 D2 D3 0 0];
01971 Info.BRICK_STATS = [];
01972 % Info.BRICK_TYPES = [];
01973 Info.BRICK_FLOAT_FACS = [];
01974 Info.BRICK_STATAUX = []; %initialization
01975 Info.BRICK_TYPES = []; %initialization: generate results in float with 3
01976 Info.TypeBytes = 4; % Force the output to be float. Why?
01977
01978 % For all terms
01979 for (i = 1:1:N_Brik),
01980 Info.BRICK_STATAUX = [Info.BRICK_STATAUX 2*i-1 4 2 dfterm_new(i) dfdenom(i)]; % 4 is for F stat defined in 3ddata.h
01981 Info.BRICK_TYPES = [Info.BRICK_TYPES 3 3]; % 3 for float; 1 for short
01982 end
01983
01984 % For 1st oder contrasts
01985 if (Contr.ord1.tot>0),
01986 for (i = 1:1:Contr.ord1.tot),
01987 Info.BRICK_STATAUX = [Info.BRICK_STATAUX 2*(i+N_Brik)-1 3 1 Contr.ord1.df(i)]; % 3 is for t stat defined in 3ddata.h
01988 Info.BRICK_TYPES = [Info.BRICK_TYPES 3 3]; % 3 for float; 1 for short
01989 end
01990 end
01991
01992 % For 2nd oder contrasts
01993 if (Contr.ord2.tot>0),
01994 for (i = 1:1:Contr.ord2.tot),
01995 Info.BRICK_STATAUX = [Info.BRICK_STATAUX 2*(i+N_Brik+Contr.ord1.tot)-1 3 1 Contr.ord2.df(i)]; % 3 is for t stat defined in 3ddata.h
01996 Info.BRICK_TYPES = [Info.BRICK_TYPES 3 3]; % 3 for float; 1 for short
01997 end
01998 end
01999
02000 % For 3rd oder contrasts
02001 if (Contr.ord3.tot>0),
02002 for (i = 1:1:Contr.ord3.tot),
02003 Info.BRICK_STATAUX = [Info.BRICK_STATAUX 2*(i+N_Brik+Contr.ord1.tot+Contr.ord2.tot)-1 3 1 Contr.ord3.df(i)]; % 3 is for t stat defined in 3ddata.h
02004 Info.BRICK_TYPES = [Info.BRICK_TYPES 3 3]; % 3 for float; 1 for short
02005 end
02006 end
02007
02008 % For 3rd oder contrasts
02009 if (Contr.ord4.tot>0),
02010 for (i = 1:1:Contr.ord4.tot),
02011 Info.BRICK_STATAUX = [Info.BRICK_STATAUX 2*(i+N_Brik+Contr.ord1.tot+Contr.ord2.tot+Contr.ord3.tot)-1 3 1 Contr.ord4.df(i)]; % 3 is for t stat defined in 3ddata.h
02012 Info.BRICK_TYPES = [Info.BRICK_TYPES 3 3]; % 3 for float; 1 for short
02013 end
02014 end
02015
02016 % ANOCVA
02017 % if (cov.do == 1),
02018 % Info.BRICK_STATAUX = [Info.BRICK_STATAUX 2*(i+N_Brik+Contr.ord1.tot+Contr.ord2.tot+Contr.ord3.tot)-1 3 1 covoutdf]; % 3 is for t stat defined in 3ddata.h
02019 % Info.BRICK_TYPES = [Info.BRICK_TYPES 3 3]; % 3 for float; 1 for short
02020 % end
02021
02022 % For all terms
02023 Info.BRICK_LABS = [];
02024 for (i = 1:1:N_Brik), %create labels for all terms
02025 Info.BRICK_LABS = [Info.BRICK_LABS cell2mat(tnames_new(i)) '~']; % Label for intensity
02026 Info.BRICK_LABS = [Info.BRICK_LABS cell2mat(tnames_new(i)) ' ' 'F' '~']; % Label for F
02027 end
02028
02029 % For 1st order contrasts
02030 if (Contr.ord1.tot>0),
02031 for (i = 1:1:Contr.ord1.tot), %create labels for contrasts
02032 Info.BRICK_LABS = [Info.BRICK_LABS Contr.ord1.label(i).nm '~']; % Label for intensity
02033 Info.BRICK_LABS = [Info.BRICK_LABS Contr.ord1.label(i).nm ' ' 't' '~']; % Label for t
02034 end
02035 end
02036
02037 % For 2nd order contrasts
02038 if (Contr.ord2.tot>0),
02039 for (i = 1:1:Contr.ord2.tot), %create labels for contrasts
02040 Info.BRICK_LABS = [Info.BRICK_LABS Contr.ord2.label(i).nm '~']; % Label for intensity
02041 Info.BRICK_LABS = [Info.BRICK_LABS Contr.ord2.label(i).nm ' ' 't' '~']; % Label for t
02042 end
02043 end
02044
02045 % For 3rd order contrasts
02046 if (Contr.ord3.tot>0),
02047 for (i = 1:1:Contr.ord3.tot), %create labels for contrasts
02048 Info.BRICK_LABS = [Info.BRICK_LABS Contr.ord3.label(i).nm '~']; % Label for intensity
02049 Info.BRICK_LABS = [Info.BRICK_LABS Contr.ord3.label(i).nm ' ' 't' '~']; % Label for t
02050 end
02051 end
02052
02053 % For 4th order contrasts
02054 if (Contr.ord4.tot>0),
02055 for (i = 1:1:Contr.ord4.tot), %create labels for contrasts
02056 Info.BRICK_LABS = [Info.BRICK_LABS Contr.ord4.label(i).nm '~']; % Label for intensity
02057 Info.BRICK_LABS = [Info.BRICK_LABS Contr.ord4.label(i).nm ' ' 't' '~']; % Label for t
02058 end
02059 end
02060
02061 % ANOCVA
02062 % if (cov.do == 1),
02063 % Info.BRICK_LABS = [Info.BRICK_LABS cov.label '~']; % Label for beta
02064 % Info.BRICK_LABS = [Info.BRICK_LABS cov.label ' ' 't' '~']; % Label for t
02065 % end
02066
02067 Info.SCENE_DATA(2) = 11; % 11 means Bucket Func type, defined in AFNI doc README.attributes
02068 Info.SCENE_DATA(3) = 1; % 1 means 3DIM_HEAD_FUNC, defined in AFNI doc README.attributes
02069
02070 Opt.Frames = []; %Because it might have been set as the frame list in the case of input files with multiple subbriks during loading
02071
02072 % Write output: reshape because the 3rd dimension is supposed to be Z in the normal situation.
02073 if (data_type == 0),
02074 [err2, ErrMessage, NewInfo] = WriteBrik(reshape(M, D1, D2, 1, Info.DATASET_RANK(2)), Info, Opt);
02075 elseif (data_type == 1), % Collapse the 2nd dimsion, which is 1 for surface data.
02076 [err2, ErrMessage, NewInfo] = WriteBrik(squeeze(reshape(M, D1, D2, 1, Info.DATASET_RANK(2))), Info, Opt);
02077 end
02078 fprintf(1, 'done in %f seconds\n', toc);
02079
02080 end % end of the big loop of running ANOVA and writing up: One slice a time
02081
02082 fprintf(1, '\nCongratulations, job is successfully done!!! Total runtime: %f minutes...', etime(clock,t0)/60);
02083 fprintf(1, '\nOutput files are %s%s.*\n\n', Opt.Prefix, Opt.View);
02084 err = 0;
02085 return;
