Home > marsbar > @mardo_5 > private > pr_estimate.m

pr_estimate

PURPOSE ^

Estimation of a General Linear Model

SYNOPSIS ^

function SPM = pr_estimate(SPM, marsY)

DESCRIPTION ^

 Estimation of a General Linear Model
 FORMAT SPM = pr_estimate(SPM, marsY)
 Inputs 
 SPM      - SPM design structure
 marsY    - marsY data object, or 2D data (Y) matrix

 Outputs
 SPM      - modified estimated design structure, with data contained as
            field marsY

 Based on spm_spm from spm2:
 @(#)spm_spm.m    2.66 Andrew Holmes, Jean-Baptiste Poline, Karl Friston 03/03/27

 There are some changes in this version
 1) The specified Vi field, can contain either
 a cell array, in which case it is standard SPM covariance components,
 or a struct array, in which case it can specify other methods of
 estimating the covariance, in particular, real AR(n) estimation
 
 2) The design will specify if the covariance should be calculated from
 the summarized time course(s), or from the component voxels, then
 averaged.  Voxel time courses are used by default, and if
 SPM.xVi.cov_calc is set to 'vox', but summary time
 courses can be used by setting SPM.xVi.cov_calc to 'summary'. 
 
 3) Normally, if the W matrix is present, the V matrix should also be
 present.  Because it is boring to calculate the V matrix and then WVW,
 if we know WVW is I, V can be present, but empty, in which case WVW is
 assumed to be I.  It just saves time.
 
 $Id: pr_estimate.m 543 2004-12-14 08:58:05Z matthewbrett $

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function SPM = pr_estimate(SPM, marsY)
0002 % Estimation of a General Linear Model
0003 % FORMAT SPM = pr_estimate(SPM, marsY)
0004 % Inputs
0005 % SPM      - SPM design structure
0006 % marsY    - marsY data object, or 2D data (Y) matrix
0007 %
0008 % Outputs
0009 % SPM      - modified estimated design structure, with data contained as
0010 %            field marsY
0011 %
0012 % Based on spm_spm from spm2:
0013 % @(#)spm_spm.m    2.66 Andrew Holmes, Jean-Baptiste Poline, Karl Friston 03/03/27
0014 %
0015 % There are some changes in this version
0016 % 1) The specified Vi field, can contain either
0017 % a cell array, in which case it is standard SPM covariance components,
0018 % or a struct array, in which case it can specify other methods of
0019 % estimating the covariance, in particular, real AR(n) estimation
0020 %
0021 % 2) The design will specify if the covariance should be calculated from
0022 % the summarized time course(s), or from the component voxels, then
0023 % averaged.  Voxel time courses are used by default, and if
0024 % SPM.xVi.cov_calc is set to 'vox', but summary time
0025 % courses can be used by setting SPM.xVi.cov_calc to 'summary'.
0026 %
0027 % 3) Normally, if the W matrix is present, the V matrix should also be
0028 % present.  Because it is boring to calculate the V matrix and then WVW,
0029 % if we know WVW is I, V can be present, but empty, in which case WVW is
0030 % assumed to be I.  It just saves time.
0031 %
0032 % $Id: pr_estimate.m 543 2004-12-14 08:58:05Z matthewbrett $
0033 
0034 %-Say hello
0035 %-----------------------------------------------------------------------
0036 Finter   = spm('FigName','Stats: estimation...'); spm('Pointer','Watch')
0037 
0038 %=======================================================================
0039 % - A N A L Y S I S   P R E L I M I N A R I E S
0040 %=======================================================================
0041 
0042 %-Initialise
0043 %=======================================================================
0044 fprintf('%-40s: %30s','Initialising parameters','...computing')    %-#
0045 xX            = SPM.xX;
0046 [nScan nBeta] = size(xX.X);
0047 
0048 %-Check confounds (xX.K) and non-sphericity (xVi)
0049 %-----------------------------------------------------------------------
0050 if ~isfield(xX,'K')
0051   xX.K  = 1;
0052 end
0053 try
0054   %-If covariance components are specified use them
0055   %---------------------------------------------------------------
0056   xVi   = SPM.xVi;
0057 catch
0058   
0059   %-otherwise assume i.i.d.
0060   %---------------------------------------------------------------
0061   xVi   = struct(    'form',  'i.i.d.',...
0062             'V',     speye(nScan,nScan));
0063 
0064 end
0065 
0066 % Work out what we are going to do
0067 have_W     = isfield(xX, 'W');
0068 have_V     = isfield(xVi, 'V');
0069 
0070 % Work out type of covariance modelling. We get Vi, cov_type (as a string):
0071 % one of 'SPM' or 'fmristat' and cov_vox, which is a flag set to 1 if all
0072 % the voxel time courses should be used to calculate the resdiduals and
0073 % covariance.
0074 if ~have_V
0075   if ~isfield(xVi, 'Vi')
0076     error('No covariance specified, and no priors to calculate it');
0077   end
0078   Vi = xVi.Vi;
0079   if iscell(Vi)
0080     cov_type = 'SPM';
0081   elseif ~isstruct(Vi)
0082     error('Vi field should be cell or struct type')
0083   elseif ~isfield(Vi, 'type')
0084     error('Vi should have field specifying type');
0085   else
0086     cov_type = Vi.type;
0087   end
0088   
0089   % Covariance calculated on summary or voxel time courses
0090   cov_vox = 1;
0091   if isfield(xVi, 'cov_calc')
0092     cov_vox = strcmpi(xVi.cov_calc, 'vox');
0093   end
0094 else cov_vox = 0; end
0095 
0096 %-Get non-sphericity V
0097 %=======================================================================
0098 if have_V
0099   %-If xVi.V is specified proceed directly to parameter estimation
0100   %---------------------------------------------------------------
0101   V     = xVi.V;
0102   str   = 'parameter estimation';
0103 else
0104   % otherwise invoke ReML selecting voxels under i.i.d assumptions
0105   %---------------------------------------------------------------
0106   V     = speye(nScan,nScan);
0107   str   = '[hyper]parameter estimation';
0108 end
0109 
0110 %-Get whitening/Weighting matrix: If xX.W exists we will save WLS
0111 % estimates. Get WVW also, which can be assumed if W is a whitening
0112 % matrix
0113 %-----------------------------------------------------------------------
0114 if have_W
0115   %-If W is specified, use it
0116   %-------------------------------------------------------
0117   W     = xX.W;
0118   if isempty(V)  % V is only inv(W*W')
0119     WVW = eye(nScan);
0120   else
0121     WVW = W*V*W';
0122   end
0123 else
0124   if have_V
0125     % otherwise make W a whitening filter W*W' = inv(V)
0126     %-------------------------------------------------------
0127     [u s] = pr_spm_svd(xVi.V);
0128     s     = spdiags(1./sqrt(diag(s)),0,nScan,nScan);
0129     W     = u*s*u';
0130     W     = W.*(abs(W) > 1e-6);
0131     xX.W  = W;
0132     WVW   = eye(nScan);
0133   else
0134     % unless xVi.V has not been estimated - requiring 2 passes
0135     %-------------------------------------------------------
0136     W     = eye(nScan);
0137     str   = 'hyperparameter estimation (1st pass)';
0138   end
0139 end
0140 
0141 %-Design space and projector matrix [pseudoinverse] for WLS
0142 %=======================================================================
0143 xX.xKXs = spm_sp('Set',pr_spm_filter(xX.K,W*xX.X));        % KWX
0144 xX.pKX  = spm_sp('x-',xX.xKXs);                % projector
0145 
0146 %-If xVi.V is not defined compute Hsqr
0147 %-----------------------------------------------------------------------
0148 if ~isfield(xVi,'V')
0149   Fcname = 'effects of interest';
0150   iX0    = [SPM.xX.iB SPM.xX.iG];
0151   xCon   = spm_FcUtil('Set',Fcname,'F','iX0',iX0,xX.xKXs);
0152   X1o    = spm_FcUtil('X1o', xCon(1),xX.xKXs);
0153   Hsqr   = spm_FcUtil('Hsqr',xCon(1),xX.xKXs);
0154   trRV   = spm_SpUtil('trRV',xX.xKXs);
0155   trMV   = spm_SpUtil('trMV',X1o);
0156 end
0157 
0158 fprintf('%s%30s\n',repmat(sprintf('\b'),1,30),'...done')               %-#
0159 
0160 %=======================================================================
0161 % - F I T   M O D E L   &   W R I T E   P A R A M E T E R    I M A G E S
0162 %=======================================================================
0163 
0164 % Select whether to work with all voxel data in ROIs, or summary data
0165 % Using all data only makes sense for intial estimation of whitening
0166 if ~have_W & cov_vox
0167   str = 'voxelwise';
0168   Y = region_data(marsY);
0169   Y = [Y{:}];
0170 else
0171   str = 'pooled';
0172   Y = summary_data(marsY);
0173 end
0174 
0175 % Eliminate columns with zero variance
0176 in_cols = any(diff(Y));
0177 if ~any(in_cols), error('No variance to estimate model'); end
0178 Y = Y(:, in_cols);
0179 
0180 fprintf('%-40s: %30s\n','Covariance estimate',['...' str])               %-#
0181 fprintf('%-40s: %30s','Model','...start')    %-#
0182 
0183 n_roi = n_regions(marsY);
0184 
0185 %-Intialise variables used in the loop
0186 %=======================================================================
0187 [n S] = size(Y);                                    % no of time courses
0188 Cy    = 0;                        % <Y*Y'> spatially whitened
0189 CY    = 0;                        % <Y*Y'> for ReML
0190 EY    = 0;                        % <Y>    for ReML
0191 %-Whiten/Weight data and remove filter confounds
0192 %-------------------------------------------------------
0193 fprintf('%s%30s',repmat(sprintf('\b'),1,30),'filtering')    %-#
0194 
0195 KWY   = pr_spm_filter(xX.K,W*Y);
0196 
0197 %-General linear model: Weighted least squares estimation
0198 %------------------------------------------------------
0199 fprintf('%s%30s',repmat(sprintf('\b'),1,30),'estimation') %-#
0200 
0201 beta  = xX.pKX*KWY;            %-Parameter estimates
0202 res   = spm_sp('r',xX.xKXs,KWY);    %-Residuals
0203 ResSS = sum(res.^2);            %-Residual SSQ
0204 clear KWY                %-Clear to save memory
0205 
0206 
0207 %-If ReML hyperparameters are needed for xVi.V
0208 %-------------------------------------------------------
0209 if ~have_V
0210   if n_roi > 1
0211     wstr = {'Pooling covariance estimate across ROIs',...
0212         'This is unlikely to be valid; A better approach',...
0213         'is to run estimation separately for each ROI'};
0214     fprintf('\n');
0215     warning(sprintf('%s\n', wstr{:}));
0216   end
0217   % Cy is whitened covariance matrix; only needed for SPM REML method
0218   if strcmp(cov_type, 'SPM')
0219     q  = diag(sqrt(trRV./ResSS'),0); % spatial whitening
0220     Y  = Y * q;
0221     Cy = Y*Y';
0222   end
0223 end % have_V
0224         
0225 %-if we are saving the WLS parameters
0226 %-------------------------------------------------------
0227 if have_W
0228 
0229   %-sample covariance and mean of Y (all voxels)
0230   %-----------------------------------------------
0231   CY         = Y*Y';
0232   EY         = sum(Y,2);
0233     
0234 end % have_W
0235 clear Y                %-Clear to save memory
0236 
0237 fprintf('\n')                                                        %-#
0238 spm_progress_bar('Clear')
0239 
0240 %=======================================================================
0241 % - P O S T   E S T I M A T I O N   C L E A N U P
0242 %=======================================================================
0243 if S == 0, warning('No time courses - empty analysis!'), end
0244 
0245 %-average sample covariance and mean of Y (over voxels)
0246 %-----------------------------------------------------------------------
0247 CY          = CY/S;
0248 EY          = EY/S;
0249 CY          = CY - EY*EY';
0250 
0251 %-If not defined, compute non-sphericity V using ReML Hyperparameters
0252 %=======================================================================
0253 if ~have_V
0254 
0255   %-Estimate of residual correlations through hyperparameters
0256   %---------------------------------------------------------------
0257   str    = 'Temporal non-sphericity (over voxels)';
0258   fprintf('%-40s: %30s\n',str,'...estimation') %-#
0259   Cy            = Cy/S;
0260   
0261   % Estimate for separable designs and covariance components
0262   %---------------------------------------------------------------
0263   if isstruct(xX.K)
0264     
0265     switch cov_type
0266      case 'SPM'
0267       % Store hyperparameters
0268       m     = length(Vi);
0269       h     = zeros(m,1);
0270      case 'fmristat'
0271       % Store AR coefficients
0272       h     = zeros(length(xX.K), Vi.order);
0273      otherwise 
0274       error(['Did not recognize covariance type: ' cov_type]);
0275     end
0276     
0277     V     = sparse(nScan,nScan); 
0278     for i = 1:length(xX.K)
0279       
0280       % extract blocks from bases
0281       %-----------------------------------------------
0282       q     = xX.K(i).row;
0283       
0284       % design space for estimation (with confounds in filter)
0285       %-----------------------------------------------
0286       Xp         = xX.X(q,:);
0287       try
0288     Xp = [Xp xX.K(i).X0];
0289       end
0290       
0291       switch cov_type
0292        case 'SPM'
0293     % REML: extract blocks from bases
0294     %-----------------------------------------------
0295     p     = [];
0296     Qp    = {};
0297     for j = 1:m
0298       if nnz(xVi.Vi{j}(q,q))
0299         Qp{end + 1} = xVi.Vi{j}(q,q);
0300         p           = [p j];
0301       end
0302     end
0303 
0304     % ReML itself
0305     %-----------------------------------------------
0306     fprintf('%-30s- %i\n','  ReML Block',i);
0307     [Vp,hp]  = pr_spm_reml(Cy(q,q),Xp,Qp);
0308     V(q,q)   = V(q,q) + Vp;
0309     h(p)     = hp;
0310        
0311        case 'fmristat'
0312     % AR estimation
0313     [h(i,:) W(q,q)]  = pr_fmristat_ar(res(q,:),Xp,Vi.order);
0314 
0315       end
0316     end
0317   else
0318     [V,h] = pr_spm_reml(Cy,xX.X,xVi.Vi);
0319   end
0320   
0321   switch cov_type
0322    case 'SPM'
0323     % normalize non-sphericity and save hyperparameters
0324     %---------------------------------------------------------------
0325     V           = V*nScan/trace(V);
0326    case 'fmristat'
0327     % Set covariance matrix to empty, we have already calculated W
0328     %---------------------------------------------------------------
0329     V           = [];
0330     SPM.xX.W    = W;
0331   end
0332   
0333   xVi.h       = h;
0334   xVi.V       = V;            % Save non-sphericity xVi.V
0335   xVi.Cy      = Cy;            %-spatially whitened <Y*Y'>
0336   SPM.xVi     = xVi;            % non-sphericity structure
0337   
0338   % If xX.W is not specified use W*W' = inv(V) to give ML estimators
0339   %---------------------------------------------------------------
0340   if ~have_W
0341     % clear everything except SPM, marsY;
0342     vnames = who;
0343     vnames = vnames(~ismember(vnames, {'SPM','marsY'}));
0344     clear(vnames{:});
0345     SPM = pr_estimate(SPM,marsY);
0346     return
0347   end
0348 end
0349 
0350 
0351 %-Use non-sphericity xVi.V to compute [effective] degrees of freedom
0352 %=======================================================================
0353 xX.V            = pr_spm_filter(xX.K,pr_spm_filter(xX.K,WVW)');    % KWVW'K'
0354 [trRV trRVRV]   = spm_SpUtil('trRV',xX.xKXs,xX.V);        % trRV (for X)
0355 xX.trRV         = trRV;                        % <R'*y'*y*R>
0356 xX.trRVRV       = trRVRV;                    %-Satterthwaite
0357 xX.erdf         = trRV^2/trRVRV;                % approximation
0358 xX.Bcov         = xX.pKX*xX.V*xX.pKX';                % Cov(beta)
0359 
0360 
0361 %-scale ResSS by 1/trRV
0362 %-----------------------------------------------------------------------
0363 ResMS           = ResSS/xX.trRV;
0364 
0365 %-Create 1st contrast for 'effects of interest' (all if not specified)
0366 %=======================================================================
0367 Fcname          = 'effects of interest';
0368 try
0369   iX0     = [xX.iB xX.iG];
0370 catch
0371   iX0     = [];
0372 end
0373 xCon            = spm_FcUtil('Set',Fcname,'F','iX0',iX0,xX.xKXs);
0374 
0375 %-Compute scaled design matrix for display purposes
0376 %-----------------------------------------------------------------------
0377 xX.nKX        = spm_DesMtx('sca',xX.xKXs.X,xX.name);
0378 
0379 %-place fields in SPM
0380 %-----------------------------------------------------------------------
0381 SPM.betas                = ones(nBeta, n_roi) + NaN;
0382 SPM.betas(:, in_cols)    = beta;    
0383 SPM.ResidualMS           = ones(1, n_roi) + NaN;
0384 SPM.ResidualMS(in_cols)  = ResMS;    
0385 
0386 SPM.xVi        = xVi;                % non-sphericity structure
0387 SPM.xVi.CY     = CY;                %-<(Y - <Y>)*(Y - <Y>)'>
0388 
0389 SPM.xX         = xX;                %-design structure
0390 
0391 SPM.xCon       = xCon;                %-contrast structure
0392 
0393 %=======================================================================
0394 %- E N D: Cleanup GUI
0395 %=======================================================================
0396 fprintf('%s%30s\n',repmat(sprintf('\b'),1,30),'...done')               %-#
0397 spm('FigName','Stats: done',Finter); spm('Pointer','Arrow')
0398 fprintf('%-40s: %30s\n','Completed',spm('time'))                     %-#
0399 fprintf('...use the results section for assessment\n\n')             %-#

Generated on Wed 11-May-2022 16:26:09 by m2html © 2003-2019