diff --git a/src/f_SCR.m b/src/f_SCR.m index 192f6802..5fac18c9 100644 --- a/src/f_SCR.m +++ b/src/f_SCR.m @@ -51,7 +51,7 @@ sts = -1; sigma = 0.3; % std for event-related and spontaneous sudomotor input function sigma_SCL = 1;% std for SCL changes -sigma_offset = settings.dcm{1}.sigma_offset; % offset for aSCR sigma to constrain the amplitude/sd trade-off +sigma_offset = settings.dcm{1}.aSCR_sigma_offset; % offset for aSCR sigma to constrain the amplitude/sd trade-off Theta_n = 7; % number of parameters for the output function (the rest is for the input function) try diff --git a/src/pspm_cfg/pspm_cfg_dcm.m b/src/pspm_cfg/pspm_cfg_dcm.m index eb507131..9c644491 100644 --- a/src/pspm_cfg/pspm_cfg_dcm.m +++ b/src/pspm_cfg/pspm_cfg_dcm.m @@ -143,12 +143,13 @@ lasttrialfiltering.strtype = 'r'; lasttrialfiltering.num = [1 1]; lasttrialfiltering.help = pspm_cfg_help_format('pspm_dcm', 'model.lasttrialcutoff'); + % constrained model constrained_model = cfg_entry; -constrained_model.name ='Constrained model'; -constrained_model.tag ='constr_model'; -constrained_model.val = {0}; -constrained_model.strtype = 'i'; +constrained_model.name = 'Maximum Dispersion (Standard Deviation) for Flexible Responses [s]'; +constrained_model.tag ='constr_model'; % is mapped to model.constrained in pspm_cfg_run_dcm +constrained_model.val = {0.3}; +constrained_model.strtype = 'r'; constrained_model.num = [1 1]; constrained_model.help = pspm_cfg_help_format('pspm_dcm', 'model.constrained'); @@ -272,7 +273,7 @@ % minimum dispersion (standard deviation) for flexible responses ascr_sigma_offset = cfg_entry; ascr_sigma_offset.name = 'Minimum Dispersion (Standard Deviation) for Flexible Responses [s]'; -ascr_sigma_offset.tag = 'ascr_sigma_offset'; +ascr_sigma_offset.tag = 'ascr_sigma_offset'; % is mapped to options.aSCR_sigma_offset in pspm_cfg_run_dcm ascr_sigma_offset.strtype = 'r'; ascr_sigma_offset.num = [1 1]; ascr_sigma_offset.val = {0.1}; diff --git a/src/pspm_cfg/pspm_cfg_run_dcm.m b/src/pspm_cfg/pspm_cfg_run_dcm.m index 33e012a1..25d51b51 100644 --- a/src/pspm_cfg/pspm_cfg_run_dcm.m +++ b/src/pspm_cfg/pspm_cfg_run_dcm.m @@ -64,8 +64,8 @@ if isstruct(filter) model.filter = filter; end -% normalization, subsession threshold -model = pspm_update_struct(model, job.data_options, {'norm', 'substhresh'}); +% normalization, subsession threshold, lasttrialcutoff +model = pspm_update_struct(model, job.data_options, {'norm', 'substhresh', 'lasttrialcutoff'}); % constrained model model.constrained = job.data_options.constr_model; % channel number @@ -77,6 +77,10 @@ else options.rf = job.resp_options.rf.datafile; end + + + + % options options = pspm_update_struct(options, job.resp_options, {'crfupdate',... 'indrf',... @@ -86,10 +90,19 @@ 'sfpost',... 'sffreq',... 'sclpre',... - 'sclpost',... - 'ascr_sigma_offset',... - 'dispwin',... - 'dispsmallwin'}); + 'sclpost'}); + +if isfield(job, 'disp_options') + options = pspm_update_struct(options, job.disp_options, {'dispwin','dispsmallwin'}); +end + +% accept both cfg tag and pspm_dcm option field name +if isfield(job.inv_options, 'aSCR_sigma_offset') % in case user orients by pspm_dcm + options.aSCR_sigma_offset = job.inv_options.aSCR_sigma_offset; +elseif isfield(job.inv_options, 'ascr_sigma_offset') + options.aSCR_sigma_offset = job.inv_options.ascr_sigma_offset; +end + % condition and event names if isfield(options, 'trlnames') options.trlnames = [options.trlnames{:}]; % collapse over sessions diff --git a/src/pspm_check_model.m b/src/pspm_check_model.m index bab423f4..3c56a3be 100644 --- a/src/pspm_check_model.m +++ b/src/pspm_check_model.m @@ -112,8 +112,12 @@ % ├.substhresh: [optional, DCM (modeltype) only] Minimum duration (in seconds) of NaN % │ periods to cause splitting up into subsessions which get evaluated % │ independently (excluding NaN values). Default as 2. -% ├.constrained: [optional, DCM (modeltype) only] Constrained model for flexible -% │ responses which have fixed dispersion (0.3 s SD) but flexible latency. +% ├.constrained: [optional, DCM (modeltype) only] Constrain dispersion of flexible responses. +% │ The value stated here refers to the SD in seconds. +% │ If set to the default of 0.3 s, dispersion is assumed to +% │ be fixed. Otherwise, for each estimated response, the upper +% │ limit of the dispersion is the minimum of this value and +% │ 1/2 the duration of the flexible response window. % └────method : [optional, SF (modeltype) only] [string/cell_array] % [string] either 'auc', 'scl', 'dcm' (default), or 'mp'. % [cell_array] a cell array of methods mentioned above. @@ -312,9 +316,9 @@ end if ~isfield(model, 'constrained') - model.constrained = 0; - elseif ~any(ismember(model.constrained, [0, 1])) - warning('ID:invalid_input', 'Constrained model must be specified as 0 or 1.'); return; + model.constrained = 0.3; + elseif ~isnumeric(model.constrained) || model.constrained < settings.aSCR_sigma_offset + warning('ID:invalid_input', 'Flexible response dispersion must be numeric and larger than the minimum value specified in the settings, which is %0.1f s.', settings.aSCR_sigma_offset); return; end if ~isfield(model, 'substhresh') diff --git a/src/pspm_dcm.m b/src/pspm_dcm.m index cdb4b140..c871e455 100644 --- a/src/pspm_dcm.m +++ b/src/pspm_dcm.m @@ -11,9 +11,14 @@ % (evoked after a specified event) responses can be modelled. % For fixed responses, delay and dispersion are assumed to be constant % (either pre-determined or estimated from the data), while for flexible -% responses, both are estimated for each individual trial. +% responses, latency (and if selected, dispersion) are estimated for each +% individual trial. % Flexible responses can for example be anticipatory, decision-related, % or evoked with unknown onset. +% If flexible response dispersion is estimated from the data, we +% recommend multiplying it with the estimated amplitude to obtain a +% valid estimate of the psychological variable (see de Vries et al. +% 2026). % PsPM implements an iterative trial-by-trial algorithm. Different from % GLM, response parameters are always estimated per trial, and the % algorithm is not informed about the condition. @@ -74,8 +79,12 @@ % │ i.e. Data are normalised during inversion but results % │ transformed back into raw data units. % │ Default: 0. -% └─.constrained: [optional] Constrained model for flexible responses -% which have fixed dispersion (0.3 s SD) but flexible latency. +% └─.constrained: [optional] Constrain dispersion of flexible responses. +% The value stated here refers to the SD in seconds. +% If set to the default of 0.3 s, dispersion is assumed to +% be fixed. Otherwise, for each estimated response, the upper +% limit of the dispersion is the minimum of this value and +% 1/2 the duration of the flexible response window. % ┌────options % ├─.crfupdate: [0/1] Re-estimate RF parameters from canonical SCRF, % │ or use pre-estimated RF parameters. This can be used when @@ -113,9 +122,6 @@ % │ Default: 2 s. % ├───.sclpost: SCR-change-free interval after last event of a trial. % │ Default: 5 s. -% ├.aSCR_sigma_offset: -% │ Minimum dispersion (standard deviation) for flexible -% │ responses, in seconds. Default: 0.1 s. % ├─.dispwin: [0/1] Display progress plot. Default: display. % ├─.dispsmallwin: [0/1] % │ Display intermediate progress plots. @@ -178,10 +184,14 @@ % modelling of anticipatory skin conductance changes. Biological % Psychology, 85(1), 163-70 % [2] Model validation and improvement: -% Staib, M., Castegnetti, G., & Bach, D. R. (2015). Optimising a +% Staib M, Castegnetti G, Bach DR (2015). Optimising a % model-based approach to inferring fear learning from skin % conductance responses. Journal of Neuroscience Methods, 255, % 131-138. +% [3] Optimisation of the forward model: +% de Vries O, Mancinelli F, Sporrer JK, Bach DR (2026). Optimal +% quantification of fear conditioning from skin conductance data. +% preprint available at https://osf.io/preprints/psyarxiv/abe4t_v1 % ● History % Introduced in PsPM 3.0 % Written in 2010-2021 by Dominik R Bach (Wellcome Centre for Human Neuroimaging, UCL) @@ -212,7 +222,7 @@ % all the below should be re-factored into pspm_options ------------------- % numeric fields num_fields = {'depth', 'sfpre', 'sfpost', 'sffreq', 'sclpre', ... - 'sclpost', 'aSCR_sigma_offset'}; + 'sclpost'}; % logical fields bool_fields = {'crfupdate', 'indrf', 'getrf', 'dispwin', ... 'dispsmallwin', 'nosave'}; @@ -718,6 +728,26 @@ return end +%% temp rf solution +% In RF-only mode, pspm_dcm_inv returns RF posterior information but no +% trial-wise session field dcm.sn. Therefore skip stats/name assembly. +if options.getrf + dcm.dcmname = model.modelfile; + dcm.modelfile = model.modelfile; + dcm.input = model; + dcm.options = options; + dcm.warnings = warnings; + dcm.modeltype = 'dcm'; + dcm.modality = settings.modalities.dcm; + + if ~options.nosave + save(model.modelfile, 'dcm'); + end + + sts = 1; + return +end + %% 7 Assemble stats & names dcm.stats = []; cTrl = 0; diff --git a/src/pspm_dcm_inv.m b/src/pspm_dcm_inv.m index 924b2f0c..45cd5d6d 100644 --- a/src/pspm_dcm_inv.m +++ b/src/pspm_dcm_inv.m @@ -28,8 +28,12 @@ % ├─.missing_data: [optional] missing epoch data, originally loaded as model.missing % │ from pspm_dcm, but calculated into .missing_data (created % │ in pspm_dcm and then transferred to pspm_dcm_inv. -% └──.constrained: [optional] constrained model for flexible responses which have -% have fixed dispersion (0.3 s SD) but flexible latency +% └─.constrained: [optional] Constrain dispersion of flexible responses. +% The value stated here refers to the SD in seconds. +% If set to the default of 0.3 s, dispersion is assumed to +% be fixed. Otherwise, for each estimated response, the upper +% limit of the dispersion is the minimum of this value and +% 1/2 the duration of the flexible response window. % ┌───────options % ├─────────.eSCR: [optional] contains the data to estimate RF from % ├─────────.aSCR: [optional] contains the data to adjust the RF to @@ -51,10 +55,6 @@ % │ scl-change-free window before first event. % ├──────.sclpost: [optional, numeric, default: 5, unit: second] % │ scl-change-free window after last event. -% ├─.aSCR_sigma_offset: -% │ [optional, numeric, default: 0.1, unit: second] -% │ minimum dispersion (standard deviation) for flexible -% │ responses. % ├──────.dispwin: [optional, bool, default as 1] % │ display progress window. % └─.dispsmallwin: [optional, bool, default as 0] @@ -109,8 +109,7 @@ try model.trlstop; catch, warning('Trial ends not defined.'); return; end try model.iti; catch, warning('ITIs not defined.'); return; end try model.norm; catch, model.norm = 0; end -try model.constrained; catch, model.constrained = 0; end -try model.constrained_upper; catch, model.constrained_upper = 0; end +try model.constrained; catch, model.constrained = 0.3; end try model.aSCR; catch, model.aSCR = 0; end @@ -131,8 +130,8 @@ end try invopt.DisplayWin = options.dispwin; catch, invopt.DisplayWin = 1; end try invopt.GnFigs = options.dispsmallwin; catch, invopt.GnFigs = 0; end -sigma_offset_temp = settings.dcm{1}.sigma_offset; -try settings.dcm{1}.sigma_offset = options.aSCR_sigma_offset; catch; end + +sigma_offset = settings.dcm{1}.aSCR_sigma_offset; % offset for aSCR sigma to constrain the amplitude/sd trade-off % set general priors and initial conditions % ------------------------------------------------------------------------- @@ -140,7 +139,6 @@ sftheta = pspm_sf_theta; sf_unit = 1./exp(sftheta(5)); sftheta = sftheta(1:3); -fixedSD = 0.3; % CRF priors generated on 27.04.2010 -- % numeric values given in log(parameter space) such that these @@ -184,11 +182,7 @@ % aSCR priors -- prior.aTheta.m = zeros(1, aSCRno); -if model.constrained - prior.aTheta.s = 100 * ones(1, aSCRno); -else - prior.aTheta.s = zeros(1, aSCRno); -end +prior.aTheta.s = zeros(1, aSCRno); prior.aTheta.a = log(0.25) * ones(1, aSCRno); % shorten variable names -- @@ -385,14 +379,8 @@ for k = 1:aSCRno u(5 + k, :) = model.flexevents(k, 1); u(5 + aSCRno + k, :) = model.flexevents(k, 2); % aSCR mean upper bound - if model.constrained - u(5 + 2 * aSCRno + k, :) = fixedSD - settings.dcm{1}.sigma_offset; % aSCR SD upper bound - elseif model.constrained_upper - u(5 + 2 * aSCRno + k, :) = model.constrained_upper - settings.dcm{1}.sigma_offset; % aSCR SD upper bound - else - u(5 + 2 * aSCRno + k, :) = diff(model.flexevents(k, :))/2 - settings.dcm{1}.sigma_offset; % aSCR SD upper bound - end - end + u(5 + 2 * aSCRno + k, :) = min(diff(model.flexevents(k, :))/2, model.constrained) - sigma_offset; % aSCR SD upper bound + end for k = 1:eSCRno u(5 + 3 * aSCRno + k, :) = model.fixevents(k); % eSCR onset end @@ -406,8 +394,9 @@ % output function parameters are now fixed for n = 1:theta_n, priors.SigmaTheta(n, n) = 0; end - % if model constrained, flexible response dispersion is fixed - if model.constrained + % if flexible response dispersion is fixed, set prior precision to + % infinity + if model.constrained == 0.3 aSCRindx = theta_n + 3 * ((1:aSCRno) - 1) + 2; for n = 1:theta_n, priors.SigmaTheta(n, n) = 0; end end @@ -604,13 +593,8 @@ u(5 + u(2, 1) + (1:u(2, 1)), :) = repmat(foo(:), 1, size(u, 2)); aSCR_ln(1:aSCRno, trl) = foo(:, 1); % save first trial for transformation of parameter values into seconds % - get aSCR SD upper bound (zero for dummy events, fixed SD for constrained models) - if model.constrained - u(5 + 2 * u(2, 1) + (1:u(2, 1)), :) = repmat(fixedSD, numel(foo), size(u, 2)) - settings.dcm{1}.sigma_offset; - elseif model.constrained_upper > 0 - u(5 + 2 * u(2, 1) + (1:u(2, 1)), :) = repmat(model.constrained_upper, numel(foo), size(u, 2)) - settings.dcm{1}.sigma_offset; - else - u(5 + 2 * u(2, 1) + (1:u(2, 1)), :) = repmat(foo(:)/2, 1, size(u, 2)) - settings.dcm{1}.sigma_offset; - end + upper_bound = min([foo(:)/2, repmat(model.constrained, [numel(foo), 1])], [], 2); + u(5 + 2 * u(2, 1) + (1:u(2, 1)), :) = repmat(upper_bound, [1, size(u, 2)]) - sigma_offset; % tidy up clear aSCR_on foo aSCR_dummy else @@ -838,12 +822,8 @@ for k = 1:aSCRno sig.G0 = aSCR_ln(k, trl); aTheta(trl).m(k) = sigm(aTheta(trl).m(k), sig); - if model.constrained - sig.G0 = fixedSD - settings.dcm{1}.sigma_offset; - else - sig.G0 = aSCR_ln(k, trl)/2 - settings.dcm{1}.sigma_offset; - end - aTheta(trl).s(k) = sigm(aTheta(trl).s(k), sig) + settings.dcm{1}.sigma_offset; + sig.G0 = min(aSCR_ln(k, trl)/2, model.constrained) - sigma_offset; + aTheta(trl).s(k) = sigm(aTheta(trl).s(k), sig) + sigma_offset; end aTheta(trl).a = newzfactor .* exp(aTheta(trl).a) ./ eSCR_unit; eTheta(trl).a = newzfactor .* exp(eTheta(trl).a) ./ eSCR_unit; @@ -916,7 +896,6 @@ %% (7) clear up % ======================================================================== -settings.dcm{1}.sigma_offset = sigma_offset_temp; dcm.invmodel = model; sts = 1; return diff --git a/src/pspm_get_rf.m b/src/pspm_get_rf.m index 52fb1d9f..1403e1bb 100644 --- a/src/pspm_get_rf.m +++ b/src/pspm_get_rf.m @@ -14,11 +14,12 @@ % SPM style onsets file with one event type, or (3) an epochs file % (see pspm_dcm or pspm_get_epochs). % * outfile : (optional) a file to write the response function to -% * channel : (optional) data channel (default: look for first SCR channel) +% * channel : (optional) data channel (default: look for first SCR channel) % * options : [struct] to be passed on to pspm_dcm % ● History % Introduced in PsPM 3.0 % Written in 2008-2015 by Dominik R Bach (Wellcome Trust Centre for Neuroimaging) +% Updated in 2026 by Bernhard A. von Raußendorf %% initialise global settings @@ -26,7 +27,8 @@ pspm_init; end sts = -1; -rf = []; +rf = []; % in the m-function can be removed +theta = []; %% check input if nargin < 1 @@ -40,28 +42,73 @@ [pth infn ext] = fileparts(fn); outfile = fullfile(pth, ['RF_', infn, ext]); end -if nargin < 4 +if nargin < 4 || isempty(channel) channel = 'scr'; end +if nargin < 5 || isempty(options) + options = struct(); +end + %% call DCM options = pspm_options(options, 'get_rf'); % options.getrf = 1; %try options.nosave, catch, options.nosave = 1; end -options.channel = channel; -[foo dcm] = pspm_dcm(fn, '', events, options); -if numel(dcm{1}.prior.posterior) == 2 +if options.invalid + warning('options invalid') % change! + return; +end + +% prepare timing for pspm_dcm -> maybe change! +if isnumeric(events) + if isvector(events) + events = {events(:)}; + else + events = {events}; + end +end + +model.datafile = {fn}; +model.timing = {events}; +model.channel = channel; +model.modelfile = [tempname, '.mat']; + +[dsts, dcm] = pspm_dcm(model, options); +if dsts < 1 || isempty(dcm) + warning('RF estimation failed in pspm_dcm.'); % change! + return; +end + +% options.channel = channel; +% [foo dcm] = pspm_dcm(fn, '', events, options); + + +if iscell(dcm) + dcm_rf = dcm{1}; +else + dcm_rf = dcm; +end + +if numel(dcm_rf.prior.posterior) == 2 % based on eSCR - theta = dcm{1}.prior.posterior(2).muTheta(1:7)'; + theta = dcm_rf.prior.posterior(2).muTheta(1:7)'; else % based on aSCR (i. e. updated RF) - theta = dcm{1}.prior.posterior(3).muTheta(1:7)'; + theta = dcm_rf.prior.posterior(3).muTheta(1:7)'; end +% if numel(dcm{1}.prior.posterior) == 2 +% % based on eSCR +% theta = dcm{1}.prior.posterior(2).muTheta(1:7)'; +% else +% % based on aSCR (i. e. updated RF) +% theta = dcm{1}.prior.posterior(3).muTheta(1:7)'; +% end + %% write response function to file -if ~isempty(outfile) +if ~isempty(outfile) % should never be empty! [pth fn ext] = fileparts(outfile); - c = clock; + c = clock; % timestamp job{1} = sprintf('function [rf, theta] = %s(td)', fn); job{2} = '%-----------------------------------------------------------------------'; job{3} = ['% Response function created by pspm_get_rf, ', date, sprintf(' %02.0f:%02.0f', c(4:5))]; @@ -82,9 +129,17 @@ job{18} = sprintf('rf = Xt(1, :);'); job{19} = sprintf('rf = rf/max(rf);'); job{20} = sprintf('rf = rf(:);'); - job = strvcat(job'); + % job = strvcat(job'); + % outfile = fullfile(pth, [fn, '.m']); + % dlmwrite(outfile, job, 'delimiter', ''); outfile = fullfile(pth, [fn, '.m']); - dlmwrite(outfile, job, 'delimiter', ''); + + fid = fopen(outfile, 'w'); + fprintf(fid, '%s\n', job{:}); + fclose(fid); end + sts = 1; return + +end diff --git a/src/pspm_init.m b/src/pspm_init.m index ab3c4e2a..26aff0d8 100644 --- a/src/pspm_init.m +++ b/src/pspm_init.m @@ -740,7 +740,7 @@ % analysis. Further modalities and models can be implemented. % % DCM for SCR filter settings -defaults.dcm{1} = struct('filter', struct('lpfreq', 5, 'lporder', 1, 'hpfreq', 0.0159, 'hporder', 1, 'down', 10, 'direction', 'bi'), 'sigma_offset', 0.3); +defaults.dcm{1} = struct('filter', struct('lpfreq', 5, 'lporder', 1, 'hpfreq', 0.0159, 'hporder', 1, 'down', 10, 'direction', 'uni'), 'aSCR_sigma_offset', 0.1); % DCM for SF filter settings defaults.dcm{2} = struct('filter', struct('lpfreq', 5, 'lporder', 1, 'hpfreq', 0.0159, 'hporder', 1, 'down', 10, 'direction', 'uni')); %% 8 TAM settings diff --git a/src/pspm_options.m b/src/pspm_options.m index 2c925507..c15c9344 100644 --- a/src/pspm_options.m +++ b/src/pspm_options.m @@ -119,8 +119,6 @@ options = autofill(options, 'overwrite', 2, [0,1] ); case 'dcm' % 2.16 pspm_dcm -- - options = autofill(options, 'aSCR_sigma_offset', 0.1, '>=', 0 ); - % minimum dispersion (standard deviation) for flexible responses (second) options = autofill(options, 'crfupdate', 0, '*Num' ); % update CRF priors to observed SCRF, or use pre-estimated priors options = autofill(options, 'crfupdate', 0, 1 ); @@ -156,8 +154,6 @@ % condition descriptions) case 'dcm_inv' % 2.17 pspm_dcm_inv -- - options = autofill(options, 'aSCR_sigma_offset', 0.1, '*Num' ); - % minimum dispersion (standard deviation) for flexible responses (second) options = autofill(options, 'crfupdate', 0, '*Num' ); % update CRF priors to observed SCRF, or use pre-estimated priors, default to use pre-estimated priors options = autofill(options, 'crfupdate', 0, 1 ); diff --git a/test/pspm_dcm_test.m b/test/pspm_dcm_test.m index b115ace7..52b29d6c 100644 --- a/test/pspm_dcm_test.m +++ b/test/pspm_dcm_test.m @@ -75,7 +75,7 @@ function invalid_input(this) % options % numeric fields num_fields = {'depth', 'sfpre', 'sfpost', 'sffreq', 'sclpre', ... - 'sclpost', 'aSCR_sigma_offset'}; + 'sclpost'}; for f = 1:numel(num_fields) fl = num_fields{f}; values = {'a', {}};