00001 function [err, Qd, s, termname, nterms, sindices, dfbothSS, modw, modwo, tnames, dfterm, dfe, Contr] = PreProc(n, NF, group, varnames, FL, Contr, cov, unbalanced)
00002 %
00003 % [err, Qd, Rd, sindices, dfboth, modw, modwo, dfx] = PreProc(n,group,varnames)
00004 %
00005 %Purpose:
00006 %
00007 % Returns QR decomposition results for the design matrix projected to null space
00008 %
00009 %Input Parameters:
00010 % n: number of datasets = total # of combinations including repeats
00011 % group: cell array of factor levels
00012 % varnames:
00013 %
00014 %Output Parameters:
00015 % err : 0 No Problem
00016 % : 1 Problems
00017 %
00018 %
00019 %
00020 %Key Terms:
00021 %
00022 %More Info :
00023 %
00024 %
00025 %
00026 %
00027 % Author : Gang Chen
00028 % Date : Tue Mar 23 13:57:52 EST 2004
00029 % SSCC/NIMH/ NIH, Bethesda MD 20892
00030
00031
00032 %Define the function name for easy referencing
00033 FuncName = 'PreProc.m';
00034
00035 %Debug Flag
00036 DBG = 1;
00037
00038 %initailize return variables
00039 err = 1;
00040
00041 % Don't worry about NaN's at this point.
00042 group = group(:); % what's it for?
00043 ng = length(group); % number of factors
00044 termlist = makemodel(ng, ng); % Generate terms for all main effects plus various interactions
00045
00046 for j=1:ng
00047 gj = group{j};
00048 if (size(gj,1) == 1), gj = gj(:); end
00049 if (size(gj,1) ~= n)
00050 error('Factor %d must have %d elements.',j,n);
00051 end
00052 if (ischar(gj)), gj = cellstr(gj); end
00053 group{j} = gj;
00054 end
00055
00056 gdum = cell(ng,1);
00057 dfvar = zeros(ng,1);
00058
00059 if (unbalanced.yes == 0), % Balanced designs
00060 vconstr = cell(ng,1);
00061 %vmean = cell(ng,1); %vmean is never used!!!
00062 end
00063
00064 for j=1:ng % for each factor
00065 gj = group{j};
00066 [gij,gnj] = grp2idx(gj); % Create index vector from a grouping variable: gij is a vector
00067 % taking integer values from 1 up to the number of unique entries in gj
00068 % gnj is a cell array of names, so that gnj(gij) reproduces gj
00069 nlevels = size(gnj,1); % levels for this factor
00070 dfvar(j) = nlevels - 1; % D. F. for this factor
00071
00072 if (unbalanced.yes == 0), % balanced
00073 if (cov.do & j==cov.marker)
00074 gdum{j} = gj;
00075 dfvar(j) = 1; % D. F. = 1
00076 vconstr{j} = zeros(0,1);
00077 % vmean{j} = 1;
00078 else
00079 gdum{j} = idummy(gij, 3);
00080 vconstr{j} = idummy(1:nlevels)';
00081 % vmean{j} = ones(1,nlevels) / nlevels; % array (1Xnlevels) of one ones, but vmean is never used in the code!!!!!!!!!1
00082 end
00083 else % Unbalanced designs
00084 if (cov.do & j==cov.marker)
00085 gdum{j} = gj;
00086 dfvar(j) = 1; % D. F. = 1
00087 vconstr{j} = zeros(0,1);
00088 % vmean{j} = 1;
00089 else
00090 gdum{j} = idummy(gij, 3);
00091 end
00092 end % if (unbalanced.yes == 0): Only for balanced designs
00093
00094 end
00095
00096 % Create dummy variable arrays for each term in the model.
00097 nterms = size(termlist,1); % Number of rows (1st dimension) in termlist
00098 [sterms,sindex] = sortrows(termlist); % Sort terms in ascending order.
00099 ncols = 1;
00100 nconstr = 0;
00101
00102 termdum = cell(nterms, 1); % cell array of dummy variables which are design matrix cols
00103 termconstr = cell(nterms,1); % constraints to make each term well defined
00104 levelcodes = cell(nterms, 1); % codes for levels of each M row
00105 tnames = cell(nterms, 1); % name for entire term, e.g. A*B
00106 dfterm0 = zeros(nterms, 1); % nominal d.f. for each term
00107 termvars = cell(nterms, 1); % list of vars in each term
00108 termlength = zeros(size(sindex)); % length of each term (number of columns)
00109
00110 %randomterm = find(termlist*randomvar > 0); % indices of random terms
00111
00112 % For each term,
00113 for j=1:nterms
00114 % Loop over elements of the term
00115 df0 = 1;
00116 tm = sterms(j,:);
00117 tdum = []; % empty term so far
00118 tconstr = 1; % empty constraints so far
00119 tn = ''; % blank name so far
00120 vars = find(tm); % Find indices of nonzero elements
00121 for varidx = 1:length(vars)
00122 % Process each variable participating in this term
00123 varnum = vars(varidx); % variable name
00124 tm(varnum) = 0; % term without this variable
00125 df0 = df0 * dfvar(varnum); % d.f. so far
00126
00127 % Combine its dummy variable with the part computed so far
00128 G = gdum{varnum}; % dummy vars for this grouping var
00129 nlevterm = size(tdum,2); % levels for term so far
00130 nlevgrp = size(G,2); % levels for this grouping var
00131 tdum = termcross(G,tdum); % combine G into term dummy vars
00132
00133 % Construct the term name and constraints matrix
00134 % if (ismember(varnum, NF) & cov.do), % for the covariate, which is the last factor
00135 if (ismember(varnum, cov.marker) & cov.do), % for the covariate, which is the last factor
00136 vconstr = ones(0,1);
00137 else
00138 vconstr = ones(1, nlevgrp);
00139 end
00140 if (isempty(tn))
00141 tn = varnames{varnum};
00142 tconstr = vconstr;
00143 else
00144 tn = [varnames{varnum} '*' tn];
00145 tconstr = [kron(vconstr,eye(size(tconstr,2)));
00146 kron(eye(length(vconstr)),tconstr)]; % Kronecker
00147 end
00148
00149 % If the rest of this term is computed, take advantage of that
00150 prevterms = sterms(1:j-1,:);
00151 oldtm = find((prevterms * tm') == sum(tm)); % same vars in old term
00152 oldtm = oldtm((prevterms(oldtm,:) * ~tm') == 0); % and no others
00153 if (length(oldtm) > 0)
00154 k = sindex(oldtm(1));
00155 tdum = termcross(termdum{k}, tdum);
00156 oconstr = termconstr{k};
00157 tconstr = [kron(tconstr, eye(size(oconstr,2)));
00158 kron(eye(size(tconstr,2)), oconstr)];
00159 tn = [tn '*' tnames{k}];
00160 df0 = df0 * dfterm0(k);
00161 break;
00162 end
00163 end
00164
00165 % Store this term's dummy variables and name
00166 k = size(tdum, 2);
00167 termlength(sindex(j),1) = k;
00168 ncols = ncols + k;
00169 sj = sindex(j);
00170 termdum{sj} = tdum;
00171 termconstr{sj} = tconstr;
00172 termvars{sj} = vars;
00173 levelcodes{sj} = fullfact(dfvar(vars)+1);
00174 if (isempty(tn)), tn = 'Constant'; end
00175 tnames{sj,1} = tn;
00176 dfterm0(sj) = df0;
00177 nconstr = nconstr + size(tconstr,1);
00178 end
00179 tnames{length(tnames)+1,1} = 'Error';
00180
00181 % Create the full design matrix
00182 dmat = ones(n, ncols); % to hold design matrix of nXncols
00183 cmat = zeros(nconstr,ncols); % to hold constraints matrix
00184 cbase = 0; % base from which to fill in cmat
00185 termname = zeros(ncols,1);
00186 termstart = cumsum([2; termlength(1:end-1)]);
00187 termend = termstart + termlength - 1;
00188 for j=1:nterms
00189 clist = termstart(j):termend(j);
00190 dmat(:, clist) = termdum{j};
00191 C = termconstr{j};
00192 nC = size(C,1);
00193 cmat(cbase+1:cbase+nC,clist) = C;
00194 termname(clist) = j;
00195 cbase = cbase + nC;
00196 end
00197
00198 [err, Qd, dfx, dmat2] = QRDecom(dmat, cmat);
00199 %dfe = n - dfx;
00200
00201 % Determine which models to compare for testing each term
00202 ssw = -ones(nterms, 1); % sum of squares with this term
00203 sswo = ssw; % sum of squares without this term
00204 dfw = ssw; % residual d.f. with this term
00205 dfwo = ssw; % residual d.f. without this term
00206
00207 if (unbalanced.yes == 1)
00208 modw = tril(ones(nterms)); % get model with this term
00209 k = nterms; % locations of model with all terms
00210 else
00211
00212 % Only apply type III for sums of squares
00213 modw = ones(nterms);
00214 %k = 1:nterms;
00215 %TnotC = termsnotcontained(termlist);
00216
00217 end % if (unbalanced.yes == 1)
00218
00219 modw = logical(modw); % get model with this term
00220 modwo = logical(modw - eye(nterms)); % get model without this term
00221
00222 dfw(1:nterms) = dfx;
00223
00224 % Fit each model separately
00225 dfboth = [dfw; dfwo];
00226 dfbothSS = dfboth; % For usage in ss.m
00227
00228 % Consider interactions before their components for type 3 ss, so
00229 % examine the terms in decreasing order of the number factors in the term
00230
00231
00232 if (unbalanced.yes == 1)
00233 sindices = [(1:size(termlist,1)), (1:size(termlist,1))]';
00234 else
00235 termsize = sum(termlist,2)'; %sum of all the elements along 2nd dimension (row)
00236 [stermsize,sindices] = sort(termsize); % sort in ascending order, output is stored in stermsize; sindices is the index for the new array
00237 sindices = [sindices(:); sindices(:)];
00238
00239 end % if (unbalanced.yes == 1)
00240
00241 % Here QR decomposition is done, which is voxel independent
00242
00243 s(length(sindices)).Qdt = [];
00244
00245 for j=length(sindices):-1:1
00246 % Find the next model index to fit
00247 k = sindices(j);
00248
00249 % Look in unsorted arrays to see if we have already fit this model
00250 if j>nterms
00251 k0 = k+nterms;
00252 else
00253 k0 = k;
00254 end
00255 if dfboth(k0)~=-1
00256 continue
00257 end
00258
00259 % Find the model with this index
00260 if (j > nterms)
00261 thismod = modwo(k, :);
00262 else
00263 thismod = modw(k, :);
00264 end
00265
00266 % Get the design matrix for this model
00267 keepterms = find(thismod);
00268 clist = ismember(termname, [0 keepterms]);
00269 X = dmat2(:,clist);
00270 C = cmat(:,clist);
00271
00272 % Fit this term
00273 [err, s(j).Qdt, dfx0] = QRDecom(X, C);
00274
00275 % Use these results for each term that requires them
00276
00277 mod0 = repmat(thismod, nterms, 1);
00278 k = find(all(modw == mod0, 2));
00279 dfw(k) = dfx0;
00280 dfboth(k) = 0;
00281
00282 k = find(all(modwo == mod0, 2));
00283 dfwo(k) = dfx0;
00284 dfboth(nterms+k) = 0;
00285 end
00286 clear mod0
00287
00288 dfterm = dfw - dfwo;
00289 dfe = n-dfx; %residual degrees of freedom
00290
00291 if (Contr.do == 1),
00292
00293 % In design matrix dmat, the first column is all one's, for the total mean. Then there are totally
00294 % FL(1).N_level + FL(2).N_level + FL(3).N_level + FL(4).N_level columns for the main effects.
00295 % Next 2nd order interactions, 3rd order interaction, and 4th order interactions.
00296
00297 % Store only those mean columns in design matrix.
00298
00299 %num_col0 = 1; % the 1st column is for grand mean (0 order)
00300
00301 num_col(1) = 0;
00302 for (i = 1:1:NF),
00303 num_col(1) = num_col(1) + FL(i).N_level; %
00304 end
00305
00306 % Ignore 1st col since it is total mean
00307 %dmat_mean = dmat(:, (num_col0 + 1):(num_col0 + num_col(1)));
00308
00309 % Get the number in the sum for each mean, which happens to be in the diagonal of X'X.
00310 % This can also be otained through the user input variables, but it is generic with
00311 % the matrix operation, especially for unbalanced design.
00312 %sum_num = diag(dmat_mean' * dmat_mean);
00313
00314
00315 if (NF > 1),
00316 num_col(2) = 0;
00317 for (i = 1:1:(NF-1)),
00318 for (j = (i+1):1:NF),
00319 num_col(2) = num_col(2) + FL(i).N_level*FL(j).N_level; %Columns for 2nd order interactions
00320 end
00321 end
00322 end
00323
00324 if (NF > 2),
00325 num_col(3) = 0;
00326 for (i = 1:1:(NF-2)),
00327 for (j = (i+1):1:(NF-1)),
00328 for (k = (j+1):1:NF),
00329 num_col(3) = num_col(3) + FL(i).N_level*FL(j).N_level*FL(k).N_level; %Columns for 3rd order interactions
00330 end
00331 end
00332 end
00333 end
00334
00335 %if (NF == 4),
00336 % num_col(4) = 1;
00337 % for (i = 1:1:NF),
00338 % num_col(4) = num_col(4)*FL(i).N_level; %Columns for 4th order interactions
00339 % end
00340 %end
00341
00342 % for every design
00343 if (Contr.ord1.tot > 0), % 1st order contrasts
00344 [err, Contr] = ContrVec(1, n, NF, group, dmat, Contr, FL, num_col);
00345 end % if (Contr1.tot > 0)
00346
00347 if (NF > 1 & Contr.ord2.tot > 0), % 2nd order contrasts
00348 [err, Contr] = ContrVec(2, n, NF, group, dmat, Contr, FL, num_col);
00349 end % if (Contr2.tot > 0)
00350
00351 if (NF > 2 & Contr.ord3.tot > 0), % 3rd order contrasts
00352 [err, Contr] = ContrVec(3, n, NF, group, dmat, Contr, FL, num_col);
00353 end % if (Contr3.tot > 0)
00354
00355 if (NF > 3 & Contr.ord4.tot > 0), % 4th order contrasts
00356 [err, Contr] = ContrVec(4, n, NF, group, dmat, Contr, FL, num_col);
00357 end % if (Contr3.tot > 0)
00358
00359 end % if (Contr.do == 1)
00360
00361 err = 0;
00362 return;
