compute multivariate statistics across ROIs FORMAT [MVres]= pr_stat_compute_mv(xCon, Xs, V, betas, ResidualMS, Y) xCon - contrast structure Xs - design matrix V - covariance matrix betas - parameter estimates ResidialMS - mean sum of squares of residuals Y - data natrix Output MVres - result structure $Id$
0001 function [MVres]= pr_stat_compute_mv(xCon, Xs, V, betas, ResidualMS, Y) 0002 % compute multivariate statistics across ROIs 0003 % FORMAT [MVres]= pr_stat_compute_mv(xCon, Xs, V, betas, ResidualMS, Y) 0004 % 0005 % xCon - contrast structure 0006 % Xs - design matrix 0007 % V - covariance matrix 0008 % betas - parameter estimates 0009 % ResidialMS - mean sum of squares of residuals 0010 % Y - data natrix 0011 % 0012 % Output 0013 % MVres - result structure 0014 % 0015 % $Id$ 0016 0017 [nBetas nROI] = size(betas); 0018 nCon = length(xCon); 0019 [trRV trRVRV] = spm_SpUtil('trRV',Xs,V); 0020 erdf = trRV^2/trRVRV; 0021 RMS = sqrt(ResidualMS); 0022 %-------------------------------------------------------------------- 0023 %- Multivariate analysis 0024 %-------------------------------------------------------------------- 0025 0026 MVres = struct('y_pre',[], 'y_obs', [], 'Pf', [], 'u', [], 'ds', [] ); 0027 0028 if nCon == 1, return, end 0029 0030 YpY = Y'*Y; 0031 0032 for ii = 1:nCon 0033 0034 xC = xCon(ii); 0035 0036 %-------------------------------------------------------------------- 0037 [NF, nu, h, d, M12, XG, sXG] = sf_model_mlm(Xs, V, nROI, xC, erdf); 0038 0039 %-------------------------------------------------------------------- 0040 %- Compute svd 0041 %-------------------------------------------------------------------- 0042 %- fprintf('%-40s\n','Computing Principal Components') 0043 0044 Z = ((NF*betas)./(ones(size(NF,1),1)*RMS)); 0045 S = Z*Z'; 0046 S = S/sum(nROI); 0047 [u s u] = svd(S,0); 0048 ds = diag(s); 0049 clear s; 0050 0051 0052 %-------------------------------------------------------------------- 0053 %- STATISTICS if any ... 0054 %-------------------------------------------------------------------- 0055 %- Fq : F values for the last q eigein values. 0056 %- P : P values.for the last q eigein values. 0057 0058 Fq = zeros(1,h); 0059 for q = 0:h-1; 0060 nu1 = d*(h-q); 0061 nu2 = d*nu - (d-1)*(4*(h-q)+2*nu)/(h-q+2); 0062 Fq(q+1) = ((nu-2)/nu) * nu2/(nu2-2)*sum(ds(q+1:h))/(h-q); 0063 end 0064 Pf = 1 - spm_Fcdf(Fq,nu1,nu2); 0065 0066 0067 %- fprintf('%-40s\n','Computing predicted and observed temporal reponse') 0068 %keyboard 0069 0070 y_pre = (pinv(XG)'* M12 * u)*diag(sqrt(ds)); % predicted temporal reponse 0071 0072 gV = (diag(1./sqrt(ds))*Z)'*u; 0073 y_obs = (Y./(ones(size(Y,1),1)*RMS)/nROI)*gV; 0074 0075 %- save results for this constrast 0076 MVres(ii).y_pre = y_pre; 0077 MVres(ii).y_obs = y_obs; 0078 MVres(ii).Pf = Pf; 0079 MVres(ii).u = u; 0080 MVres(ii).ds = ds; 0081 MVres(ii).df = [nu1 nu2]; 0082 0083 end 0084 0085 0086 0087 0088 0089 0090 %=================================================================== 0091 function [NF,nu,h,d,M12,XG,sXG] = sf_model_mlm(Xs, V, nROI, xC, erdf); 0092 % Set sub-space of interest and the related matrix of normalisation. 0093 % FORMAT [NF,nu,h,d,M12,XG] = mm_model(); 0094 %- nu, h, d : degrees of freedom 0095 %- NF : matrix of normalisation 0096 %=================================================================== 0097 0098 0099 %-------------------------------------------------------------------- 0100 %- SET, COMPUTE,NORMALIZE SPACES OF INTEREST 0101 %-------------------------------------------------------------------- 0102 %- set X10 and XG 0103 %- XG= X -PG(X), PG projection operator on XG (cf. eq 1, 2) 0104 %-------------------------------------------------------------------- 0105 sX1o = spm_sp('set',spm_FcUtil('X1o',xC,Xs)); 0106 sXG = spm_sp('set',spm_FcUtil('X0',xC,Xs)); 0107 X1o = spm_sp('oP',sX1o,Xs.X); 0108 XG = spm_sp('r',sXG,Xs.X); 0109 0110 %- Compute Normalized effexts : M1/2=X'G*V*XG (cf eq 3) 0111 %-------------------------------------------------------------------- 0112 % warning off; 0113 up = spm_sp('ox',sX1o); ; %- PG=up*up' 0114 qi = up'*Xs.X; 0115 sigma = up'*V*up; 0116 M12 = (chol(sigma)*qi)'; 0117 M_12 = pinv(M12); 0118 0119 %- Compute NF : normalise factor (cf eq 4) 0120 %-------------------------------------------------------------------- 0121 NF = M_12*spm_sp('X',Xs)'*spm_sp('r',sXG,spm_sp('X',Xs)); 0122 0123 %- degrees of freedom 0124 %- nROI : number of ROI (corresponds to the number of Resels) 0125 %-------------------------------------------------------------------- 0126 d = nROI*(4*log(2)/pi)^(3/2); 0127 h = sX1o.rk; %-rank of the sub-space of interest. 0128 nu = erdf;