https://github.com/jdiedrichsen/pcm_toolbox
Tip revision: 4e290a8b2c0d0820f868b7bcb60a3da7bb30e6ee authored by Jörn Diedrichsen on 26 April 2023, 01:59:24 UTC
Update pcm_estimateRegression.m
Update pcm_estimateRegression.m
Tip revision: 4e290a8
pcm_fitModelIndivid.m
function [T,theta_hat,G_pred,INFO]=pcm_fitModelIndivid(Y,M,partitionVec,conditionVec,varargin);
% function [T,theta_hat,G_pred]=pcm_fitModelIndivid(Y,M,partitionVec,conditionVec,varargin);
% Fits pattern component model(s) specified by M to data from a number of
% subjects.
% The model parameters are all individually fit.
%==========================================================================
% INPUT:
% Y: {#Subjects}.[#Conditions x #Voxels]
% Observed/estimated beta regressors from each subject.
% Preferably multivariate noise-normalized beta regressors.
%
% M: {#Models} Cell array with structure that defines model(s). Each
% may contain the fields
% .type: Type of the model to be fitted
% 'fixed': Fixed structure without parameter (except scale for each subject)
% 'component': G is a sum of linear components
% 'feature': G=A*A', with A a linear sum of weighted feature components
% 'nonlinear': Nonlinear model with own function to return derivatives
% 'freechol': Free model in Cholesky form
% .numGparams: Scalar that defines the number of parameters
% included in model.
% .theta0: Vector of starting theta parameters to calculate predicted
% model G.
% for more fields see the manual for model specification.
%
% partitionVec: {#Subjects} Cell array with partition assignment vector
% for each subject. Rows of partitionVec{subj} define
% partition assignment of rows of Y{subj}.
% Commonly these are the scanning run #s for beta
% regressors.
% If a single vector is provided, it is assumed to be the
% same for all subjects
%
% conditionVec: {#Subjects} Cell array with condition assignment vector
% for each subject. Rows of conditionVec{subj} define
% condition assignment of rows of Y{subj}.
% If a single vector is provided, it is assumed to me the
% same for all subjects.
% If the (elements of) conditionVec are matrices, it is
% assumed to be the design matrix Z, allowing the
% specification individualized models.
%
%--------------------------------------------------------------------------
% OPTION:
% 'runEffect': How to deal with effects that may be specific to different
% imaging runs:
% 'random': Models variance of the run effect for each subject
% as a seperate random effects parameter.
% 'fixed': Consider run effect a fixed effect, will be removed
% implicitly using ReML.
% 'none': No modeling of the run Effect (not recommended for fMRI data)
%
% 'isCheckDeriv: Check the derivative accuracy of theta params. Done using
% 'checkderiv'. This function compares input to finite
% differences approximations. See function documentation.
%
% 'fitAlgorithm': Either 'NR' or 'minimize' - provides over-write for
% model specific algorithms
% 'maxIteration': Number of max minimization iterations. Default is 1000.
%
% 'likeThres': Threshold for minimal change in the log-likelihood to
% stop iterations
% 'verbose': Optional flag to show display message in the command
% line. Default is 1. Setting to 2 gives more detailed
% feedback on pcm_NR
%
% 'S': Optional specific covariance structure of the noise
% {#Subjects} = cell array of N_s x N_s matrices
% or
% S(#Subjects).S and .invS: Structure of the N_s x N_s
% normal and inverse covariances matrices
%
% 'theta0': Cell array of starting values (same format as theta{m})
% 'fitScale': Fit a additional scale parameter for each subject?
% (0/1). Default is set to 0.
%--------------------------------------------------------------------------
% OUTPUT:
% T: Structure with following subfields:
% SN: Subject number
% likelihood: log-likelihood
% scale: Scale parameter (if fitscale = 1)-exp(theta_s)
% noise: Noise parameter- exp(theta_eps)
% run: Run parameter (if run = 'random')
% iterations: Number of interations for model fit
% time: Elapsed time in sec
%
% theta{m} Cell array of estimated model parameters, each a
% #params x #numSubj matrix
% G_pred{m} Cell array of estimated G-matrices under the model
runEffect = 'fixed';
isCheckDeriv = 0;
maxIteration = []; % MaxIteration - otherwise algorithm-specific defaults
likeThres = []; % if log-likelihood decreases less that this value, it stops (for pcm_NR)
verbose = 1; % 1: Indicating the subject 2: detailed feedback
S = [];
fitAlgorithm = [];
theta0 = [];
fitScale = 0;
pcm_vararginoptions(varargin,{'runEffect','isCheckDeriv','maxIteration',...
'verbose','S','fitAlgorithm','theta0','fitScale','likeThres'});
numSubj = numel(Y);
% Set options for fitting algorithm
fitOPT.verbose = verbose;
if (~isempty(maxIteration))
fitOPT.maxIteration = maxIteration;
end
if (~isempty(likeThres))
fitOPT.likeThres = likeThres;
end
% Determine number of models
if (~iscell(M))
M={M};
end
numModels = numel(M);
% Preallocate output structure
T.SN = [1:numSubj]';
% Determine optimal algorithm for each of the models
if (~isempty(fitAlgorithm))
for m=1:numModels
M{m}.fitAlgorithm = fitAlgorithm;
end
end
M = pcm_optimalAlgorithm(M);
% Set up all parameters for the upcoming fit
[Z,B,X,YY,S,N,P,G_hat,noise0,run0]=...
pcm_setUpFit(Y,partitionVec,conditionVec,'runEffect',runEffect,'S',S);
% Loop over subject and provide inidivdual fits
for s = 1:numSubj
% Now loop over models
for m = 1:length(M)
if (verbose)
if isfield(M{m},'name');
fprintf('Fitting Subj: %d model:%s\n',s,M{m}.name);
else
fprintf('Fitting Subj: %d model:%d\n',s,m);
end
end
tic;
% Get starting guess for theta0 is not provided
if (~isfield(M{m},'theta0'))
M{m}.theta0=pcm_getStartingval(M{m},G_hat(:,:,s));
end
% If naive noise ceiling model, use G_hat
if strcmp(M{m}.type,'freedirect')
G0 = pcm_makePD(mean(G_hat,3));
M{m}.numGparams=0;
M{m}.Gc = G0;
else
G0 = pcm_calculateG(M{m},M{m}.theta0);
end
g0 = G0(:);
% If fitting scale parameterm get initial estimate using OLS
if (fitScale)
g_hat = G_hat(:,:,s);
g_hat = g_hat(:);
scaling = (g0'*g_hat)/(g0'*g0);
if ((scaling<10e-5)||~isfinite(scaling)); scaling = 10e-5; end % Enforce positive scaling
scale0(s,m) = log(scaling);
end
if (numel(theta0)<m || size(theta0{m},2)<s)
th0m = [M{m}.theta0;noise0(s)];
if(fitScale)
th0m = [th0m;scale0(s,m)];
end
if (strcmp(runEffect,'random'))
th0m = [th0m;run0(s)];
end
theta0{m}(:,s)=th0m;
end;
% Set options for likelihood function
OPT.fitScale = fitScale;
OPT.runEffect=B{s};
if (~isempty(S))
OPT.S = S(s);
end;
% Now do the fitting, using the preferred optimization routine
switch (M{m}.fitAlgorithm)
case 'minimize' % Use minimize to find maximum liklhood estimate runEffect',B{s});
fcn = @(x) pcm_likelihoodIndivid(x,YY{s},M{m},Z{s},X{s},P(s),OPT);
if (isempty(maxIteration))
maxIteration=1000;
end
[theta_hat{m}(:,s),fX,T.iterations(s,m)] = ...
minimize(theta0{m}(:,s), fcn, maxIteration);
T.likelihood(s,m) = -fX(end); %invert the sign
case 'NR'
fcn = @(x) pcm_likelihoodIndivid(x,YY{s},M{m},Z{s},X{s},P(s),OPT);
[theta_hat{m}(:,s),T.likelihood(s,m),T.iterations(s,m),T.reg(s,m),INFO.regH{s,m},INFO.thetaH{s,m}]= ...
pcm_NR(theta0{m}(:,s),fcn,fitOPT);
end
G_pred{m}(:,:,s) = pcm_calculateG(M{m},theta_hat{m}(1:M{m}.numGparams,s));
T.noise(s,m) = exp(theta_hat{m}(M{m}.numGparams+1,s));
if (fitScale)
T.scale(s,m) = exp(theta_hat{m}(M{m}.numGparams+2,s));
end
if strcmp(runEffect,'random')
T.run(s,m) = exp(theta_hat{m}(M{m}.numGparams+2+fitScale,s));
end
T.time(s,m) = toc;
% This is an optional check if the dervivate calculation is correct
if (isCheckDeriv)
d = pcm_checkderiv(fcn,theta_hat{m}(:,s)-0.001,0.00001);
fprintf('discrepency :%d\n',d);
keyboard;
end
end % for each model
end % for each subject
INFO.theta0=theta0;