pr_estimate

 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

 Originally written by Jean-Baptiste Poline
 $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 % Originally written by Jean-Baptiste Poline
0013 % $Id$
0014 
0015 %-Say hello
0016 %-----------------------------------------------------------------------
0017 Finter   = spm('FigName','Stats: estimation...'); spm('Pointer','Watch');
0018     
0019 %-------------------------------------------------------------------------
0020 %- set the methods
0021 
0022 COV_estim = 'assumed';           % covariance is assumed to be imposed by filter K
0023 GLM_resol = 'OLS';               % ordinary least square
0024  
0025 
0026 %-------------------------------------------------------------------------
0027 %- get the design structure, the filter and the data
0028 xX = SPM.xX;  
0029 if ~isfield(xX,'X'), 
0030   error('The design does not contain a design matrix');
0031 end
0032 
0033 % allow matrix or object to be passed as input data
0034 marsY = marsy(marsY);
0035 Y     = summary_data(marsY);
0036 n_roi = size(Y,2);  %- Y is a time by n_roi matrix
0037 
0038 % Remove columns with no variance
0039 in_cols = any(diff(Y));
0040 if ~any(in_cols), error('No variance to estimate model'); end
0041 Y = Y(:, in_cols);
0042 
0043 % We are going to ignore AR(1) options
0044 if mars_struct('isthere', xX, 'xVi', 'Form')
0045   if ~strcmp(xX.xVi.Form, 'none')
0046     warning(['Sorry, we are going to ignore autocorrelation option: ' ...
0047          xX.xVi.Form]);
0048   end
0049 end
0050 
0051 %----------------------------------------------------------------------------------
0052 %- Estimation of the covariance structure of the Ys
0053 fprintf('\nEstimating covariance...');
0054 switch COV_estim
0055 
0056     case {'AR(p)'}
0057         %- compute the temporal cov of Y (V) with AR(p)
0058         if ~isfield(xX,'xVi')
0059            xX.xVi = struct(    'Vi', speye(size(xX.X,1)),...
0060                     'Form',    'AR(p)'); 
0061         end
0062         % xX.xVi = estimate_cov(Y,xX);
0063 
0064 
0065     case {'assumed'}
0066         if ~isfield(xX,'xVi')
0067            xX.xVi = struct(    'Vi', speye(size(xX.X,1)),...
0068                     'Form',    'none'); 
0069         end
0070         %- else, the covariance structure is supposed to be
0071         %- stored in xX.xVi
0072 
0073     otherwise
0074         warning('COV_estim does not exist');
0075 end
0076 fprintf('Done\n');
0077 
0078 switch GLM_resol
0079 
0080     case {'OLS'}
0081         fprintf('Using OLS\n');
0082         %- no filter already defined
0083         if ~isfield(xX,'K')
0084            xX.K  = speye(size(xX.X,1));
0085         end
0086         % else assume that the filter is xX.K
0087     
0088         KVi    = pr_spm_filter('apply', xX.K, xX.xVi.Vi); 
0089         V      = pr_spm_filter('apply', xX.K, KVi'); 
0090         Y      = pr_spm_filter('apply', xX.K, Y);
0091         fprintf('Setting filter...');
0092         KXs    = spm_sp('Set', pr_spm_filter('apply', xX.K, xX.X));
0093         fprintf('Done.\n');
0094         clear KVi;
0095 
0096     case {'MaxLik'} 
0097         'MaxLik' 
0098         %- compute the inverse filter -  put it in K ?
0099         %- filter data and design
0100         %- V = speye(size(xX.X,1));
0101 
0102     otherwise
0103         warning('GLM_resol does not exist');
0104 end
0105 
0106 %- compute GLM
0107 fprintf('Computing estimates...');
0108 if ~spm_sp('isspc',KXs), Xs = spm_sp('set',KXs); else Xs = KXs;  end
0109 
0110 [trRV trRVRV] = spm_SpUtil('trRV',Xs,V); 
0111 beta          = spm_sp('x-', Xs, Y);                 %-Parameter estimates
0112 res           = spm_sp('r', Xs, Y);                  %-Residuals
0113 ResMS         = sum(res.^2)./trRV;                   %-Res variance estimation
0114 xX.erdf       = trRV^2/trRVRV;
0115 
0116 fprintf('Done.\n');
0117 
0118 % fill up design related stuff
0119 xX.V     = V;                                     %-V matrix
0120 xX.xKXs  = KXs;                            %-Filtered design matrix
0121 xX.pKX   = spm_sp('x-',xX.xKXs);
0122 xX.pKXV  = xX.pKX*xX.V;                %-for contrast variance weight
0123 xX.Bcov  = xX.pKXV*xX.pKX';            %-Variance of est. param.
0124 [xX.trRV,xX.trRVRV] ...                %-Variance expectations
0125          = spm_SpUtil('trRV',xX.xKXs,xX.V);
0126 xX.nKX   = spm_DesMtx('sca',xX.xKXs.X,xX.Xnames);% scaled design matrix for display
0127 
0128 % Put back into design
0129 nBeta                    = size(xX.X, 2);
0130 SPM.betas                = ones(nBeta, n_roi) + NaN;
0131 SPM.betas(:, in_cols)    = beta;    
0132 SPM.ResidualMS           = ones(1, n_roi) + NaN;
0133 SPM.ResidualMS(in_cols)  = ResMS;    
0134 
0135 SPM.xX    = xX;
0136 SPM.marsY = marsY;
0137 
0138 %-Default F-contrasts (in contrast structure)
0139 %=======================================================================
0140 F_iX0 = [];
0141 if isfield(SPM, 'F_iX0')
0142   F_iX0 = SPM.F_iX0;
0143 end
0144 if isempty(F_iX0)
0145   F_iX0 = struct(    'iX0',        [],...
0146             'name',        'all effects');
0147 elseif ~isstruct(F_iX0)
0148   F_iX0 = struct(    'iX0',        F_iX0,...
0149             'name',        'effects of interest');
0150 end
0151 
0152 %-Create Contrast structure array
0153 %-----------------------------------------------------------------------
0154 xCon  = spm_FcUtil('Set',F_iX0(1).name,'F','iX0',F_iX0(1).iX0,xX.xKXs);
0155 for i = 2:length(F_iX0)
0156     xcon = spm_FcUtil('Set',F_iX0(i).name,'F','iX0',F_iX0(i).iX0,xX.xKXs);
0157     xCon = [xCon xcon];
0158 end
0159 SPM.xCon = xCon;
0160 
0161 %=======================================================================
0162 %- E N D: Cleanup GUI
0163 %=======================================================================
0164 spm_progress_bar('Clear')
0165 spm('FigName','Stats: done',Finter); spm('Pointer','Arrow')
0166 fprintf('%-40s: %30s\n','Completed',spm('time'))                     %-#
0167 fprintf('...use the results section for assessment\n\n')             %-#
0168 
0169 return