Home > marsbar > @mardo_99 > private > pr_stat_compute_mv.m

pr_stat_compute_mv

PURPOSE ^

compute multivariate statistics across ROIs

SYNOPSIS ^

function [MVres]= pr_stat_compute_mv(xCon, Xs, V, betas, ResidualMS, Y)

DESCRIPTION ^

 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$

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

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;

Generated on Thu 17-Mar-2016 01:14:46 by m2html © 2005