private function to compute mutlivariate statistics across ROIs FORMAT [MVres] = pr_stat_compute_mv(SPM,Ic) Input SPM - SPM design structure Ic - indices into contrast structure (xCon in SPM) Output MVres - mulitvariate result structure $Id: pr_stat_compute_mv.m 92 2004-01-07 08:37:16Z matthewbrett $
0001 function [MVres] = pr_stat_compute_mv(SPM,Ic) 0002 % private function to compute mutlivariate statistics across ROIs 0003 % FORMAT [MVres] = pr_stat_compute_mv(SPM,Ic) 0004 % 0005 % Input 0006 % SPM - SPM design structure 0007 % Ic - indices into contrast structure (xCon in SPM) 0008 % 0009 % Output 0010 % MVres - mulitvariate result structure 0011 % 0012 % $Id: pr_stat_compute_mv.m 92 2004-01-07 08:37:16Z matthewbrett $ 0013 0014 %-Get contrast definitions (if available) 0015 %----------------------------------------------------------------------- 0016 try 0017 xCon = SPM.xCon; 0018 catch 0019 xCon = []; 0020 end 0021 0022 %-set all contrasts by default 0023 %----------------------------------------------------------------------- 0024 if nargin < 2 0025 Ic = 1:length(xCon); 0026 end 0027 if any(Ic > length(xCon)) 0028 error('Indices too large for contrast structure'); 0029 end 0030 0031 0032 % Get relevant fields from design 0033 xCon = xCon(Ic); 0034 Xs = SPM.xX.xKXs; 0035 V = SPM.xX.V; 0036 betas = SPM.betas; 0037 ResidualMS = SPM.ResidualMS; 0038 Y = summary_data(SPM.marsY); 0039 0040 % setup calculation 0041 [nBetas nROI] = size(betas); 0042 nCon = length(xCon); 0043 [trRV trRVRV] = spm_SpUtil('trRV',Xs,V); 0044 erdf = trRV^2/trRVRV; 0045 RMS = sqrt(ResidualMS); 0046 0047 %-------------------------------------------------------------------- 0048 %- Multivariate analysis 0049 %-------------------------------------------------------------------- 0050 0051 MVres = struct('y_pre',[], 'y_obs', [], 'Pf', [], 'u', [], 'ds', [] ); 0052 0053 if nCon == 1, return, end 0054 0055 YpY = Y'*Y; 0056 0057 for ii = 1:nCon 0058 0059 xC = xCon(ii); 0060 0061 %-------------------------------------------------------------------- 0062 [NF, nu, h, d, M12, XG, sXG] = sf_model_mlm(Xs, V, nROI, xC, erdf); 0063 0064 %-------------------------------------------------------------------- 0065 %- Compute svd 0066 %-------------------------------------------------------------------- 0067 %- fprintf('%-40s\n','Computing Principal Components') 0068 0069 Z = ((NF*betas)./(ones(size(NF,1),1)*RMS)); 0070 S = Z*Z'; 0071 S = S/sum(nROI); 0072 [u s u] = svd(S,0); 0073 ds = diag(s); 0074 clear s; 0075 0076 0077 %-------------------------------------------------------------------- 0078 %- STATISTICS if any ... 0079 %-------------------------------------------------------------------- 0080 %- Fq : F values for the last q eigein values. 0081 %- P : P values.for the last q eigein values. 0082 0083 Fq = zeros(1,h); 0084 for q = 0:h-1; 0085 nu1 = d*(h-q); 0086 nu2 = d*nu - (d-1)*(4*(h-q)+2*nu)/(h-q+2); 0087 Fq(q+1) = ((nu-2)/nu) * nu2/(nu2-2)*sum(ds(q+1:h))/(h-q); 0088 end 0089 Pf = 1 - spm_Fcdf(Fq,nu1,nu2); 0090 0091 0092 %- fprintf('%-40s\n','Computing predicted and observed temporal reponse') 0093 %keyboard 0094 0095 y_pre = (pinv(XG)'* M12 * u)*diag(sqrt(ds)); % predicted temporal reponse 0096 0097 gV = (diag(1./sqrt(ds))*Z)'*u; 0098 y_obs = (Y./(ones(size(Y,1),1)*RMS)/nROI)*gV; 0099 0100 %- save results for this constrast 0101 MVres(ii).y_pre = y_pre; 0102 MVres(ii).y_obs = y_obs; 0103 MVres(ii).Pf = Pf; 0104 MVres(ii).u = u; 0105 MVres(ii).ds = ds; 0106 MVres(ii).df = [nu1 nu2]; 0107 0108 end 0109 0110 0111 0112 0113 0114 0115 %=================================================================== 0116 function [NF,nu,h,d,M12,XG,sXG] = sf_model_mlm(Xs, V, nROI, xC, erdf); 0117 % Set sub-space of interest and the related matrix of normalisation. 0118 % FORMAT [NF,nu,h,d,M12,XG] = mm_model(); 0119 %- nu, h, d : degrees of freedom 0120 %- NF : matrix of normalisation 0121 %=================================================================== 0122 0123 0124 %-------------------------------------------------------------------- 0125 %- SET, COMPUTE,NORMALIZE SPACES OF INTEREST 0126 %-------------------------------------------------------------------- 0127 %- set X10 and XG 0128 %- XG= X -PG(X), PG projection operator on XG (cf. eq 1, 2) 0129 %-------------------------------------------------------------------- 0130 sX1o = spm_sp('set',spm_FcUtil('X1o',xC,Xs)); 0131 sXG = spm_sp('set',spm_FcUtil('X0',xC,Xs)); 0132 X1o = spm_sp('oP',sX1o,Xs.X); 0133 XG = spm_sp('r',sXG,Xs.X); 0134 0135 %- Compute Normalized effexts : M1/2=X'G*V*XG (cf eq 3) 0136 %-------------------------------------------------------------------- 0137 % warning off; 0138 up = spm_sp('ox',sX1o); ; %- PG=up*up' 0139 qi = up'*Xs.X; 0140 sigma = up'*V*up; 0141 M12 = (chol(sigma)*qi)'; 0142 M_12 = pinv(M12); 0143 0144 %- Compute NF : normalise factor (cf eq 4) 0145 %-------------------------------------------------------------------- 0146 NF = M_12*spm_sp('X',Xs)'*spm_sp('r',sXG,spm_sp('X',Xs)); 0147 0148 %- degrees of freedom 0149 %- nROI : number of ROI (corresponds to the number of Resels) 0150 %-------------------------------------------------------------------- 0151 d = nROI*(4*log(2)/pi)^(3/2); 0152 h = sX1o.rk; %-rank of the sub-space of interest. 0153 nu = erdf;