00001 function [err, fstat, intensity_new, dfterm_new, dfdenom, tnames_new, LC] = SumsOFSquares(y, NF, FL, ntot, nterms, Qd, s, sindices, dfbothSS, modw, modwo, tnames, dfterm, dfe, dsgn, N_Brik, Contr)
00002 %
00003 % [err,] = ss.m ()
00004 %
00005 %Purpose:
00006 %
00007 %
00008 %
00009 %Input Parameters:
00010 %
00011 % y: values from all factor combinations at one voxel in vector format
00012 % ntot: total number of files = length of vector y
00013 % tername: term names such as
00014 % termlist:
00015 % dmat: design matrix in the converted regression model
00016 % cmat: constraints matrix
00017 %
00018 %Output Parameters:
00019 % err : 0 No Problem
00020 % : 1 Problems
00021 % fstat: a vector of F values for those terms at a voxel
00022 % intensity_new: intensity for those terms at a voxel, which are defined as the sqare root of MS terms
00023 % dfterm: vector of degree of freedom for those terms (numerators) at each voxel (but the same for all voxels).
00024 % This is why I don't differentiate them among the voxels in anova.m
00025 % dfdenom: vector
00026 %
00027 %Key Terms:
00028 %
00029 %More Info :
00030 %
00031 %
00032 %
00033 %
00034 % Author : Gang Chen
00035 % Date : Tue Mar 23 14:05:05 EST 2004
00036 % SSCC/NIMH/ National Institutes of Health, Bethesda MD 20892
00037
00038
00039 %Define the function name for easy referencing
00040 FuncName = 'ss.m';
00041
00042 %Debug Flag
00043 DBG = 1;
00044
00045 %initailize return variables
00046 err = 1;
00047
00048 dfboth = dfbothSS;
00049 y0 = y; % y gets mutated down below. Keep the original y so that the contrasts can be done later at the end.
00050
00051 % Fit the full model and compute the residual sum of squares
00052 mu = mean(y);
00053 y = y - mu; % two passes to improve accuracy
00054 mu2 = mean(y);
00055 %mu = mu + mu2; % THIS CAN BE DELETED FOR MY CASE since it is not used from this point!!!!!
00056 y = y - mu2;
00057 sst = sum(y.^2); % SSTO
00058
00059 [ssx, y] = dofit(y, Qd); % y: predicted value; ssx:
00060 sse = sst - ssx; %SS(Error) = SS(Total) - SS(Rest)
00061
00062 % Fit each model, get its residual SS and d.f.
00063 ssw(1:nterms) = ssx; % for full model we already know the results
00064
00065 for j=length(sindices):-1:1
00066 % Find the next model index to fit
00067 k = sindices(j);
00068
00069 % Look in unsorted arrays to see if we have already fit this model
00070 if j>nterms
00071 k0 = k+nterms;
00072 else
00073 k0 = k;
00074 end
00075 if dfboth(k0)~=-1
00076 continue
00077 end
00078
00079 % Find the model with this index
00080 if (j > nterms)
00081 thismod = modwo(k, :);
00082 else
00083 thismod = modw(k, :);
00084 end
00085
00086 % Fit this model
00087
00088 % [ssx0,dfx0] = QRfit(X, y, C);
00089 ssx0 = dofit(y, s(j).Qdt);
00090
00091 % Use these results for each term that requires them
00092 mod0 = repmat(thismod, nterms, 1);
00093 k = find(all(modw == mod0, 2));
00094 ssw(k) = ssx0;
00095 k = find(all(modwo == mod0, 2));
00096 sswo(k) = ssx0;
00097 end
00098 clear mod0
00099
00100
00101 % Compute the sum of squares attributed to each term
00102 ssterm = ssw - sswo;
00103 ssterm(dfterm==0) = 0; % make this exact
00104
00105
00106 % ========================================================
00107 % !!!The ones for dfterm are the same for all voxles: is there a way to move them out of this routine?
00108 % Also the switches should stay in PreProc.m instead of being here!
00109 %=========================================================
00110
00111 switch NF
00112 case 1,
00113 fstat = repmat(0, size(ssterm)); % only ss term: SSTR -- treatment sum of squares
00114
00115 case 2,
00116 fstat = repmat(0, size(ssterm));
00117
00118 case 3,
00119 switch dsgn
00120 case 1,
00121 fstat = repmat(0, [1 N_Brik]);
00122
00123 case 2, % 7 terms: 1 (A); 2 (B); 3 (C); 4 (AB), 5 (AC), 6 BC, 7 ABC
00124 fstat = repmat(0, [1 N_Brik]); % same as fstat = repmat(0, [1 N_Brik]);
00125
00126 case {3, 4,} % BXC(A): only 5 terms - 1 (A); 2 (B); 3 C(A); 4 (AB), 5 BC(A)
00127
00128 ssterm(3) = ssterm(3) + ssterm(5); % SSC(A) = SSC + SSAC
00129 dfterm(3) = dfterm(3) + dfterm(5); % DF C(A) = DF C + DF AC
00130
00131 ssterm(6) = ssterm(6) + ssterm(7); % SSBC(A) = SSC + SSAC
00132 dfterm(6) = dfterm(6) + dfterm(7);
00133
00134 fstat = repmat(0, [1 N_Brik]);
00135
00136 end % switch dsgn
00137
00138 case 4,
00139
00140 switch dsgn
00141 % Order of the 15 terms: 1 (A); 2 (B); 3 (C); 4 (D); 5 (AB); 6 (AC); 7 (AD); 8 (BC); 9 (BD); 10 (CD)
00142 % 11 (ABC); 12 (ABD); 13 (ACD); 14 (BCD); 15 (ABCD)
00143 case {1, 2,}
00144 fstat = repmat(0, [1 N_Brik]);
00145 case {3, 4,} % only 11 terms in nesting case without AD, ABD, ACD, and ABCD: 1 (A); 2 (B); 3 (C); 4 (D); 5 (AB); 6 (AC);
00146 % 7 (BC); 8 (BD); 9 (CD); 10 (ABC); 11 (BCD);
00147 ssterm(4) = ssterm(4) + ssterm(7); %SSD(A) = SSD + SSAD
00148 dfterm(4) = dfterm(4) + dfterm(7); %DF(D(A)) = DF(D) + DF(AD)
00149
00150 ssterm(9) = ssterm(9) + ssterm(12); %SSBD(A) = SSB + SSABD
00151 dfterm(9) = dfterm(9) + dfterm(12); %DF(BD(A)) = DF(BD) + DF(ABD)
00152
00153 ssterm(10) = ssterm(10) + ssterm(13); %SSCD(A) = SSC + SSACD
00154 dfterm(10) = dfterm(10) + dfterm(13); %DF(CD(A)) = DF(CD) + DF(ACD)
00155
00156 ssterm(14) = ssterm(14) + ssterm(15); %SSBCD(A) = SSBCD + SSABCD
00157 dfterm(14) = dfterm(14) + dfterm(15); %DF(BCD(A)) = DF(BCD)+ DF(ABCD)
00158
00159 fstat = repmat(0, [1 N_Brik]);
00160
00161 case 5,
00162 % There are only 9 terms: 1 (A); 2 (B); 3 (C); 4 (D); 5 (AB); 6 (AC); 7 (BC); 8 (CD); 9 (ABC)
00163
00164 ssterm(4) = ssterm(4) + ssterm(7) + ssterm(9) + ssterm(12); %SSD(AB) = SSD + SSAD + SSBD + SSABD
00165 dfterm(4) = dfterm(4) + dfterm(7) + dfterm(9) + dfterm(12); %DF(D(AB)) = (DF(A)+ DF(B) + DF(AB) + 1)DF(D)
00166
00167 ssterm(10) = ssterm(10) + ssterm(13) + ssterm(14) + ssterm(15); %SSCD(AB) = SSCD + SSACD + SSBCD + SSABCD
00168 dfterm(10) = dfterm(10) + dfterm(13) + dfterm(14) + dfterm(15); %DF(CD(AB)) = (DF(A)+ DF(B) + DF(AB) + 1)DF(D)
00169
00170 fstat = repmat(0, [1 N_Brik]);
00171 end
00172
00173 case 5,
00174
00175 switch dsgn %Order of the 2^5 - 1 = 31 terms: 1 (A); 2 (B); 3 (C); 4 (D); 5 (E); 6 (AB); 7 (AC); 8 (AD); 9 (AE);
00176 % 10 (BC); 11 (BD); 12 (BE); 13 (CD); 14 (CE); 15 (DE) 16 (ABC); 17 (ABD); 18 (ABE); 19 (ACD); 20 (ACE); 21 (ADE);
00177 % 22 (BCD); 23 (BCE); 24 (BDE); 25 (CDE); 26 (ABCD); 27 (ABCE); 28 (ABDE) 29 (ACDE); 30 (BCDE); 31 (ABCDE)
00178
00179 case {1, 2,}
00180 fstat = repmat(0, [1 N_Brik]); % pure cross design
00181
00182 case {3,4}, % only 23 terms
00183
00184 ssterm(5) = ssterm(5) + ssterm(9); %SSE(A) = SSE + SSAE
00185 dfterm(5) = dfterm(5) + dfterm(9); %DF(E(A)) = DF(E) + DF(AE)
00186
00187 ssterm(12) = ssterm(12) + ssterm(18); %SSBE(A) = SSBE + SSABE
00188 dfterm(12) = dfterm(12) + dfterm(18); %DF(BE(A)) = DF(BE) + DF(ABE)
00189
00190 ssterm(14) = ssterm(14) + ssterm(20); %SSCE(A) = SSCE + SSACE
00191 dfterm(14) = dfterm(14) + dfterm(20); %DF(CE(A)) = DF(CE) + DF(ACE)
00192
00193 ssterm(15) = ssterm(15) + ssterm(21); %SSDE(A) = SSDE + SSADE
00194 dfterm(15) = dfterm(15) + dfterm(21); %DF(DE(A)) = DF(DE) + DF(ADE)
00195
00196 ssterm(23) = ssterm(23) + ssterm(27); %SSBCE(A) = SSBCE + SSABCE
00197 dfterm(23) = dfterm(23) + dfterm(27); %DF(BCE(A)) = DF(BCE) + DF(ABCE)
00198
00199 ssterm(24) = ssterm(24) + ssterm(28); %SSBDE(A) = SSBDE + SSABDE
00200 dfterm(24) = dfterm(24) + dfterm(28); %DF(BDE(A)) = DF(BDE) + DF(ABDE)
00201
00202 ssterm(25) = ssterm(25) + ssterm(29); %SSCDE(A) = SSCDE + SSACDE
00203 dfterm(25) = dfterm(25) + dfterm(29); %DF(CDE(A)) = DF(CDE) + DF(ACDE)
00204
00205 ssterm(30) = ssterm(30) + ssterm(31); %SSBCDE(A) = SSBCDE + SSABCDE
00206 dfterm(30) = dfterm(30) + dfterm(31); %DF(BCDE(A)) = DF(BCDE) + DF(ABCDE)
00207
00208 fstat = repmat(0, [1 N_Brik]);
00209 end % switch dsgn
00210
00211 end % switch NF
00212
00213 % Compute the mean square for each term
00214 msterm = ssterm ./ max(1, dfterm');
00215 mse = sse * (dfe>0) / max(1, dfe);
00216 intensity = sqrt(msterm); %Intensity for each term
00217
00218 % Equated computed and expected mean squares, then solve for
00219 % variance component estimates
00220
00221
00222 switch NF
00223 case 1,
00224 dfterm_new = dfterm; tnames_new = tnames; msterm_new = msterm; intensity_new = intensity;
00225 msdenom = repmat(mse, size(msterm)); dfdenom = repmat(dfe, size(msterm));
00226
00227 case 2, %totally 3 terms: A, B, abd AXB.
00228 dfterm_new = dfterm; tnames_new = tnames; msterm_new = msterm; intensity_new = intensity;
00229
00230 switch dsgn % Allocate denominator and its d.f. for each F ratio
00231 case 1,
00232 msdenom = repmat(mse, size(msterm));
00233 dfdenom = repmat(dfe, size(msterm));
00234 case 2,
00235 msdenom = [msterm(3), mse, mse];
00236 dfdenom = [dfterm(3), dfe, dfe];
00237 case 3,
00238 msdenom = [msterm(3), msterm(3), mse];
00239 dfdenom = [dfterm(3), dfterm(3), dfe];
00240 end % Close swtich dsgn
00241
00242 case 3,
00243 switch dsgn
00244 case 1,
00245 dfterm_new = dfterm; tnames_new = tnames; msterm_new = msterm; intensity_new = intensity;
00246 msdenom = repmat(mse, size(msterm)); dfdenom = repmat(dfe, size(msterm));
00247
00248 case 2, % 7 terms: 1 (A); 2 (B); 3 (C); 4 (AB), 5 (AC), 6 BC, 7 ABC
00249 dfterm_new = dfterm; tnames_new = tnames; msterm_new = msterm; intensity_new = intensity;
00250 msdenom = [msterm(5), msterm(6), mse, msterm(7), mse, mse, mse];
00251 dfdenom = [dfterm(5), dfterm(6), dfe, dfterm(7), dfe, dfe, dfe];
00252
00253 case 3, % 5 terms: 1 A; 2 B; 3 C(A); 4 (AB), 5 BC(A)
00254 msterm_new = [msterm(1:4), msterm(6)]; % Throw out those four which do not exist for nesting: AC and ABC.
00255 intensity_new = [intensity(1:4), intensity(6)]; % Throw out those four which do not exist for nesting.
00256 dfterm_new = [dfterm(1:4)', dfterm(6)'];
00257 tnames_new = [tnames(1:4); tnames(6)]; % Only preserve those valid for nesting. Semicolon for coloumn catenation
00258 msdenom = [msterm(3), msterm(6), mse, msterm(6), mse,0,0]; % pad 2 extra 0's for potential complaints
00259 dfdenom = [dfterm(3), dfterm(6), dfe, dfterm(6), dfe,0,0];
00260
00261 case 4, % 5 terms: 1 A; 2 B; 3 C(A); 4 (AB), 5 BC(A)
00262 msterm_new = [msterm(1:4), msterm(6)]; % Throw out those four which do not exist for nesting: AC and ABC.
00263 intensity_new = [intensity(1:4), intensity(6)]; % Throw out those four which do not exist for nesting.
00264 dfterm_new = [dfterm(1:4)', dfterm(6)'];
00265 tnames_new = [tnames(1:4); tnames(6)]; % Only preserve those valid for nesting. Semicolon for coloumn catenation
00266 msdenom = [msterm(4), mse, msterm(6), mse, mse,0,0]; % 2 extra 0's are for error-prone problem down below for contrasts
00267 dfdenom = [dfterm(4), dfe, dfterm(6), dfe, dfe,0,0];
00268 end
00269
00270 case 4,
00271 switch dsgn %Order of the 2^4 - 1 = 15 terms: 1 (A); 2 (B); 3 (C); 4 (D); 5 (AB); 6 (AC); 7 (AD); 8 (BC); 9 (BD); 10 (CD)
00272 % 11 (ABC); 12 (ABD); 13 (ACD); 14 (BCD); 15 (ABCD)
00273 case 1,
00274 dfterm_new = dfterm; tnames_new = tnames; msterm_new = msterm; intensity_new = intensity;
00275 msdenom = repmat(mse, size(msterm)); dfdenom = repmat(dfe, size(msterm));
00276 case 2,
00277 dfterm_new = dfterm; tnames_new = tnames; msterm_new = msterm; intensity_new = intensity;
00278 msdenom = [msterm(7), msterm(9), msterm(10), mse, msterm(12), msterm(13), mse, msterm(14), mse, mse, msterm(15), mse, mse, mse, mse];
00279 dfdenom = [dfterm(7), dfterm(9), dfterm(10), dfe, dfterm(12), dfterm(13), dfe, dfterm(14), dfe, dfe, dfterm(15), dfe, dfe, dfe, dfe];
00280 case 3, % only 11 terms in nesting case without AD, ABD, ACD, and ABCD: 1 (A); 2 (B); 3 (C); 4 (D); 5 (AB); 6 (AC);
00281 % 7 (BC); 8 (BD); 9 (CD); 10 (ABC); 11 (BCD);
00282 msterm_new = [msterm(1:6), msterm(8:11), msterm(14)]; % Throw out those four which do not exist for nesting: AD, ABD, ACD, and ABCD.
00283 intensity_new = [intensity(1:6), intensity(8:11), intensity(14)]; % Throw out those four which do not exist for nesting.
00284 dfterm_new = [dfterm(1:6)', dfterm(8:11)', dfterm(14)'];
00285 tnames_new = [tnames(1:6); tnames(8:11); tnames(14)]; % Only preserve those valid for nesting. Semicolon for coloumn catenation
00286 msdenom = [msterm(4), msterm(9), msterm(10), mse, msterm(9), msterm(10), msterm(14), mse, mse, msterm(14), mse];
00287 dfdenom = [dfterm(4), dfterm(9), dfterm(10), dfe, dfterm(9), dfterm(10), dfterm(14), dfe, dfe, dfterm(14), dfe];
00288 case 4,
00289 msterm_new = [msterm(1:6), msterm(8:11), msterm(14)]; % Throw out those four which do not exist for nesting: AD, ABD, ACD, and ABCD.
00290 intensity_new = [intensity(1:6), intensity(8:11), intensity(14)]; % Throw out those four which do not exist for nesting.
00291 dfterm_new = [dfterm(1:6)', dfterm(8:11)', dfterm(14)'];
00292 tnames_new = [tnames(1:6); tnames(8:11); tnames(14)]; % Only preserve those valid for nesting. Semicolon for coloumn catenation
00293 msdenom = [msterm(6), msterm(8), mse, msterm(10), msterm(11), mse, mse, msterm(14), mse, mse, mse]; % denominator MS
00294 dfdenom = [dfterm(6), dfterm(8), dfe, dfterm(10), dfterm(11), dfe, dfe, dfterm(14), dfe, dfe, dfe]; % denominator DF
00295 case 5, % only 9 terms in nesting case without AD, BD, ABD, ACD, BCD and ABCD: 1 (A); 2 (B); 3 (C); 4 (D); 5 (AB); 6 (AC); 7 (BC); 8 (CD); 9 (ABC)
00296 msterm_new = [msterm(1:6), msterm(8), msterm(10:11)]; % Throw out those four which do not exist for nesting: AD, BD, ABD, ACD, BCD and ABCD.
00297 intensity_new = [intensity(1:6), intensity(8), intensity(10:11)]; % Throw out those 6 which do not exist for nesting.
00298 dfterm_new = [dfterm(1:6)', dfterm(8), dfterm(10:11)'];
00299 tnames_new = [tnames(1:6); tnames(8); tnames(10:11)]; % Only preserve those valid for nesting. Semicolon for coloumn catenation
00300 msdenom = [msterm(4), msterm(4), msterm(10), mse, msterm(4), msterm(10), msterm(10), mse, msterm(10),0,0,0,0,0,0];
00301 % denominator MS: 6 extra 0's are for error-prone problem down below for contrasts
00302 dfdenom = [dfterm(4), dfterm(4), dfterm(10), dfe, dfterm(4), dfterm(10), dfterm(10), dfe, dfterm(10),0,0,0,0,0,0]; % denominator DF
00303 end % Close switch dsgn
00304
00305 case 5,
00306 switch dsgn %Order of the 2^5 - 1 = 31 terms: 1 (A); 2 (B); 3 (C); 4 (D); 5 (E); 6 (AB); 7 (AC); 8 (AD); 9 (AE);
00307 % 10 (BC); 11 (BD); 12 (BE); 13 (CD); 14 (CE); 15 (DE) 16 (ABC); 17 (ABD); 18 (ABE); 19 (ACD); 20 (ACE); 21 (ADE);
00308 % 22 (BCD); 23 (BCE); 24 (BDE); 25 (CDE); 26 (ABCD); 27 (ABCE); 28 (ABDE) 29 (ACDE); 30 (BCDE); 31 (ABCDE)
00309 case 1, % 31 terms
00310 dfterm_new = dfterm; tnames_new = tnames; msterm_new = msterm; intensity_new = intensity;
00311 msdenom = repmat(mse, size(msterm)); dfdenom = repmat(dfe, size(msterm));
00312 case 2, % 31 terms
00313 dfterm_new = dfterm; tnames_new = tnames; msterm_new = msterm; intensity_new = intensity;
00314 msdenom = [msterm(9), msterm(12), msterm(14), msterm(15), mse, msterm(18), msterm(20), msterm(21), mse, msterm(23), ...
00315 msterm(24), mse, msterm(25), mse, mse, msterm(27), msterm(28), mse, msterm(29), mse, mse, msterm(30), mse, mse, mse, ...
00316 msterm(31), mse, mse, mse, mse, mse];
00317 dfdenom = [dfterm(9), dfterm(12), dfterm(14), dfterm(15), dfe, dfterm(18), dfterm(20), dfterm(21), dfe, dfterm(23), ...
00318 dfterm(24), dfe, dfterm(25), dfe, dfe, dfterm(27), dfterm(28), dfe, dfterm(29), dfe, dfe, dfterm(30), dfe, dfe, dfe, ...
00319 dfterm(31), dfe, dfe, dfe, dfe, dfe];
00320 case 3, % 23 terms
00321 msterm_new = [msterm(1:8), msterm(10:17), msterm(19), msterm(22:26), msterm(30)];
00322 intensity_new = [intensity(1:8), intensity(10:17), intensity(19), intensity(22:26), intensity(30)];
00323 dfterm_new = [dfterm(1:8)', dfterm(10:17)', dfterm(19), dfterm(22:26)', dfterm(30)];
00324 tnames_new = [tnames(1:8); tnames(10:17); tnames(19); tnames(22:26); tnames(30)];
00325 msdenom = [msterm(5), msterm(12), msterm(14), msterm(15), mse, msterm(12), msterm(14), msterm(15), msterm(23), msterm(24), ...
00326 mse, msterm(25), mse, mse, msterm(23), msterm(24), msterm(25), msterm(30), mse, mse, mse, msterm(30), mse];
00327 dfdenom = [dfterm(5), dfterm(5), dfterm(5), dfterm(5), dfe, dfterm(12), dfterm(14), dfterm(15), dfterm(23), dfterm(24), ...
00328 dfe, dfterm(25), dfe, dfe, dfterm(23), dfterm(24), dfterm(25), dfterm(30),dfe, dfe, dfe, dfterm(30), dfe]; % denominator DF
00329
00330 case 4, % 23 terms
00331 msterm_new = [msterm(1:8), msterm(10:17), msterm(19), msterm(22:26), msterm(30)];
00332 intensity_new = [intensity(1:8), intensity(10:17), intensity(19), intensity(22:26), intensity(30)];
00333 dfterm_new = [dfterm(1:8)', dfterm(10:17)', dfterm(19), dfterm(22:26)', dfterm(30)];
00334 tnames_new = [tnames(1:8); tnames(10:17); tnames(19); tnames(22:26); tnames(30)];
00335 msdenom = [msterm(8), msterm(11), msterm(13), mse, msterm(15), msterm(17), msterm(19), mse, msterm(22), mse, ...
00336 msterm(24), mse, msterm(25), mse, msterm(26), mse, mse, mse, msterm(30), mse, mse, mse, mse];
00337 dfdenom = [dfterm(8), dfterm(11), dfterm(13), dfe, dfterm(15), dfterm(17), dfterm(19), dfe, dfterm(23), dfe, ...
00338 dfterm(24), dfe, dfterm(25), dfe, dfterm(26), dfe, dfe, dfe, dfterm(30), dfe, dfe, dfe, dfe]; % denominator DF
00339
00340 end
00341
00342 end % Close swtich NF
00343
00344
00345 % Compute an F statistic for each term
00346
00347 t = (msdenom>0); %Only calculate it when denominator is greater than 0
00348 fstat(t) = msterm_new(t) ./ msdenom(t);
00349
00350 %===================================================
00351 % Contrast tests
00352 if (Contr.do == 0),
00353 LC = 0; % assign this so that output argument LC would not be empty
00354 else
00355
00356
00357 if (NF == 1),
00358 % Get t values for contrast tests within each factor (1st order)
00359
00360 % dmat(:, num_col0+1:num_col1+num_col0) is the matrix for 1st order contrasts
00361 if (Contr.ord1.tot > 0),
00362 for (i = 1:1:Contr.ord1.tot),
00363 LC.t1(i).t = 0; % initializtion in case it is assigned later on due to denominator of 0
00364 LC.t1(i).value = Contr.ord1.cnt(i).vec * y0; % intensity for this 1st order contrast
00365 tmp = msdenom(Contr.ord1.cnt(i).idx1)*Contr.ord1.cnt(i).scalar;
00366 if (tmp > 0), LC.t1(i).t = LC.t1(i).value/sqrt(tmp); end
00367 end
00368 end
00369 end %if (NF == 1)
00370
00371 if (NF == 2),
00372 % Get t values for contrast tests within each factor (1st order)
00373
00374 % dmat(:, num_col0+1:num_col1+num_col0) is the matrix for 1st order contrasts
00375 if (Contr.ord1.tot > 0),
00376 for (i = 1:1:Contr.ord1.tot),
00377 LC.t1(i).t = 0; % initializtion in case it is assigned later on due to denominator of 0
00378 LC.t1(i).value = Contr.ord1.cnt(i).vec * y0; % intensity for this 1st order contrast
00379 tmp = msdenom(Contr.ord1.cnt(i).idx1)*Contr.ord1.cnt(i).scalar;
00380 if (tmp > 0), LC.t1(i).t = LC.t1(i).value/sqrt(tmp); end
00381 end
00382 end
00383 end %if (NF == 2)
00384
00385 if (NF == 3),
00386 % Get t values for contrast tests within each factor (1st order)
00387
00388 % dmat(:, num_col0+1:num_col1+num_col0) is the matrix for 1st order contrasts
00389 if (Contr.ord1.tot > 0),
00390 for (i = 1:1:Contr.ord1.tot),
00391 LC.t1(i).t = 0; % initializtion in case it is assigned later on due to denominator of 0
00392 LC.t1(i).value = Contr.ord1.cnt(i).vec * y0; % intensity for this 1st order contrast
00393 tmp = msdenom(Contr.ord1.cnt(i).idx1)*Contr.ord1.cnt(i).scalar;
00394 if (tmp > 0), LC.t1(i).t = LC.t1(i).value/sqrt(tmp); end
00395 end
00396 end
00397
00398 % Get t values for contrast tests (2nd order)
00399
00400 % dmat(:, num_col0+num_col1+1:num_col1+num_col0+num_col2) is the matrix for 2nd order contrasts
00401 if (Contr.ord2.tot > 0), % 7 terms: 1 (A); 2 (B); 3 (C); 4 (AB), 5 (AC), 6 BC, 7 ABC
00402 for (i = 1:1:Contr.ord2.tot),
00403 LC.t2(i).t = 0; % initializtion in case it is assigned later on due to denominator of 0
00404 switch Contr.ord2.cnt(i).idx1
00405 case 1,
00406 switch Contr.ord2.cnt(i).idx2
00407 case 2, what = msdenom(4); % MSAB
00408 case 3, what = msdenom(5) * (dsgn == 1 | dsgn == 2) + msdenom(3) * (dsgn == 4); % MSAC
00409 % For design type 4 -- BXC(A): the denominator for this contrast C(A) is MSBC(A), the one for main effect of C(A)!
00410 end
00411 case 2,
00412 if (Contr.ord2.cnt(i).idx2 == 3),
00413 what = msdenom(6) * (dsgn == 1 | dsgn == 2) + msdenom(5)* (dsgn == 3 | dsgn == 4); % MSBC
00414 else fprintf('\nSomething is wrong in the contrast coding!\n'); fprintf(2,'Halted: Ctrl+c to exit'); pause; end
00415 end %switch Contr.ord2.cnt(i).idx1
00416 LC.t2(i).value = Contr.ord2.cnt(i).vec * y0; % intensity for this 2nd order contrast
00417 tmp = what*Contr.ord2.cnt(i).scalar;
00418 if (tmp > 0), LC.t2(i).t = LC.t2(i).value/sqrt(tmp); end
00419 end %for (i = 1:1:Contr.ord2.tot),
00420 end %if (Contr.ord2.tot > 0),
00421
00422 % dmat(:, num_col0+num_col1+1:num_col1+num_col0+num_col2) is the matrix for 2nd order contrasts
00423 if (Contr.ord3.tot > 0),
00424 for (i = 1:1:Contr.ord3.tot),
00425 LC.t3(i).t = 0; % initializtion in case it is assigned later on due to denominator of 0
00426 switch Contr.ord3.cnt(i).idx1
00427 case 1,
00428 switch Contr.ord3.cnt(i).idx2
00429 case 2,
00430 if (Contr.ord3.cnt(i).idx3 == 3),
00431 what = msdenom(7) * (dsgn == 1 | dsgn == 2); % MSABC
00432 else fprintf('\nSomething is wrong in the contrast coding!\n'); fprintf(2,'Halted: Ctrl+c to exit'); pause; end
00433 case 3, fprintf('\nSomething is wrong in the contrast coding!\n'); fprintf(2,'Halted: Ctrl+c to exit'); pause;
00434 end
00435 case 2, fprintf('\nSomething is wrong in the contrast coding!\n'); fprintf(2,'Halted: Ctrl+c to exit'); pause;
00436 case 3, % Less likely occur
00437 fprintf('\nSomething is wrong in the contrast coding!\n'); fprintf(2,'Halted: Ctrl+c to exit'); pause;
00438 end
00439 LC.t3(i).value = Contr.ord3.cnt(i).vec * y0; % intensity for this 2nd order contrast
00440 tmp = what*Contr.ord3.cnt(i).scalar;
00441 if (tmp > 0), LC.t3(i).t = LC.t3(i).value/sqrt(tmp); end
00442 end %for (i = 1:1:Contr.ord3.tot),
00443 end
00444 end % If (NF == 3)
00445
00446 if (NF == 4),
00447 % Get t values for contrast tests within each factor (1st order)
00448
00449 % dmat(:, num_col0+1:num_col1+num_col0) is the matrix for 1st order contrasts
00450 if (Contr.ord1.tot > 0),
00451 for (i = 1:1:Contr.ord1.tot),
00452 LC.t1(i).t = 0; % initializtion in case it is assigned later on due to denominator of 0
00453 LC.t1(i).value = Contr.ord1.cnt(i).vec * y0; % intensity for this 1st order contrast
00454 tmp = msdenom(Contr.ord1.cnt(i).idx1)*Contr.ord1.cnt(i).scalar;
00455 if (tmp > 0), LC.t1(i).t = LC.t1(i).value/sqrt(tmp); end
00456 end
00457 end
00458
00459 % Get t values for contrast tests within each factor (2nd order)
00460 % dmat(:, num_col0+num_col1+1:num_col1+num_col0+num_col2) is the matrix for 2nd order contrasts
00461 if (Contr.ord2.tot > 0),
00462 for (i = 1:1:Contr.ord2.tot),
00463 LC.t2(i).t = 0; % initializtion in case it is assigned later on due to denominator of 0
00464 % switch dsgn
00465 switch Contr.ord2.cnt(i).idx1
00466 case 1,
00467 switch Contr.ord2.cnt(i).idx2
00468 case 2, what = msdenom(5); % MSAB
00469 case 3, what = msdenom(6); % MSAC
00470 case 4, what = msdenom(7) * (dsgn == 1 | dsgn == 2); % MSAD
00471 end
00472 case 2,
00473 switch Contr.ord2.cnt(i).idx2
00474 case 3, what = msdenom(8) * (dsgn == 1 | dsgn == 2) + msdenom(7)* (dsgn == 3 | dsgn == 4 | dsgn == 5); % MSBC
00475 case 4, what = msdenom(9) * (dsgn == 1 | dsgn == 2) + msdenom(8)* (dsgn == 3 | dsgn == 4); % Less likely to occur: MSBD
00476 end
00477 case 3, % Less likely occur
00478 switch Contr.ord2.cnt(i).idx2
00479 case 4, what = msdenom(10)* (dsgn == 1 | dsgn == 2) + msdenom(9)* (dsgn == 3 | dsgn == 4) + msdenom(7)*(dsgn == 5); % Less likely occur: MSCD
00480 end
00481 end
00482 % end
00483 % Can i remove the following loop by doing it in a matrix operation fashion?
00484 %for (j = 1:1:Contr2.cnt(i).NT), % each term of this contrast
00485 % LC2(i).value = LC2(i).value + Contr2.cnt(i).coef(j) * dmat(:, Contr2.cnt(i).code(j).pos)' * y0;
00486 %scalar = scalar + Contr2.cnt(i).coef(j) * Contr2.cnt(i).coef(j);
00487 %end
00488
00489 LC.t2(i).value = Contr.ord2.cnt(i).vec * y0; % intensity for this 2nd order contrast
00490 tmp = what*Contr.ord2.cnt(i).scalar;
00491 if (tmp > 0), LC.t2(i).t = LC.t2(i).value/sqrt(tmp); end
00492 end
00493 end
00494
00495 % dmat(:, num_col0+num_col1+1:num_col1+num_col0+num_col2) is the matrix for 2nd order contrasts
00496 if (Contr.ord3.tot > 0),
00497 for (i = 1:1:Contr.ord3.tot),
00498 LC.t3(i).t = 0; % initializtion in case it is assigned later on due to denominator of 0
00499 % if (dsgn == 3),
00500 switch Contr.ord3.cnt(i).idx1
00501
00502 case 1,
00503 switch Contr.ord3.cnt(i).idx2
00504
00505 case 2,
00506 switch Contr.ord3.cnt(i).idx3
00507 case 3, what = msdenom(11) * (dsgn == 1 | dsgn == 2) + msdenom(10) * (dsgn == 3 | dsgn == 4) + msdenom(9) * (dsgn == 5); % MSABC
00508 case 4, what = msdenom(12) * (dsgn == 1 | dsgn == 2); % MSABD not exist for (dsgn == 3 | dsgn == 4 | dsgn == 5)
00509 end
00510
00511 case 3,
00512 if (Contr.ord3.cnt(i).idx3 == 4), what = msdenom(13) * (dsgn == 1 | dsgn == 2); % MSACD
00513 else fprintf('\nSomething is wrong in the contrast coding!\n'); fprintf(2,'Halted: Ctrl+c to exit'); pause; end
00514
00515 case 4, fprintf('\nSomething is wrong in the contrast coding!\n'); fprintf(2,'Halted: Ctrl+c to exit'); pause;
00516 end
00517
00518 case 2,
00519 switch Contr.ord3.cnt(i).idx2
00520 case 3,
00521 if (Contr.ord3.cnt(i).idx3 == 4), what = msdenom(14) * (dsgn == 1 | dsgn == 2) + msdenom(11) * (dsgn == 3 | dsgn == 4); % MSBCD
00522 else fprintf('\nSomething is wrong in the contrast coding!\n'); fprintf(2,'Halted: Ctrl+c to exit'); pause;end
00523 case 4,
00524 fprintf('\nSomething is wrong in the contrast coding!\n'); fprintf(2,'Halted: Ctrl+c to exit'); pause;
00525 end
00526 case 3, % Less likely occur
00527 fprintf('\nSomething is wrong in the contrast coding!\n');
00528 fprintf(2,'Halted: Ctrl+c to exit'); pause;
00529 case 4,
00530 fprintf('\nSomething is wrong in the contrast coding!\n');
00531 fprintf(2,'Halted: Ctrl+c to exit'); pause;
00532 end
00533 % end
00534 LC.t3(i).value = Contr.ord3.cnt(i).vec * y0; % intensity for this 3rd order contrast
00535 tmp = what*Contr.ord3.cnt(i).scalar;
00536 if (tmp > 0), LC.t3(i).t = LC.t3(i).value/sqrt(tmp); end
00537 end
00538 end
00539 end % if (NF == 4)
00540
00541
00542 if (NF == 5),
00543 % Get t values for contrast tests within each factor (1st order)
00544
00545 % dmat(:, num_col0+1:num_col1+num_col0) is the matrix for 1st order contrasts
00546 if (Contr.ord1.tot > 0),
00547 for (i = 1:1:Contr.ord1.tot),
00548 LC.t1(i).t = 0; % initializtion in case it is assigned later on due to denominator of 0
00549 LC.t1(i).value = Contr.ord1.cnt(i).vec * y0; % intensity for this 1st order contrast
00550 tmp = msdenom(Contr.ord1.cnt(i).idx1)*Contr.ord1.cnt(i).scalar;
00551 if (tmp > 0), LC.t1(i).t = LC.t1(i).value/sqrt(tmp); end
00552 end
00553 end
00554
00555 % Get t values for contrast tests within each factor (2nd order)
00556 % dmat(:, num_col0+num_col1+1:num_col1+num_col0+num_col2) is the matrix for 2nd order contrasts
00557 %Order of the 2^5 - 1 = 31 terms: 1 (A); 2 (B); 3 (C); 4 (D); 5 (E); 6 (AB); 7 (AC); 8 (AD); 9 (AE);
00558 % 10 (BC); 11 (BD); 12 (BE); 13 (CD); 14 (CE); 15 (DE) 16 (ABC); 17 (ABD); 18 (ABE); 19 (ACD); 20 (ACE); 21 (ADE);
00559 % 22 (BCD); 23 (BCE); 24 (BDE); 25 (CDE); 26 (ABCD); 27 (ABCE); 28 (ABDE) 29 (ACDE); 30 (BCDE); 31 (ABCDE)
00560
00561 if (Contr.ord2.tot > 0),
00562 for (i = 1:1:Contr.ord2.tot),
00563 LC.t2(i).t = 0; % initializtion in case it is assigned later on due to denominator of 0
00564 % switch dsgn
00565 switch Contr.ord2.cnt(i).idx1
00566 case 1,
00567 switch Contr.ord2.cnt(i).idx2
00568 case 2, what = msdenom(6); % denominator for MSAB
00569 case 3, what = msdenom(7); % MSAC
00570 case 4, what = msdenom(8); % MSAD
00571 case 5, what = msdenom(9) * (dsgn == 1 | dsgn == 2); % MSAE
00572 end
00573 case 2,
00574 switch Contr.ord2.cnt(i).idx2
00575 case 3, what = msdenom(10) * (dsgn == 1 | dsgn == 2) + msdenom(9) * (dsgn == 3 | dsgn == 4); % MSBC
00576 case 4, what = msdenom(11) * (dsgn == 1 | dsgn == 2) + msdenom(10) * (dsgn == 3 | dsgn == 4); % Less likely to occur: MSBD
00577 case 5, what = msdenom(12) * (dsgn == 1 | dsgn == 2) + msdenom(11) * (dsgn == 3 | dsgn == 4); % Less likely to occur: MSBE
00578 end
00579 case 3, % Less likely occur
00580 switch Contr.ord2.cnt(i).idx2
00581 case 4, what = msdenom(13)* (dsgn == 1 | dsgn == 2) + msdenom(12) * (dsgn == 3 | dsgn == 4); % Less likely occur: MSCD
00582 case 5, what = msdenom(14)* (dsgn == 1 | dsgn == 2) + msdenom(13) * (dsgn == 3 | dsgn == 4); % Less likely occur: MSCE
00583 end
00584 case 4, % Less likely occur
00585 switch Contr.ord2.cnt(i).idx2
00586 case 5, what = msdenom(15)* (dsgn == 1 | dsgn == 2) + msdenom(14) * (dsgn == 3 | dsgn == 4); % Less likely occur: MSDE
00587 end
00588
00589 end
00590 % end
00591 % Can i remove the following loop by doing it in a matrix operation fashion?
00592 %for (j = 1:1:Contr2.cnt(i).NT), % each term of this contrast
00593 % LC2(i).value = LC2(i).value + Contr2.cnt(i).coef(j) * dmat(:, Contr2.cnt(i).code(j).pos)' * y0;
00594 %scalar = scalar + Contr2.cnt(i).coef(j) * Contr2.cnt(i).coef(j);
00595 %end
00596
00597 LC.t2(i).value = Contr.ord2.cnt(i).vec * y0; % intensity for this 2nd order contrast
00598 tmp = what*Contr.ord2.cnt(i).scalar;
00599 if (tmp > 0), LC.t2(i).t = LC.t2(i).value/sqrt(tmp); end
00600 end
00601 end
00602
00603 % dmat(:, num_col0+num_col1+1:num_col1+num_col0+num_col2) is the matrix for 2nd order contrasts
00604 if (Contr.ord3.tot > 0),
00605 for (i = 1:1:Contr.ord3.tot),
00606 LC.t3(i).t = 0; % initializtion in case it is assigned later on due to denominator of 0
00607 % if (dsgn == 3),
00608 switch Contr.ord3.cnt(i).idx1
00609
00610 case 1,
00611 switch Contr.ord3.cnt(i).idx2
00612
00613 case 2,
00614 switch Contr.ord3.cnt(i).idx3
00615 case 3, what = msdenom(16) * (dsgn == 1 | dsgn == 2) + msdenom(15) * (dsgn == 3 | dsgn == 4); % MSABC
00616 case 4, what = msdenom(17) * (dsgn == 1 | dsgn == 2) + msdenom(16) * (dsgn == 3 | dsgn == 4); % MSABD
00617 case 5, what = msdenom(18) * (dsgn == 1 | dsgn == 2); % MSABE
00618 end % switch Contr.ord3.cnt(i).idx3
00619
00620 case 3,
00621 switch Contr.ord3.cnt(i).idx3
00622 case 4, what = msdenom(19) * (dsgn == 1 | dsgn == 2) + msdenom(17) * (dsgn == 3 | dsgn == 4); % MSACD
00623 case 5, what = msdenom(20) * (dsgn == 1 | dsgn == 2); % MSACE
00624 end % switch Contr.ord3.cnt(i).idx3
00625
00626 case 4,
00627 if (Contr.ord3.cnt(i).idx3 == 5), what = msdenom(21) * (dsgn == 1 | dsgn == 2); % ADE
00628 else fprintf('\nSomething is wrong in the contrast coding!\n'); fprintf(2,'Halted: Ctrl+c to exit'); pause;
00629 end
00630 end % switch Contr.ord3.cnt(i).idx2
00631
00632 case 2,
00633 switch Contr.ord3.cnt(i).idx2
00634 case 3,
00635 switch Contr.ord3.cnt(i).idx3
00636 case 4, what = msdenom(22) * (dsgn == 1 | dsgn == 2) + msdenom(18) * (dsgn == 3 | dsgn == 4); % MSBCD
00637 case 5, what = msdenom(23) * (dsgn == 1 | dsgn == 2) + msdenom(19) * (dsgn == 3 | dsgn == 4); % MSBCE
00638 end % switch Contr.ord3.cnt(i).idx3
00639 case 4,
00640 if (Contr.ord3.cnt(i).idx3 == 5),
00641 if (dsgn == 1 | dsgn == 2),
00642 what = msdenom(24);
00643 elseif (dsgn == 3)
00644 what = msdenom(20);
00645 end
00646 % what = msdenom(24) * (dsgn == 1 | dsgn == 2) + msdenom(20) * (dsgn == 3); % BDE
00647 else fprintf('\nSomething is wrong in the contrast coding!\n'); fprintf(2,'Halted: Ctrl+c to exit'); pause; end
00648 end % switch Contr.ord3.cnt(i).idx2
00649 case 3,
00650 if (Contr.ord3.cnt(i).idx2 == 4 & Contr.ord3.cnt(i).idx3 == 5),
00651 if (dsgn == 1 | dsgn == 2),
00652 what = msdenom(25);
00653 elseif (dsgn == 3),
00654 what = msdenom(21);
00655 end
00656 % what = msdenom(25) * (dsgn == 1 | dsgn == 2) + msdenom(21) * (dsgn == 3); % MSCDE
00657 else
00658 fprintf('\nSomething is wrong in the contrast coding!\n');
00659 fprintf(2,'Halted: Ctrl+c to exit'); pause;
00660 end
00661 case 4,
00662 fprintf('\nSomething is wrong in the contrast coding!\n');
00663 fprintf(2,'Halted: Ctrl+c to exit'); pause;
00664 end % switch Contr.ord3.cnt(i).idx1
00665 % end
00666 LC.t3(i).value = Contr.ord3.cnt(i).vec * y0; % intensity for this 3rd order contrast
00667 tmp = what*Contr.ord3.cnt(i).scalar;
00668 if (tmp > 0), LC.t3(i).t = LC.t3(i).value/sqrt(tmp); end
00669 end % for (i = 1:1:Contr.ord3.tot)
00670 end % if (Contr.ord3.tot > 0)
00671
00672 if (Contr.ord4.tot > 0),
00673 for (i = 1:1:Contr.ord4.tot),
00674 LC.t4(i).t = 0; % initializtion in case it is assigned later on due to denominator of 0
00675 switch Contr.ord4.cnt(i).idx1
00676
00677 case 1,
00678 switch Contr.ord4.cnt(i).idx2
00679 case 2,
00680 switch Contr.ord4.cnt(i).idx3
00681 case 3,
00682 switch Contr.ord4.cnt(i).idx4
00683 case 4,
00684 if (dsgn == 1 | dsgn == 2),
00685 what = msdenom(26);
00686 elseif (dsgn == 3),
00687 what = msdenom(22);
00688 end
00689 % what = msdenom(26) * (dsgn == 1 | dsgn == 2) + msdenom(22) * (dsgn == 3); % MSABCD
00690 case 5,
00691 if (dsgn == 1 | dsgn == 2),
00692 what = msdenom(27);
00693 end
00694 % what = msdenom(27) * (dsgn == 1 | dsgn == 2); % MSABCE
00695 end % switch Contr.ord4.cnt(i).idx4
00696 case 4,
00697 if (Contr.ord4.cnt(i).idx4 == 5),
00698 if (dsgn == 1 | dsgn == 2),
00699 what = msdenom(28);
00700 end
00701 % what = msdenom(28) * (dsgn == 1 | dsgn == 2); % ABDE
00702 else fprintf('\nSomething is wrong in the contrast coding!\n'); fprintf(2,'Halted: Ctrl+c to exit'); pause;
00703 end
00704 case 5, fprintf('\nSomething is wrong in the contrast coding!\n');
00705 fprintf(2,'Halted: Ctrl+c to exit'); pause;
00706 end % switch Contr.ord4.cnt(i).idx3
00707
00708 case 3,
00709 switch Contr.ord4.cnt(i).idx3
00710 case 4,
00711 if (Contr.ord4.cnt(i).idx4 == 5),
00712 if (dsgn == 1 | dsgn == 2),
00713 what = msdenom(29);
00714 end
00715 % what = msdenom(29) * (dsgn == 1 | dsgn == 2); % ACDE
00716 else fprintf('\nSomething is wrong in the contrast coding!\n'); fprintf(2,'Halted: Ctrl+c to exit'); pause;
00717 end
00718 case 5, fprintf('\nSomething is wrong in the contrast coding!\n');
00719 fprintf(2,'Halted: Ctrl+c to exit'); pause;
00720 end % switch Contr.ord4.cnt(i).idx3
00721
00722 case 4,
00723 fprintf('\nSomething is wrong in the contrast coding!\n');
00724 fprintf(2,'Halted: Ctrl+c to exit'); pause;
00725 end % switch Contr.ord4.cnt(i).idx2
00726
00727 case 2,
00728 switch Contr.ord4.cnt(i).idx2
00729 case 3,
00730 switch Contr.ord4.cnt(i).idx3
00731 case 4,
00732 if (Contr.ord4.cnt(i).idx4 == 5),
00733 if (dsgn == 1 | dsgn == 2),
00734 what = msdenom(30);
00735 end
00736 % what = msdenom(30) * (dsgn == 1 | dsgn == 2) + msdenom(23) * (dsgn == 3); % BCDE
00737 else fprintf('\nSomething is wrong in the contrast coding!\n'); fprintf(2,'Halted: Ctrl+c to exit'); pause;
00738 end
00739 case 5, fprintf('\nSomething is wrong in the contrast coding!\n');
00740 fprintf(2,'Halted: Ctrl+c to exit'); pause;
00741 end % switch Contr.ord4.cnt(i).idx3
00742 case 4,
00743 fprintf('\nSomething is wrong in the contrast coding!\n'); fprintf(2,'Halted: Ctrl+c to exit'); pause;
00744 end % switch Contr.ord4.cnt(i).idx2
00745 case {3,4,5,}
00746 fprintf('\nSomething is wrong in the contrast coding!\n');
00747 fprintf(2,'Halted: Ctrl+c to exit'); pause;
00748 end % switch Contr.ord4.cnt(i).idx1
00749 LC.t4(i).value = Contr.ord4.cnt(i).vec * y0; % intensity for this 3rd order contrast
00750 tmp = what*Contr.ord4.cnt(i).scalar;
00751 if (tmp > 0), LC.t4(i).t = LC.t4(i).value/sqrt(tmp); end
00752 end % for (i = 1:1:Contr.ord4.tot)
00753 end % if (Contr.ord4.tot > 0)
00754
00755 end % if (NF == 5)
00756
00757
00758 end %if (Contr.do == 0)
00759
00760 err = 0;
00761 return;
