Non-linear model defaults

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

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};

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

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')

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;

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,16 +130,15 @@
 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
 % -------------------------------------------------------------------------
 % SF priors --
 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