Description of TrackingAnalysis

0001 % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0002 % Script for Performing Slice-Based
0003 % Multiple Regression on fMRI Time Series
0004 %
0005 % Written by:
0006 % John Darrell Van Horn, Ph.D., M.Eng.
0007 % Dartmouth Brain Imaging Center
0008 % Dartmouth College
0009 % Hanover, NH 03753
0010 %
0011 % Phone: (603) 646-2909
0012 % Email: John.D.Van.Horn@dartmouth.edu
0013 %
0014 %
0015 % Last modified: 11-29-2001
0016 %
0017 % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0018 
0019 % **********************
0020 % Clear Matlab Workspace
0021 clear;
0022 clc;
0023 
0024 % ****************
0025 % Close All Figure
0026 close all;
0027 
0028 % ***************************************
0029 % Define Endian Type for Image Data Input
0030 ENDIANIN='ieee-le';   % Use this when running on Intel-Based Processors and
0031                       % with data written using an Intel Processor
0032                       % e.g. A DELL Linux Box
0033 
0034 %ENDIANIN='ieee-be';  % Use this when running on SUN-Based Processors and
0035                       % with data written using a SUN Processor
0036                       % e.g. dexter.dartmouth.edu
0037 
0038 % ****************************************
0039 % Define Endian Type for Image Data Output
0040 ENDIANOUT='ieee-le';  % Use this when running on Intel-Based Processors and
0041                       % with data written using an Intel Processor
0042                       % e.g. A DELL Linux Box
0043 
0044 %ENDIANOUT='ieee-be'; % Use this when running on SUN-Based Processors and
0045                       % with data written using a SUN Processor
0046                       % e.g. dexter.dartmouth.edu
0047 
0048 % *************************
0049 % Start Matlab Elapsed Time
0050 tic;
0051 
0052 % *********************
0053 % Set Masking Threshold
0054 THRESH=550.0;
0055 
0056 % **************************************
0057 % Flag for Normalizing the Design Matrix
0058 NORMDESMTX=0;
0059 
0060 % ***************************
0061 % fMR Data Normalization Mode
0062 DATA_NORM_MODE = 3;
0063 
0064 % ***************************
0065 % Flag for Outlier Protection
0066 OUTLIER_CORRECTION=1;
0067 OUTLIER_THRESH=4.0;
0068 
0069 % ************************************************
0070 % Flag for The Existence of a Matlab Contrast File
0071 CONTRAST_FILE = 0;
0072 
0073 % ************************************
0074 % Critical Probability Threshold Value
0075 pcrit=0.05;
0076 
0077 % ***************************
0078 % Set A Small Tolerance Value
0079 TINY=1e-5;
0080 
0081 % ********************
0082 % F_indiv Scatterplots
0083 F_INDIV_PLOTS=0;
0084 
0085 % ***************
0086 % Local Max Plots
0087 LOCAL_MAX_PLOTS=1;
0088 
0089 % *********************
0090 % Image File Dimensions
0091 TR=2.0;
0092 RUNS=4;
0093 VOLS_PER_RUN=126;
0094 XDIM=64;
0095 YDIM=64;
0096 ZDIM=25;
0097 XVOXDIM=3.75;
0098 YVOXDIM=3.75;
0099 ZVOXDIM=5.5;
0100 N=RUNS*VOLS_PER_RUN;
0101 NUMVOX=XDIM*YDIM*ZDIM;
0102 
0103 % ************************************************************
0104 % If Christian Buchel's Structural Equation Modeling Routine's
0105 % are Available set SEM=1, else SEM=0
0106 SEM=1;
0107 
0108 % *************************
0109 % Select Design Matrix File
0110 [DesMtxFileName, CWD]=uigetfile('*.mat','Select Design Matrix Matlab File');
0111 
0112 % *********************
0113 % Load Design Matrix, X
0114 disp('Loading Design Matrix...');
0115 close;
0116 eval(['load ' [CWD DesMtxFileName] ';']);
0117 DesMtx=DesMtxB;  
0118 clear DesMtxA DesMtxB;
0119 D=size(DesMtx);  %For the rest of the script the Design Matrix is called DesMtx
0120 VARIABLES=D(1);
0121 TIMEPOINTS=D(2);
0122 if(length(D)==3)
0123   SLICES=D(3);
0124 elseif(length(D)==2)
0125   SLICES=1;
0126 else
0127   disp(['This ain`t no stinkin design matrix!']);
0128   exit;
0129 end
0130 
0131 % ********************************************************************
0132 % If Design Matrix is not 3D Slice-based (i.e. is only a single slice)
0133 % Replicate it to create a Design volume
0134 if(SLICES==1)
0135      tmp=DesMtx;
0136      DesMtx=zeros(VARIABLES,TIMEPOINTS,ZDIM);
0137      for i=1:ZDIM,
0138      DesMtx(:,:,i)=tmp(:,:);
0139      end
0140      D=size(DesMtx);
0141      VARIABLES=D(1);
0142      TIMEPOINTS=D(2);
0143      SLICES=ZDIM;
0144 end
0145 
0146 % *************************************************
0147 % Create Selection Vector and Define Masking Vector
0148 % for Removing Effects of Confounding Variables
0149 % e.g. those variables that will be regressed out
0150 % for the reduced model
0151 S=zeros(1,VARIABLES);
0152 S([37:end])=1;
0153 
0154 % **********************
0155 % Load a Contrast Matrix
0156 if(CONTRAST_FILE==1),
0157 
0158   % Select Contrast Matrix Matlab File
0159   [ContMtxFileName, CWD]=uigetfile('*.mat','Select Contrast Matrix Matlab File');
0160   eval(['if(exist(' [CWD ContMtxFileName] '))']);
0161   eval(['load ' [CWD ContMtxFileName] ';']);
0162 else
0163   disp(['No contrast matrix exists by that name.']);
0164   disp(['Carrying on anyway......']);
0165   C=[(1-eye(VARIABLES-sum(S),VARIABLES-sum(S))) ones(VARIABLES-sum(S),sum(S))];
0166 end
0167 
0168 % ***************************
0169 % Select Image List Text File
0170 [ImageFileList, DataDir]=uigetfile('*.txt');
0171 
0172 % ***********************
0173 % Read in Image Filenames
0174 disp('Reading in Filenames...');
0175 close;
0176 drawnow;
0177 fname=[DataDir ImageFileList];
0178 fid=fopen(fname,'r',ENDIANIN);
0179 
0180 for i=1:N,
0181   F=fgetl(fid);
0182   FILES(i,[1:length(F)])=F;
0183 end
0184 fclose(fid);
0185 
0186 % *************************************
0187 % Read in Functional Image Data Volumes
0188 disp('Loading Functional Image Data...');
0189 Y=zeros(N,NUMVOX);
0190 for i=1:N,
0191     fid=fopen(FILES(i,:),'r',ENDIANIN);
0192         U=fread(fid,NUMVOX,'int16');
0193         Y(i,:)=U';
0194         fclose(fid);
0195 end
0196 
0197 % *************************
0198 % Load or Create Image Mask
0199 if(exist('Mask.img'))
0200   disp('Loading mask image from disk...');
0201   fid=fopen('Mask.img','r',ENDIANIN);
0202   mask=fread(fid,NUMVOX,'int8');
0203   fclose(fid);
0204 else
0205   disp('Generating Image Mask...');
0206   I=(Y>THRESH);
0207   mask=prod(I);
0208   clear I;
0209 
0210   % Save Mask Image to PWD
0211   fid=fopen(['Mask.img'],'w',ENDIANOUT);
0212   fwrite(fid,mask,'int8');
0213   fclose(fid);
0214   str=['!makeheader4 Mask.hdr ', num2str(XDIM), ' ', num2str(YDIM), ' ', num2str(ZDIM),...
0215   ' ', num2str(XVOXDIM), ' ', num2str(YVOXDIM), ' ', num2str(ZVOXDIM), ' 1 8 1 0 1 0 2'];
0216   eval(str);
0217 end
0218 V=find(mask>0);
0219 
0220 if(length(V)>0)
0221   disp(['Analysing ' num2str(sum(length(V))) ' voxels above threshold...']);
0222 else
0223   disp(['No voxels found above threshold.']);
0224   clear;
0225   quit;
0226 end
0227 
0228 % ********************************
0229 % Standard Normalize Design Matrix
0230 if (NORMDESMTX==1)
0231    disp('Normalizing Design Matrix...');
0232    Z=stdnormalize(DesMtx);
0233 else
0234    Z=DesMtx;
0235 end
0236 
0237 % ******************************
0238 % Normalize fMR Image Data by...
0239 switch DATA_NORM_MODE
0240   case 1     % 1) Proportional Scaling
0241 
0242      disp('Normalizing Image Data by Proportional Scaling...');
0243      ImageMeans=zeros(1,N);
0244      ImageMeans=mean(Y(:,V)');
0245      for i=1:N,
0246     for j=1:length(V),
0247             Y(i,V(j))=Y(i,V(j))./ImageMeans(i);
0248         end
0249      end
0250 
0251   case 2     % 2) Standard Normalization Over All Scans
0252      disp('Normalizing Image Data by Standard Normalization...');
0253      ImageMeans=zeros(1,N);
0254      ImageMeans=mean(Y(:,V)');
0255      ImageStds=std(Y(:,V)');
0256      for i=1:N,
0257         for j=1:length(V),
0258             Y(i,V(j))=(Y(i,V(j))-ImageMeans(i))./ImageStds(i);
0259         end
0260      end
0261 
0262   case 3     % 3) Standard Normalization of Scans Within Run
0263      disp('Normalizing Image Data by Standard Normalization Within Run...');
0264 
0265      for i=1:RUNS,
0266      ImageMeans=zeros(1,N);
0267      ImageMeans=mean(Y([(i-1)*VOLS_PER_RUN+1:i*VOLS_PER_RUN],V)');
0268      ImageStds=std(Y([(i-1)*VOLS_PER_RUN+1:i*VOLS_PER_RUN],V)');
0269         for k=((i-1)*VOLS_PER_RUN+1):(i*VOLS_PER_RUN),
0270            for l=1:length(V),
0271             Y(k,V(l))=(Y(k,V(l))-ImageMeans(i))./ImageStds(i);
0272            end
0273         end
0274      end
0275 
0276   otherwise
0277      disp(['Unknown MR Data Normalization Method']);
0278      exit;
0279 end
0280 
0281      save Y;
0282 
0283 % **************************************
0284 % Perform Slice-Wise Pseudo-Inversion of
0285 % the Design Matrix for Use Below
0286 PinvDesMtx = zeros(VARIABLES,VARIABLES,SLICES);
0287 for i=1:SLICES,
0288     PinvDesMtx(:,:,i)=pinv(Z(:,:,i)*Z(:,:,i)');
0289 end
0290 
0291 
0292 % **************************************
0293 % Perform Slice-wise Multiple Regression
0294 disp('Performing Slice-wise Multiple Regression...');
0295 Beta=zeros(VARIABLES,NUMVOX);
0296 H=zeros(1,NUMVOX);
0297 E=zeros(1,NUMVOX);
0298 T=zeros(1,NUMVOX);
0299 Lambda=zeros(1,NUMVOX);
0300 Lambda_reduced=zeros(1,NUMVOX);
0301 F_omnibus=zeros(1,NUMVOX);
0302 F_reduced=zeros(1,NUMVOX);
0303 F_reduced_mask=zeros(1,NUMVOX);
0304 F_indiv=zeros(VARIABLES-sum(S),NUMVOX);
0305 
0306 if(OUTLIER_CORRECTION==1)
0307      disp(' ');
0308      disp('Outlier Correction in effect...');
0309      disp(' ');
0310 end
0311 
0312 zlast=0;
0313 for i=1:length(V),
0314 
0315   % When enabled, inspect y vector and set outliers equal to zero
0316   % and set values in the data equal to mean of the rest of the series
0317   if(OUTLIER_CORRECTION==1)  
0318     q=Y(:,V(i));
0319     q=(q-mean(q))./std(q);
0320     for j=1:N
0321       if(abs(q(j))>OUTLIER_THRESH)
0322          Y(j,V(i))=(sum(Y(:,V(i)))-Y(j,V(i)))./(N-1);
0323       end
0324     end
0325   end
0326 
0327   y=Y(:,V(i));
0328   y=y-mean(y);
0329 
0330   z=floor(V(i)./(XDIM*YDIM))+1;
0331   if (z~=zlast)
0332      x=Z(:,:,z);
0333      zlast=z;
0334      disp(['Now analyzing data on slice ' num2str(z) '...']);
0335   end
0336 
0337   b=PinvDesMtx(:,:,z)*x*y;
0338   Beta(:,V(i))=b;
0339   yprime_full=x'*Beta(:,V(i));
0340   T(V(i))=(y-mean(y))'*(y-mean(y));
0341   E(V(i))=(y-yprime_full)'*(y-yprime_full);
0342   H(V(i))=T(V(i))-E(V(i));
0343 
0344   Lambda(V(i))=det(E(V(i)))./det(T(V(i)));
0345   if((Lambda(V(i))<0.0)|(Lambda(V(i))>1.0))
0346      disp(['Lambda out of range on slice ' num2str(z) '.']);
0347      Lambda(V(i))=1.0;
0348   end
0349 
0350   F_omnibus(V(i))=((N-VARIABLES-1)./VARIABLES).*((1-Lambda(V(i)))./Lambda(V(i)));
0351   
0352   yprime_reduced=x'*(Beta(:,V(i)).*S');
0353   EplusH_reduced=(y-yprime_reduced)'*(y-yprime_reduced);
0354   Lambda_reduced(V(i))=det(E(V(i)))./det(EplusH_reduced);
0355   if((Lambda_reduced(V(i))<0.0)|(Lambda_reduced(V(i))>1.0))
0356      disp(['Lambda_reduced out of range on slice ' num2str(z) '.']);
0357      Lambda_reduced(V(i))=1.0;
0358   end
0359   F_reduced(V(i))=((N-VARIABLES-1)./(VARIABLES-sum(S))).*...
0360 ((1-Lambda_reduced(V(i)))./Lambda_reduced(V(i)));
0361 
0362 
0363 [CONTRASTS A]=size(C);
0364 
0365 
0366   for j=1:CONTRASTS,
0367      if(sum(C(j,:))>0)                           % e.g. Sum of coeffs GT 0.0, assume a Reduced Model Vector
0368        yprime_j=x'*(Beta(:,V(i)).*C(j,:)');
0369        EplusH_j=(y-yprime_j)'*(y-yprime_j);
0370        Lambda_j=det(E(V(i)))./det(EplusH_j);
0371      elseif(abs(sum(C(j,:)))<TINY)               % e.g. Abs sum of coeffs <= TINY, Assume a Contrast Model Vector
0372        yprime_j=x'*(Beta(:,V(i)).*C(j,:)');
0373        H_j=yprime_j'*yprime_j;
0374        E_j=y'*y-yprime_j'*yprime_j;
0375        EplusH_j=(E_j+H_j);
0376        Lambda_j=det(E_j)./det(EplusH_j);
0377      else
0378        disp(['Contrast ' num2str(j) ' appears to be invalid.  Please check.']);
0379      end
0380 
0381      if((Lambda_j<0.0)|(Lambda_j>1.0))
0382         disp(['Lambda_j ' num2str(Lambda_j) ' out of range on slice ' num2str(z) '.']);
0383         Lambda_j=1.0;
0384      end
0385      
0386      F_indiv(j,V(i))=((N-VARIABLES-1)./1).*((1-Lambda_j)./Lambda_j);
0387 
0388      if(F_INDIV_PLOTS==1)
0389        if(Lambda_j<0.5)
0390          subplot(1,2,1);
0391          plot(x(j,:),(y-yprime_j),'.');
0392          hold on;
0393          y_reg=x(j,:).*Beta(j,V(i));
0394          plot(x(j,:),y_reg,'r');
0395          axis square;
0396          xlabel('X_j');
0397          ylabel('Y_i');
0398          title(['F_indiv #' num2str(j) ' = ' num2str(F_indiv(j,V(i))) ' (Lambda_j = ' num2str(Lambda_j) ') on slice ' num2str(z) ' at voxel ' num2str(V(i))]);
0399          hold off;
0400          subplot(1,2,2);
0401          plot((y-yprime_j),'g.');
0402          hold on;
0403          plot(y_reg,'r');
0404      xlabel('Time');
0405          ylabel('Y_i');
0406          title(['Beta_j = ' num2str(Beta(j,V(i))) ', R-sqr = ' num2str(1-Lambda_j)]);
0407          hold off;
0408          drawnow;
0409          pause(0.5);
0410          print VoxelPlots.eps -dpsc -append;
0411       end
0412      end
0413 
0414   end
0415 
0416   % Periodically display program elapsed and projected time info to stdout
0417   if(mod(i,XDIM*YDIM)==0)
0418     pctcomplete = floor(100.*i./length(V));
0419     disp([ num2str(pctcomplete) ' percent complete...']);
0420     disp(['Program Elapsed Running Time = ' num2str(floor(toc/60)) ':' num2str(60*mod((toc/60),1)) ]);
0421     runtime=toc./(pctcomplete/100);
0422     timeleft=runtime-toc;
0423     disp(['Projected to run for ' num2str(floor(timeleft/60)) ':' num2str(60*mod((timeleft/60),1)) ' more minutes']);
0424   end
0425 
0426 end
0427 
0428 % **********************************
0429 % Critical F Value for Reduced Model
0430 Fcrit=Fcritical(pcrit,(VARIABLES-sum(S)),(N-VARIABLES-1));
0431 Rsqr_crit=(Fcrit*(VARIABLES-sum(S))./(N-VARIABLES-1))./(1+(Fcrit*(VARIABLES-sum(S))./(N-VARIABLES-1)));
0432 disp([' ']);
0433 disp(['Critical Reduced Model F(' num2str(VARIABLES-sum(S)) ',' num2str(N-VARIABLES-1) ')=' num2str(Fcrit) ' at p=' num2str(pcrit)]);
0434 disp(['Model R-Squared at Critical F = ' num2str(Rsqr_crit)]);
0435 disp([' ']);
0436 
0437 % Threshold F_Reduced Image at Critical F
0438 %F_reduced_mask(V)=((F_reduced(V)==Fcrit)|(F_reduced(V)>Fcrit));
0439 %for i=1:CONTRASTS,
0440 %      F_indiv(i,:)=F_indiv(i,:).*F_reduced_mask;
0441 %end
0442 
0443 % *******************
0444 % Save Results to PWD
0445 disp('Saving Results Images...');
0446 
0447 % ************************
0448 % Beta Parameter Estimates
0449 for i=1:VARIABLES,
0450 fid=fopen(['Beta_' num2str(i) '.img'],'w',ENDIANOUT);
0451 fwrite(fid,Beta(i,:),'float32');
0452 fclose(fid);
0453 str=['!makeheader4 ', ['Beta_' num2str(i) '.hdr '] , num2str(XDIM),...
0454 ' ', num2str(YDIM), ' ', num2str(ZDIM),...
0455 ' ', num2str(XVOXDIM), ' ', num2str(YVOXDIM), ' ', num2str(ZVOXDIM),...
0456 ' 1 32 32676 -32676 1 0 16'];
0457 eval(str);
0458 end
0459 
0460 % **************
0461 % Total Variance
0462 fid=fopen(['TotalVariance.img'],'w',ENDIANOUT);
0463 fwrite(fid,T,'float32');
0464 fclose(fid);
0465 str=['!makeheader4 TotalVariance.hdr ', num2str(XDIM), ' ', num2str(YDIM), ' ', num2str(ZDIM),...
0466 ' ', num2str(XVOXDIM), ' ', num2str(YVOXDIM), ' ', num2str(ZVOXDIM),...
0467 ' 1 32 32676 0 1 0 16'];
0468 eval(str);
0469 
0470 % ****************************
0471 % Reduced Model Total Variance
0472 fid=fopen(['ReducedModelVariance.img'],'w',ENDIANOUT);
0473 fwrite(fid,EplusH_reduced,'float32');
0474 fclose(fid);
0475 str=['!makeheader4 ReducedModelVariance.hdr ', num2str(XDIM), ' ', num2str(YDIM), ' ', num2str(ZDIM),...
0476 ' ', num2str(XVOXDIM), ' ', num2str(YVOXDIM), ' ', num2str(ZVOXDIM),...
0477 ' 1 32 32676 0 1 0 16'];
0478 eval(str);
0479 
0480 % *************
0481 % Wilk's Lambda
0482 fid=fopen(['Lambda.img'],'w',ENDIANOUT);
0483 fwrite(fid,Lambda,'float32');
0484 fclose(fid);
0485 str=['!makeheader4 Lambda.hdr ', num2str(XDIM), ' ', num2str(YDIM), ' ', num2str(ZDIM),...
0486 ' ', num2str(XVOXDIM), ' ', num2str(YVOXDIM), ' ', num2str(ZVOXDIM),...
0487 ' 1 32 1.0 0.0 1 0 16'];
0488 eval(str);
0489 
0490 % *********
0491 % Omnibus F
0492 fid=fopen(['F_' num2str(VARIABLES) '_' num2str(N-VARIABLES-1) '.img'],'w',ENDIANOUT);
0493 fwrite(fid,F_omnibus,'float32');
0494 fclose(fid);
0495 str=['!makeheader4 ', ['F_' num2str(VARIABLES) '_' num2str(N-VARIABLES-1) '.hdr '],...
0496 num2str(XDIM), ' ', num2str(YDIM), ' ', num2str(ZDIM),...
0497 ' ', num2str(XVOXDIM), ' ', num2str(YVOXDIM), ' ', num2str(ZVOXDIM),...
0498 ' 1 32 32676 0 1 0 16'];
0499 eval(str);
0500 
0501 % ******************************************************
0502 % Wilks Reduced Model: Corrected for Confounding Effects
0503 fid=fopen(['Lambda_reduced.img'],'w',ENDIANOUT);
0504 fwrite(fid,Lambda_reduced,'float32');
0505 fclose(fid);
0506 str=['!makeheader4 ', ['Lambda_reduced.hdr '],...
0507 num2str(XDIM), ' ', num2str(YDIM), ' ', num2str(ZDIM),...
0508 ' ', num2str(XVOXDIM), ' ', num2str(YVOXDIM), ' ', num2str(ZVOXDIM),...
0509 ' 1 32 1.0 0.0 1 0 16'];
0510 eval(str);
0511 
0512 % **************************************************
0513 % F Reduced Model: Corrected for Confounding Effects
0514 fid=fopen(['F_reduced_' num2str(VARIABLES-sum(S)) '_' num2str(N-VARIABLES-1) '.img'],'w',ENDIANOUT);
0515 fwrite(fid,F_reduced,'float32');
0516 fclose(fid);
0517 str=['!makeheader4 ', ['F_reduced_' num2str(VARIABLES-sum(S)) '_' num2str(N-VARIABLES-1) '.hdr '],...
0518 num2str(XDIM), ' ', num2str(YDIM), ' ', num2str(ZDIM),...
0519 ' ', num2str(XVOXDIM), ' ', num2str(YVOXDIM), ' ', num2str(ZVOXDIM),...
0520 ' 1 32 32676 0 1 0 16'];
0521 eval(str);
0522 
0523 % ***************************************
0524 % F Reduced Model: Thresholded Mask Image
0525 fid=fopen(['F_reduced_mask.img'],'w',ENDIANOUT);
0526   fwrite(fid,F_reduced_mask,'int8');
0527   fclose(fid);
0528   str=['!makeheader4 F_reduced_mask.hdr ', num2str(XDIM), ' ', num2str(YDIM), ' ', num2str(ZDIM),...
0529   ' ', num2str(XVOXDIM), ' ', num2str(YVOXDIM), ' ', num2str(ZVOXDIM), ' 1 8 1 0 1 0 2'];
0530   eval(str);
0531 
0532 % *******************************************
0533 % Individual Regressor and CONTRASTS F Values
0534 for j=1:CONTRASTS,
0535    fid=fopen(['F_regressor_' num2str(j) '_' num2str(1) '_' num2str(N-VARIABLES-1) '.img'],'w',ENDIANOUT);
0536    fwrite(fid,F_indiv(j,:),'float32');
0537    fclose(fid);
0538    str=['!makeheader4 ', ['F_regressor_' num2str(j) '_' num2str(1) '_' num2str(N-VARIABLES-1) '.hdr '],...
0539 num2str(XDIM), ' ', num2str(YDIM), ' ', num2str(ZDIM),...
0540 ' ', num2str(XVOXDIM), ' ', num2str(YVOXDIM), ' ', num2str(ZVOXDIM),...
0541 ' 1 32 32676 0 1 0 16'];
0542    eval(str);
0543 end
0544 
0545 % ***************************************************
0546 % For Display, Plot Results of Maximum Omnibus F-test
0547 MaxF_reduced=find(F_reduced==max(F_reduced));
0548 z=floor(MaxF_reduced./(XDIM*YDIM))+1;
0549 x(:,:)=Z(:,:,z);
0550 yprime=x'*Beta(:,MaxF_reduced);
0551 
0552 figure;
0553 y=Y(:,MaxF_reduced)-mean(Y(:,MaxF_reduced));
0554 plot(y,'r.');
0555 hold
0556 plot(yprime,'b');
0557 hold off
0558 xlabel('Time');
0559 ylabel('Normalized MR Units');
0560 title(['Overall Fitted Model at Voxel ' num2str(MaxF_reduced) ' on slice ' num2str(z)]);
0561 
0562 % **********************************
0563 % Generate Plot of Corrected Results
0564 Ycorrected=(Y(:,MaxF_reduced)-mean(Y(:,MaxF_reduced)))-x'*(Beta(:,MaxF_reduced).*S');
0565 yprime=x'*Beta(:,MaxF_reduced)-x'*(Beta(:,MaxF_reduced).*S');
0566 
0567 figure;
0568 plot(Ycorrected,'r.');
0569 hold
0570 plot(yprime,'b');
0571 hold off
0572 title(['MR Signal Corrected for Run to Run Effects at Voxel ' num2str(MaxF_reduced) ' on slice ' num2str(z)]);
0573 xlabel('Time');
0574 ylabel('Normalized MR Units');
0575 
0576 
0577 % *******************************************************************************
0578 % Find Local Maxima for Reduced Model Image and Write Results Tables to Text File
0579 % and Time Series to a Matlab file
0580 
0581 if(LOCAL_MAX_PLOTS==1)
0582  !rm LocalMaxPlots.eps;
0583 end
0584 
0585 outfile='TrackingLocalMax.txt';
0586 eval(['!rm ' outfile]);
0587 disp(' ');
0588 disp('Now computing parameter image FWHM and identifiying local maxima...');
0589 
0590 G=[];
0591 R=[];
0592 prob=1.0;
0593 
0594 for i=1:length(S),
0595    if(S(i)==1)
0596       [Sx, Sy, Sz, U, Resels]=ComputeResels2(Beta(i,:), mask, XDIM, YDIM, ZDIM);
0597       G=[G;[Sx,Sy,Sz]];
0598       R=[R;Resels];
0599    end
0600 end
0601 Sxmean=mean(G(:,1));
0602 Symean=mean(G(:,2));
0603 Szmean=mean(G(:,3));
0604 
0605 R=length(V)./((Sxmean*Symean*Szmean).^(2/3));
0606 %B=(4/3)*pi*sqrt(det(S./2));
0607 Reselsmean=mean(R);
0608 LocalMax=FindLocalMax(F_reduced, XDIM, YDIM, ZDIM, Sxmean, Symean, Szmean, Reselsmean, Fcrit);
0609 disp([num2str(length(LocalMax)) ' local maxima identified']);
0610 PeakReducedModelTimeseries=zeros(N,length(LocalMax));
0611 
0612 if(LOCAL_MAX_PLOTS==1)  % Make a new figure if local max plots are desired
0613          figure;
0614          !rm LocalMaxPlots.*
0615 end
0616 
0617 [P Q]=size(LocalMax);
0618 
0619 for i=1:P,
0620 
0621       xlocation=LocalMax(i,2);
0622       ylocation=LocalMax(i,3);
0623       Index=LocalMax(i,5);
0624       %zlocation=floor(Index./(XDIM*YDIM))+1;
0625       zlocation=LocalMax(i,4);
0626       
0627       ANOVATables(outfile,xlocation,ylocation,zlocation,...
0628       N,VARIABLES,Beta(:,Index),T(Index),...
0629       F_omnibus(Index),(VARIABLES-sum(S)),F_reduced(Index),...
0630       Lambda(Index),Lambda_reduced(Index),EplusH_reduced,F_indiv(:,Index));
0631       
0632       x=Z(:,:,zlocation);
0633       yprime=x'*(Beta(:,Index).*S');
0634       y=Y(:,Index)-mean(Y(:,Index));
0635       PeakReducedModelTimeseries(:,i)=(y-yprime);
0636       yint=x'*(Beta(:,Index).*(1-S)');
0637 
0638       if(LOCAL_MAX_PLOTS==1)
0639          
0640          subplot(1,2,1);
0641          plot((y-yprime),'r.');
0642          hold on;
0643          plot(yint,'b');
0644          xlabel('Time');
0645          ylabel('Y');
0646          title(['Local Max F Reduced ' num2str(i) ' = ' num2str(F_reduced(Index)) ' (Lambda Reduced = ' num2str(Lambda_reduced(Index)) ') on voxel ' num2str([xlocation ylocation zlocation])]);
0647          hold off;
0648          %drawnow;
0649          
0650          subplot(1,2,2);
0651          plot((y-yprime),yint,'g.');
0652          %hold on;
0653          %b=linspace(min(y-yprime),max(y-yprime),10);
0654          %plot(b,b,'r');
0655      xlabel('Y');
0656          ylabel('Y Predicted');
0657          title(['R-sqr = ' num2str(1-Lambda_reduced(Index))]);
0658          axis square;
0659          hold off;
0660          drawnow;
0661          %pause(0.5);
0662          print LocalMaxPlots.eps -dpsc -append;
0663       end
0664 
0665 end
0666 
0667 B=Beta(:,LocalMax(:,5));
0668 save PeakReducedModelTimeseries PeakReducedModelTimeseries LocalMax Z B C;
0669 
0670 % ***************************************************************************
0671 % Fill Structure for Christian Buechel's Structural Equation Modeling Routine
0672 %
0673 % Struct array Data
0674 % -----------------
0675 % Data(k).X      - n x m matrix of time-series of m regions
0676 % Data(k).L      - 3 x m matrix for centers of VOI for m regions
0677 % Data(k).Rname     - m x l matrix containing string descriptors (length l)
0678 %           for m regions
0679 % Data(k).useit     - index vector, indicating which regions of X are used
0680 %           (max(useit) <= m)
0681 % Data(k).FlagOb - Vector of size(useit) indicating whether a variable
0682 %           is observed (=1) or latent (=0)
0683 % Data(k).mod     - Flag to indicate whether backwards modulatory connections
0684 %           should be introduced for all (driving) connections
0685 % Data(k).RT     - Repetition time for experiment (to correct for
0686 %           autocorrelation)
0687 % Data(k).U      - 2 x g matrix coding g unilateral connections
0688 %           (U(1,:) = sources and U(2,:) = destination)
0689 % Data(k).B      - 2 x f matrix coding f bidirectional connections
0690 %           between (B(1,:) and B(2,:)
0691 %
0692 %
0693 % Cell array C
0694 % --------------
0695 % C{h}(f)     - struct array cantaining f constraints
0696 % C{h}(f).value     - constrain the following paths to value,
0697 %           if value = Inf impose equality constraint
0698 % C{h}(f).conn     - 2 x n matrix, (1,n) is index into 'Data(k)',
0699 %           (2,n) denotes which paths within "Data"
0700 %           !note, counting starts with U and then B
0701 %
0702 % Struct Misc
0703 % -----------
0704 % Misc.Output      - Degree of output wanted: 0 = none, 1 = some, 2 = all tables,
0705 % Misc.Descr     - String to describe the analysis
0706 % Misc.random     - Use random starting estimates
0707 %
0708 %
0709 % OUTPUT:
0710 %
0711 % Struct array G
0712 % --------------
0713 % G(h).value     - constrain the following paths to value,
0714 %           if value = Inf impose equality constraint
0715 % G(h).conn     - 2 x n matrix, (1,n) is index into Data array (k),
0716 %           (2,n) denotes which paths within "Data"
0717 %           !note, counting starts with U and then B
0718 % G(h).Free      - Matrix of free parameters
0719 % G(h).xfix      - fixed parameters
0720 % G(h).SEM     - see SEM
0721 % G(h).x0        - starting parameters for optimisation
0722 % G(h).chi_sq     - chi squared goodness of fit
0723 % G(h).df     - degrees of freedom
0724 % G(h).p     - p value
0725 % G(h).RMSEA     - RMSEA fit index
0726 % G(h).ECVI     - scale AIC fit index
0727 %
0728 % Struct array SEM
0729 % ----------------
0730 % SEM(k).Ustring    - Strings describing uni-directional connections
0731 % SEM(k).Bstring    - Strings describing bi-directional connections
0732 % SEM(k).Rname        - Strings describing regions (only updated
0733 %                  if modulatory connections are modelled)
0734 % SEM(k).Constring    - Strings describing all
0735 % SEM(k).Fil         - Filter to calculate covariances with latent variables
0736 % SEM(k).Cov        - covariance matrix
0737 % SEM(k).df         - degrees of freedom
0738 % SEM(k).ConX         - matrix of unidirectional paths
0739 % SEM(k).ConZ         - matrix of bidirectional paths
0740 % SEM(k).A        - asymmetric path coefficients
0741 % SEM(k).S        - symmetric path coefficients
0742 % SEM(k).AL        - Modification indices for asymmetric path coefficients (*)
0743 % SEM(k).SL        - Modification indices for symmetric path coefficients
0744 % SEM(k).f        - fit index from optimisation
0745 % SEM(k).Res        - residual covariances
0746 % SEM(k).Est         - estimated covariances
0747 % SEM(k).L         - modification indices for each parameter of
0748 %
0749 if (SEM==1)
0750 
0751 Rname=[];
0752 for i=1:P,
0753   if(i<10)
0754    reg=['00' num2str(i)];
0755   elseif((i>9)&(i<100))
0756    reg=['0' num2str(i)];
0757   elseif((i>99)&(i<1000))
0758    reg=[num2str(i)];
0759   end
0760   Rname=[Rname; reg];
0761 end
0762 
0763      st = struct('X', PeakReducedModelTimeseries(:,1:P)',...
0764          'L', LocalMax(:,2:4),...
0765          'Rname',Rname,...
0766          'useit', [1:P],...
0767          'FlagOb', ones(size([1:P])),...
0768          'mod', 1,...
0769          'RT', TR,...
0770                  'U', [1 2;...
0771               1 2],...             
0772          'B', [1 2;...
0773                    1 2]);
0774      C1 =   struct('value',{[]},...
0775           'conn',{[]});
0776      C2 =   struct(    'value',{Inf,-.5},...
0777           'conn',{[1 2;...
0778                  1 1],...
0779                [3;...
0780                 1]});
0781      clear C;
0782      Data(1)=st;
0783      C{1}=C1;
0784      C{2}=C2;
0785      Misc =   struct('Output',2,...
0786           'Descr' ,'Reduced Model Regression Local Max Data',...
0787           'random',0);
0788 
0789      save SEM Data C Misc;
0790 
0791 end
0792 
0793 % *******
0794 % The End
0795 disp('All Done.');
0796 disp(['Total Program Running Time = ' num2str(floor(toc/60)) ':' num2str(60*mod(toc/60,1)) ]);
0797 
0798 
0799 
0800 
0801 % %%%%%%%%% %
0802 % JVH, 2001 %
0803 % %%%%%%%%% %
0804 
0805 
0806 
0807 
0808 
0809 
0810
TrackingAnalysis

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SOURCE CODE
