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$
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$ 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') %-#