Package {rsDCM}

Type: Package
Title: Robust and Sparse Dynamic Causal Modelling for Functional MRI
Version: 0.1.0
Date: 2026-09-24
Maintainer: Godfred Arhin <arhin0122@gmail.com>
Description: Provides a robust and sparse method for group-level Dynamic Causal Modelling (DCM) of functional magnetic resonance imaging (fMRI) data: Student-t weighting of subjects for robustness, combined with a nonlocal product-moment (pMOM) spike-and-slab prior for sparse selection of group-level effects (<doi:10.48550/arXiv.2609.06379>). The package also provides an R implementation of single-subject DCM for fMRI using variational Laplace inversion (Friston et al., 2003 <doi:10.1016/S1053-8119(03)00202-7>), including the bilinear neural state equation and the Buxton-Friston hemodynamic response model, ported from the 'SPM25' (version 25.01.02) toolbox for 'MATLAB'.
License: GPL-2
Encoding: UTF-8
LazyData: true
Depends: R (≥ 4.0.0)
Imports: Matrix, expm, MASS, methods, stats, utils
Suggests: testthat (≥ 3.0.0), knitr, rmarkdown, R.matlab, withr
URL: https://github.com/Kay202/rsDCM
BugReports: https://github.com/Kay202/rsDCM/issues
VignetteBuilder: knitr
Config/testthat/edition: 3
Config/roxygen2/version: 8.0.0
NeedsCompilation: no
Packaged: 2026-09-25 04:09:30 UTC; nsanyal
Author: Godfred Arhin [aut, cre, trl], Nilotpal Sanyal [aut], SPM25 Authors [ctb, cph] (Original MATLAB SPM25 (v25.01.02) implementation; see LICENSE.note)
Repository: CRAN
Date/Publication: 2026-10-06 07:40:02 UTC

rsDCM: Robust and Sparse Dynamic Causal Modelling for Functional MRI

Description

rsDCM provides a robust and sparse method for group-level Dynamic Causal Modelling (DCM) of fMRI data, together with an R implementation of single-subject DCM ported from the MATLAB SPM25 toolbox.

Details

The package's main method is rsdcm (with the convenience wrapper rsdcm_fit): a group-level model that weights subjects by a Student-t likelihood (robustness to outliers) and selects group-level effects with a nonlocal product-moment (pMOM) spike-and-slab prior (sparsity), with ReML-estimated between-subject variance components. The “rs” is robust and sparse; it does not denote resting-state fMRI.

For single-subject inversion, the main entry point is dcm_estimate.

Group-level modelling (robust and sparse)

rsdcm

Robust, sparse group DCM (Student-t + pMOM).

rsdcm_fit

Assemble rsdcm inputs from a list of fitted DCMs.

rsdcm_options

Get or set package numerical options.

Single-subject DCM estimation (SPM25 port)

dcm_estimate

Full DCM inversion (main entry point).

dcm_nlsi_GN

Variational Laplace / Gauss-Newton inversion.

dcm_int

Bilinear-system integrator (forward model).

dcm_fmri_priors

Construct fMRI DCM priors.

dcm_fx_fmri, dcm_gx_fmri

Neural-state and BOLD observation equations.

dcm_bireduce, dcm_kernels

Bilinear reduction and Volterra kernels.

dcm_evidence, dcm_log_evidence, dcm_log_evidence_reduce

Model evidence and Bayesian model reduction.

Parametric Empirical Bayes (Gaussian group modelling)

dcm_peb_prepare, dcm_peb_run

Prepare and fit a second- or third-level PEB.

dcm_peb_of_pebs

Third-level PEB-of-PEBs over a directory of subject PEBs.

dcm_peb_design, dcm_peb_files, dcm_peb_load

Design-matrix and subject-PEB loading helpers.

Datasets

toy_dcm

Three-region single-subject DCM specification.

narps_dcm

48-subject NARPS DCM summaries used by rsdcm.

Lower-level numerical helpers (matrix and vector utilities, numerical differentiation, matrix exponentials, and error-covariance bases) are also exported and documented; see, for example, dcm_vec, dcm_inv, dcm_diff, and dcm_Ce.

Acknowledgment

The single-subject DCM routines are a derivative work of the MATLAB SPM25 (version 25.01.02) toolbox, distributed under GPL-2 by the Wellcome Centre for Human Neuroimaging. See the LICENSE.note file in the package source for the list of ported routines.

Author(s)

Maintainer: Godfred Arhin arhin0122@gmail.com [translator]

Authors:

Other contributors:

References

Arhin, G., Sanyal, N. (2026). Robust and Sparse Group Dynamic Causal Modeling via Student-t Parametric Empirical Bayes and Nonlocal Priors. arXiv:2609.06379. doi:10.48550/arXiv.2609.06379

Botvinik-Nezer, R., Holzmeister, F., Camerer, C.F., et al. (2020). Variability in the analysis of a single neuroimaging dataset by many teams. Nature, 582(7810), 84-88. doi:10.1038/s41586-020-2314-9

Friston, K.J., Harrison, L., Penny, W. (2003). Dynamic causal modelling. NeuroImage, 19(4), 1273-1302. doi:10.1016/S1053-8119(03)00202-7

Friston, K.J., Mattout, J., Trujillo-Barreto, N., Ashburner, J., Penny, W. (2007). Variational free energy and the Laplace approximation. NeuroImage, 34(1), 220-234. doi:10.1016/j.neuroimage.2006.08.035

Valerio, D., Peres, A., Bergstrom, F., Seidel, P., Almeida, J. (2025). Neural and behavioral similarity-driven tuning curves for manipulable objects. Imaging Neuroscience, 3. doi:10.1162/imag_a_00482

See Also

Useful links:

Error-covariance basis (AR or FAST)

Description

Construct a list of covariance basis matrices for the noise model in variational inversion. Mirrors SPM25's spm_Ce.

Usage

dcm_Ce(t, v = NULL, a = NULL)

Arguments

t

Either a numeric vector of session lengths (defaults to AR basis) or a string "ar" / "fast" selecting the basis type.

v

Session lengths (when t is a string).

a

AR coefficient (or sampling interval for FAST).

Value

A list of sparse Matrix objects.

Univariate Normal CDF (SPM-style)

Description

Vectorized Normal CDF with NaN for non-positive variances.

Usage

dcm_Ncdf(x, u = 0, v = 1)

Arguments

x

Quantile.

u

Mean.

v

Variance.

Value

Numeric vector of CDF values.

Examples

dcm_Ncdf(0)                     # 0.5
round(dcm_Ncdf(c(-1.96, 0, 1.96)), 4)

# Note the third argument is a VARIANCE, not a standard deviation
dcm_Ncdf(2, u = 0, v = 4)       # == pnorm(2, mean = 0, sd = 2)

AR(p) autocorrelation / precision matrix

Description

Returns either a banded precision (for q != 0) or a Toeplitz autocorrelation derived from the AR coefficients a.

Usage

dcm_Q(a, n, q = 0)

Arguments

a

AR coefficient vector.

n

Matrix dimension.

q

Output flag: 0 for autocorrelation, otherwise precision.

Value

Sparse Matrix or numeric matrix.

Bilinear reduction of a non-linear DCM

Description

Reduces a non-linear state equation to bilinear form dx/dt = M0 x + sum_i u_i M1[[i]] x, also returning the lead-field expansion. Mirrors SPM25's spm_bireduce.

Usage

dcm_bireduce(M, P)

Arguments

M

Model list (with f, g, x, u, ...).

P

Parameter structure.

Value

List with M0, M1, L1, L2.

Stack a list of equally-shaped block matrices into a column-vectorized matrix

Description

Stack a list of equally-shaped block matrices into a column-vectorized matrix

Usage

dcm_blocks_to_mat(blocks, row_major = FALSE)

Arguments

blocks

List of numeric matrices, all the same shape.

row_major

If TRUE, vectorize by rows.

Value

A matrix whose columns are the vectorized blocks.

Concatenate a list of matrix blocks

Description

Concatenate a list of matrix blocks

Usage

dcm_cat(x, d = NULL)

Arguments

x

A list (or matrix of lists) of numeric blocks.

d

Optional dimension along which to concatenate.

Value

A single matrix.

Compute a stable identifier for a numeric structure

Description

Used by SPM25 to fingerprint a dataset; returned in the DCM result.

Usage

dcm_data_id(...)

Arguments

...

Numeric data (any shape, possibly nested).

Value

Numeric scalar.

Polynomial detrend (column-wise)

Description

Polynomial detrend (column-wise)

Usage

dcm_detrend(x, p = 0)

Arguments

x

Matrix or vector.

p

Polynomial order. p = 0 just centres the columns.

Value

Detrended matrix of the same shape.

Examples

# p = 0 centres each column
x <- cbind(1:10, (1:10) * 2 + 5)
round(colMeans(dcm_detrend(x, 0)), 10)

# p = 1 removes a linear trend, so perfectly linear columns go to ~0
round(dcm_detrend(x, 1), 8)

Block-diagonal layout from a list

Description

Block-diagonal layout from a list

Usage

dcm_diag(...)

Arguments

...

A list and an optional offset K.

Value

A list (or matrix) with the inputs on the diagonal.

High-order numerical Jacobian

Description

Forward-difference numerical Jacobian of a (possibly nested-output) function with respect to one or more of its arguments. Mirrors SPM25's dcm_diff.

Usage

dcm_diff(...)

Arguments

...

Function, its arguments, the index (or vector of indices) of the argument(s) to differentiate with respect to, and optionally a list of projection matrices.

Value

A list containing the Jacobian J (or higher-order derivatives) and the function value f0.

Local linearisation integration step

Description

Computes a one-step Bayesian update dx = (expm(t*J) - I) * J^{-1} * f via an augmented matrix exponential. Mirrors SPM25's spm_dx.

Usage

dcm_dx(dfdx, f, t = Inf, Q = NULL)

Arguments

dfdx

Jacobian df/dx of the model.

f

Current residual / gradient.

t

Step length (or list with a single regulariser).

Q

Optional skew matrix.

Value

Update dx matching the structure of f.

Euclidean column normalization

Description

Euclidean column normalization

Usage

dcm_en(X, p = NULL)

Arguments

X

Matrix.

p

Optional polynomial detrend order applied first.

Value

X with each non-zero column scaled to unit Euclidean norm.

Examples

X <- matrix(c(3, 4, 0, 0, 5, 12), nrow = 2)
dcm_en(X)
sqrt(colSums(dcm_en(X)^2))  # 1, 0, 1 (the all-zero column is left alone)

Estimate a Dynamic Causal Model for fMRI

Description

Inverts a fully-specified DCM using variational Laplace inversion. This is the main user-facing entry point of the package. Mirrors SPM25's dcm_estimate. Progress is reported via message() and can be silenced with suppressMessages().

Usage

dcm_estimate(P, save = FALSE)

Arguments

P

A DCM list, or a path to an .RData/.rds file containing one.

save

Logical. If TRUE and P is a file path, the estimated DCM is saved back to that path.

Value

The estimated DCM (a list) with posterior fields populated. The main fields of interest are Ep (posterior expectations, with Ep$A, Ep$B, Ep$C the connectivity estimates), Cp (posterior covariance), Pp (posterior probabilities), and F (the negative free energy, used for model comparison).

Supported model variants

This release implements only the deterministic, single-state fMRI DCM. The two-state (options$two_state), stochastic (options$stochastic) and spectral / cross-spectral-density (options$induced) variants are not yet implemented: switching any of them on causes dcm_estimate to stop with an informative error. These variants are planned for a future update. Non-linear DCM (a non-empty d array) is supported.

See Also

dcm_nlsi_GN for the underlying inversion, dcm_fmri_priors for the priors it constructs.

Examples

data(toy_dcm)
str(toy_dcm, max.level = 1)


# Full inversion of the bundled three-region model (takes about a minute).
# Progress reporting goes through message(), so it can be silenced.
fit <- suppressMessages(dcm_estimate(toy_dcm))
round(fit$Ep$A, 3)   # posterior connectivity estimates
round(fit$Pp$A, 3)   # posterior probability each connection is non-zero
fit$F                # negative free energy, for model comparison

Approximate model evidence (AIC, BIC)

Description

AIC and BIC penalties for an estimated DCM, plus per-region cost terms. Mirrors SPM25's spm_dcm_evidence.

Usage

dcm_evidence(DCM)

Arguments

DCM

An estimated DCM.

Value

List with per-region cost, AIC penalty, BIC penalty, and overall AIC and BIC.

Scaled-and-squared matrix exponential

Description

SPM25-style matrix exponential used as a fallback to expm.

Usage

dcm_expm(J, x = NULL)

Arguments

J

Square matrix.

x

Optional vector to multiply by expm(J) on the right.

Value

Matrix or vector.

Examples

# expm of a diagonal matrix is just exp() of the diagonal
J <- diag(c(-1, -2))
round(dcm_expm(J), 8)
round(diag(exp(c(-1, -2))), 8)

# Supplying x returns expm(J) %*% x without forming the product yourself
round(dcm_expm(J, c(1, 1)), 8)

Locate or describe entries in a structured parameter

Description

Given a structured parameter object X (with named fields), returns the linear indices of one or more named fields, or the structural location of a given linear index.

Usage

dcm_fieldindices(X, ...)

Arguments

X

A structured parameter (list with numeric fields).

...

One or more field names or numeric indices.

Value

Integer indices or a character description.

Find indices of free parameters in a DCM

Description

Find indices of free parameters in a DCM

Usage

dcm_find_pC(...)

Arguments

...

Either (DCM, fields) or (pC, pE, fields).

Value

List with i (free indices), plus pC, pE, Np.

Recover priors from a DCM container

Description

Recover priors from a DCM container

Usage

dcm_find_rC(DCM)

Arguments

DCM

A DCM list.

Value

List with pC, pE.

DCM mode generator

Description

Maps a vector of mode parameters to a connectivity matrix. Mirrors SPM25's spm_dcm_fmri_mode_gen.

Usage

dcm_fmri_mode_gen(Ev, modes, Cv = NULL)

Arguments

Ev

Numeric vector of mode parameters.

modes

Numeric matrix of mode columns.

Cv

Optional covariance of Ev.

Value

If Cv is NULL, the connectivity matrix Ep. Otherwise a list with Ep and propagated covariance Cp.

Construct fMRI DCM priors

Description

Builds the prior expectations pE and prior covariance pC for a deterministic fMRI DCM. Mirrors SPM25's dcm_fmri_priors.

Usage

dcm_fmri_priors(A, B, C, D, options = list())

Arguments

A

Connectivity adjacency matrix.

B

Modulatory adjacency array.

C

Driving-input adjacency matrix.

D

Non-linear adjacency array.

options

List of model options (stochastic, induced, two_state, backwards, precision, decay).

Value

List with pE (prior expectation), x (initial state template), and pC (prior covariance).

Model variants

Only the deterministic, single-state model is supported end-to-end in this release. Branches for the two-state (options$two_state) and spectral (options$induced) variants exist but are experimental and are not wired through dcm_estimate, which rejects them. They are planned for a future update.

See Also

dcm_estimate, which builds these priors for you.

Examples

# Priors for the bundled three-region (deterministic) model
data(toy_dcm)
pr <- dcm_fmri_priors(toy_dcm$a, toy_dcm$b, toy_dcm$c,
                      D = NULL, options = toy_dcm$options)
names(pr)
pr$pE$A          # prior expectation of the endogenous connections
dim(pr$x)        # 3 regions x 5 hemodynamic states

Coerce a function-like object to a function

Description

Accepts a function, a function name (as character), or a one-line body and returns a real function. Mirrors SPM25's spm_funcheck.

Usage

dcm_funcheck(f)

Arguments

f

A function or character string.

Value

A function.

fMRI neural and hemodynamic state equation

Description

Computes dx/dt for the bilinear neural state equation coupled to the Buxton-Friston hemodynamic model. Mirrors SPM25's spm_fx_fmri.

Usage

dcm_fx_fmri(x, u, P, M)

Arguments

x

State matrix (rows = regions, cols = state variables).

u

Driving inputs at the current time.

P

List of model parameters (A, B, C, D, transit, decay, epsilon, ...).

M

Model structure list (optional).

Value

Matrix dx/dt of the same shape as x.

fMRI state equation with analytic Jacobians

Description

Same as dcm_fx_fmri but additionally returns analytic Jacobians dfdx and dfdu. Mirrors SPM25's spm_fx_fmri when called with nargout > 1.

Usage

dcm_fx_fmri2(x, u, P, M)

Arguments

x

State matrix (rows = regions, cols = state variables).

u

Driving inputs at the current time.

P

List of model parameters (A, B, C, D, transit, decay, epsilon, ...).

M

Model structure list (optional).

Value

List with f, dfdx, D = 1, dfdu.

fMRI BOLD observation equation

Description

Computes the BOLD signal from the hemodynamic state variables. Mirrors SPM25's spm_gx_fmri.

Usage

dcm_gx_fmri(x, u, P, M)

Arguments

x

State matrix (rows = regions, cols = state variables).

u

Driving inputs at the current time.

P

List of model parameters (A, B, C, D, transit, decay, epsilon, ...).

M

Model structure list (optional).

Value

List with the BOLD prediction g and Jacobian dgdx.

Integrate a bilinear DCM

Description

Forward-integrates the model defined by M$f (state equation) and M$g (observation equation) under the input U, returning the predicted observations at the requested sample points.

Usage

dcm_int(P, M, U)

Arguments

P

Parameter structure.

M

Model list. Required fields: f, x; optional g, ns, delays, l.

U

List with input matrix u and time step dt.

Details

M$f and M$g may be given as functions or as the names of functions. M$n must be the length of the flattened state (length(dcm_vec(M$x)), i.e. regions x hidden states), while M$l is the number of observed outputs (regions) and M$m the number of inputs.

Value

Numeric matrix of predicted observations (rows = samples).

See Also

dcm_estimate, which calls this as its forward model.

Examples

# Simulate the BOLD response of a two-region model to a boxcar input.
n <- 2L
pri <- dcm_fmri_priors(A = matrix(1, n, n),
                       B = array(0, c(n, n, 1)),
                       C = matrix(c(1, 0), n, 1),
                       D = array(0, c(n, n, 0)),
                       options = list())

U <- list(u = matrix(c(rep(1, 16), rep(0, 16)), ncol = 1), dt = 1)
M <- list(f = "dcm_fx_fmri", g = "dcm_gx_fmri", x = pri$x,
          m = ncol(U$u), n = length(pri$x), l = nrow(pri$x), ns = 32)

# The priors put C at zero, so start from them and switch on a driving
# input to region 1 and a connection from region 1 to region 2.
P <- pri$pE
P$C[1, 1] <- 1
P$A[2, 1] <- 0.4

y <- dcm_int(P, M, U)
dim(y)                  # 32 samples x 2 regions
round(y[seq(1, 32, 4), ], 3)

Inverse of an ill-conditioned matrix

Description

Computes solve(A + tol*I) with an automatically chosen tolerance, matching SPM25's spm_inv. Includes a fast path for diagonal Matrix objects.

Usage

dcm_inv(A, TOL = NULL)

Arguments

A

Square numeric matrix (dense or Matrix).

TOL

Optional tolerance. If NULL, set automatically.

Value

Inverse matrix.

Examples

A <- matrix(c(2, 1, 1, 2), 2, 2)
dcm_inv(A)

# Unlike solve(), a singular matrix is regularised rather than an error
singular <- matrix(1, 2, 2)
dcm_inv(singular)

Volterra kernels of a bilinear system

Description

Computes the first- and second-order Volterra kernels of a bilinear system specified in M0/M1 form. Mirrors SPM25's spm_kernels.

Usage

dcm_kernels(...)

Arguments

...

Either (M0, M1, N, dt), (M0, M1, L1, N, dt), or (M0, M1, L1, L2, N, dt).

Value

List with kernels K0, K1, K2, and the intermediate response H1.

Total number of numeric entries in a nested structure

Description

Returns the length dcm_vec(X) would produce, without allocating it.

Usage

dcm_length(X)

Arguments

X

A numeric, logical, or list (possibly nested).

Value

Integer length.

Examples

P <- list(A = matrix(1:4, 2), C = c(5, 6))
dcm_length(P)            # 6
length(dcm_vec(P))       # same, but allocates the vector

Bayesian model reduction (full-rank)

Description

Computes the change in log-evidence and the reduced posterior when replacing the original prior with a reduced prior. Mirrors SPM25's dcm_log_evidence.

Usage

dcm_log_evidence(qE, qC, pE, pC, rE = NULL, rC = NULL, ...)

Arguments

qE

Posterior expectation under the original priors.

qC

Posterior covariance.

pE

Original prior expectation.

pC

Original prior covariance.

rE

Reduced prior expectation.

rC

Reduced prior covariance.

...

Passed to rE if it is a function.

Value

List with F, sE, sC.

Bayesian model reduction (subspace projection)

Description

Reduced-rank version of dcm_log_evidence. Mirrors SPM25's dcm_log_evidence_reduce.

Usage

dcm_log_evidence_reduce(qE, qC, pE, pC, rE, rC, TOL = 1e-08)

Arguments

qE

Posterior expectation under the original priors.

qC

Posterior covariance.

pE

Original prior expectation.

pC

Original prior covariance.

rE

Reduced prior expectation.

rC

Reduced prior covariance.

TOL

Tolerance.

Value

List with F, sE, sC.

Log-determinant of a (semi-)definite matrix

Description

Robust log-determinant suitable for sparse, dense, or rank-deficient covariance matrices. Mirrors SPM25's spm_logdet, with fast paths for 1x1 and diagonal inputs.

Usage

dcm_logdet(C)

Arguments

C

Square matrix.

Value

Numeric log-determinant (or NaN if non-positive).

Examples

C <- diag(c(1, 2, 4))
dcm_logdet(C)          # log(1 * 2 * 4)
log(prod(c(1, 2, 4)))  # same

# Zero rows/columns are dropped rather than sending the result to -Inf
dcm_logdet(diag(c(1, 2, 0)))

Variational Laplace inversion of a non-linear system

Description

Performs Gauss-Newton optimisation of the variational free energy for a non-linear forward model with Gaussian priors. Mirrors SPM25's dcm_nlsi_GN.

Usage

dcm_nlsi_GN(M, U, Y)

Arguments

M

Model specification (list with IS/f/g, priors pE, pC, hyperpriors hE, hC, ...).

U

Input structure passed through to M$IS.

Y

Data (or list with $y, $Q, $X0, $dt).

Details

Most users should call dcm_estimate instead, which assembles M, U and Y from a DCM specification and calls this function. Use dcm_nlsi_GN directly only to invert a non-linear model that is not an fMRI DCM.

Progress is reported per Gauss-Newton iteration as EM:(+) k F: ..., where (+) marks an accepted step and (-) a rejected one. Set M$noprint <- 1 to silence it.

Value

List with posterior expectation Ep, covariance Cp, log-precision estimate Eh, free energy F, and components.

See Also

dcm_estimate for the user-facing entry point.

Examples

# dcm_nlsi_GN is the Gauss-Newton optimiser that dcm_estimate() calls after
# assembling M, U and Y from a DCM specification. For a runnable inversion
# see the example in ?dcm_estimate; dcm_nlsi_GN returns the same posterior
# fields (Ep, Cp, Eh, F).

Build a group-level (between-subject) design matrix

Description

Build a group-level (between-subject) design matrix

Usage

dcm_peb_design(n_subj, covariates = NULL, Xnames = NULL)

Arguments

n_subj

Number of subjects (rows).

covariates

Optional data frame of between-subject covariates with n_subj rows. NULL gives an intercept-only design.

Xnames

Optional character vector of column names for the design.

Value

A numeric design matrix with n_subj rows.

Examples

# Intercept only: tests the group mean
dcm_peb_design(4)

# With a between-subject covariate
dcm_peb_design(4, covariates = data.frame(age = c(21, 34, 46, 58)))

Locate subject PEB files in a directory

Description

Finds files named like PEB_sub*.rds or PEB_sub*.mat and extracts the numeric subject identifier from each filename.

Usage

dcm_peb_files(peb_dir)

Arguments

peb_dir

Directory to search.

Value

List with files (full paths) and subjects (integer identifiers), in the order returned by list.files.

Examples

# Set up a directory with two dummy subject PEB files
d <- tempfile(); dir.create(d)
saveRDS(list(), file.path(d, "PEB_sub-01.rds"))
saveRDS(list(), file.path(d, "PEB_sub-02.rds"))

found <- dcm_peb_files(d)
basename(found$files)
found$subjects

unlink(d, recursive = TRUE)

Load subject PEBs from .rds or .mat files

Description

Loads every subject PEB found by dcm_peb_files and normalises MATLAB structs to native R-list shape.

Usage

dcm_peb_load(peb_dir, subjects = NULL)

Arguments

peb_dir

Directory containing the subject PEB files.

subjects

Optional integer vector; keep only these subject numbers.

Details

Reading .mat files requires the R.matlab package. It is listed in Suggests, so install it only if you have MATLAB PEBs to read; .rds input needs nothing extra.

Value

List with pebs (named list of PEB structures), subjects, and files.

See Also

dcm_peb_of_pebs, which calls this.

Run a third-level PEB-of-PEBs over a directory of subject PEBs

Description

Loads subject PEBs (.rds or .mat), builds the third-level design, runs the PEB-of-PEBs, and returns a tidy summary. Optionally saves the group PEB.

Usage

dcm_peb_of_pebs(
  peb_dir,
  subjects = NULL,
  covariates = NULL,
  Xnames = NULL,
  save_path = NULL,
  M = list(Q = "single", alpha = 1, beta = 16, maxit = 64),
  field = "all",
  verbose = TRUE
)

Arguments

peb_dir

Directory containing the subject PEB files.

subjects

Optional integer vector; restrict to these subject numbers.

covariates

Optional data frame of between-subject covariates, one row per subject in the order returned by dcm_peb_files. NULL gives an intercept-only design (the group mean).

Xnames

Optional custom design column names.

save_path

Optional path to save the group PEB to, as .rds. NULL (the default) does not write anything to disk.

M

Third-level model options passed to dcm_peb_prepare.

field

Parameter blocks to model; "all" is typical for PEB inputs.

verbose

Logical. Report progress via message().

Details

Nothing is written to disk unless save_path is supplied.

Value

List with group_PEB, save_path (NULL if not saved), subjects, X (the design), and summary (a tidy data frame of parameter estimates).

See Also

dcm_peb_prepare, dcm_peb_run.

Prepare a PEB (second- or third-level) model

Description

Gathers the first-level posterior densities, selects the parameter indices implied by field, projects them onto a rank-reduced subspace, and builds the second-level priors, hyperpriors and precision components. Mirrors the preparation half of SPM25's spm_dcm_peb.

Usage

dcm_peb_prepare(P, M = list(), field = c("A", "B"))

Arguments

P

List of estimated DCMs (or, for a PEB-of-PEBs, a list of PEBs). Each element must carry M$pE, M$pC, Ep, Cp and F.

M

Second-level model specification. Recognised fields: X (design matrix, subjects x covariates), W, Q (covariance component option, see Details), bE, bC, pC, hE, hC, alpha, beta, Xnames, maxit. Defaults to an intercept-only design.

field

Character vector of parameter blocks to model (e.g. c("A", "B")), the string "all" for every named block, or numeric indices into dcm_vec(M$pE).

Details

M$Q selects the form of the second-level precision components: "single" (one component, the default), "fields" (one per requested field), "all" (one per parameter), "none" (no components), or a list of numeric matrices for manual specification.

Value

A list ("prep") consumed by dcm_peb_run, holding the projected densities, design, priors, precision components and sizes.

See Also

dcm_peb_run to fit the prepared model.

Fit a prepared PEB model by variational Laplace

Description

Runs the variational Laplace (Fisher scoring) loop on the structure built by dcm_peb_prepare, and returns the group-level PEB together with the first-level DCMs updated under the empirical priors. Mirrors the estimation half of SPM25's spm_dcm_peb.

Usage

dcm_peb_run(prep, verbose = TRUE)

Arguments

prep

The list returned by dcm_peb_prepare.

verbose

Logical. Report free energy per iteration via message(). Silence with suppressMessages() or verbose = FALSE.

Details

Progress is reported as VL Iter k: F=... dF=... [t=...], where t is the log step size of the Fisher scoring update. The loop stops when the step size collapses or the free-energy increase falls below 1e-4.

Value

List with two elements: PEB, the group-level result (Ep group parameter expectations, Cp their covariance, Eh/Ch log-precision estimates, F free energy, and the Pnames/Xnames/Snames labels), and P, the input DCMs with Ep, Cp, M$pE, M$pC and F updated under the empirical priors.

See Also

dcm_peb_prepare, dcm_peb_of_pebs.

Logistic sigmoid

Description

Standard SPM logistic function 1/(1 + exp(-x)), vectorized.

Usage

dcm_phi(x)

Arguments

x

Numeric input.

Value

Numeric output of the same shape.

Pseudo-inverse via SVD with truncation

Description

Pseudo-inverse via SVD with truncation

Usage

dcm_pinv(A, TOL = NULL)

Arguments

A

Numeric matrix.

TOL

Singular-value tolerance.

Value

Pseudo-inverse of A.

Examples

# Left inverse of a tall matrix
A <- matrix(c(1, 2, 3, 4, 5, 7), nrow = 3)
round(dcm_pinv(A) %*% A, 8)   # ~ 2 x 2 identity

Sparse identity-like matrix

Description

Build an identity-like sparse matrix of any shape, with optional shifted diagonals. Mirrors SPM25's spm_speye.

Usage

dcm_speye(m, n = m, k = 0, c = 0)

Arguments

m, n

Output dimensions.

k

Diagonal offset (0 for main).

c

Cyclic-fill flag (0/1/2) as in SPM25.

Value

A sparse Matrix.

Examples

dcm_speye(3)           # 3 x 3 sparse identity
dcm_speye(3, 3, 1)     # ones on the first super-diagonal
dcm_speye(2, 4)        # non-square is fine

Truncated SVD with sparsity awareness

Description

Computes a thin SVD and drops singular components below a relative tolerance U. Mirrors SPM25's spm_svd.

Usage

dcm_svd(X, U = NULL)

Arguments

X

Numeric matrix.

U

Relative tolerance for retaining singular values.

Value

List with U, S, V.

Trace of a matrix product

Description

Computes trace(A %*% B) efficiently as sum(t(A) * B).

Usage

dcm_trace(A, B)

Arguments

A, B

Numeric matrices.

Value

Numeric scalar.

Reshape a flat vector back into a template structure

Description

Inverse of dcm_vec: splits a flat numeric vector across the shape implied by one or more template arguments.

Usage

dcm_unvec(vX, ...)

Arguments

vX

Flat numeric vector.

...

One or more templates whose structure is copied; vX is distributed across them in order.

Value

A structure (or list of structures) matching the templates.

See Also

dcm_vec for the forward operation.

Examples

# Round-trip: flatten a structure and rebuild it
P <- list(A = matrix(1:4, 2), C = c(5, 6))
identical(dcm_unvec(dcm_vec(P), P), P)

# The template supplies the shape; the vector supplies the values
tmpl <- list(A = matrix(0, 2, 2), C = numeric(2))
dcm_unvec(1:6, tmpl)

Flatten a nested numeric structure into a vector

Description

Recursively walks a numeric, logical, or nested list structure and returns a single flat numeric vector. Internal helper; mirrors SPM25's spm_vec.

Usage

dcm_vec(X, ...)

Arguments

X

A numeric, logical, or list (possibly nested).

...

Additional structures to flatten and concatenate.

Value

A numeric vector.

See Also

dcm_unvec for the inverse operation.

Examples

# A DCM parameter structure flattens in field order
P <- list(A = matrix(1:4, 2), C = c(5, 6))
dcm_vec(P)

# Nested lists are walked recursively
dcm_vec(list(a = 1, b = list(c = 2:3, d = 4)))

Zero-fill a structure

Description

Returns a copy of X with the same shape but every numeric entry replaced by 0.

Usage

dcm_zeros(X)

Arguments

X

A nested structure.

Value

A structure with zeros, same shape as X.

Examples

dcm_zeros(list(A = matrix(1:4, 2), C = c(5, 6)))

Krylov-based action of matrix exponential

Description

Computes w = expm(t*A) %*% v without forming expm(t*A). Used as a fallback in dcm_dx for very large systems.

Usage

expv(t, A, v, tol = 1e-07, m = NULL)

Arguments

t

Step length.

A

Square matrix.

v

Vector.

tol

Tolerance.

m

Krylov subspace dimension.

Value

Vector expm(t*A) %*% v.

Subject-level DCM summaries for a group analysis (NARPS)

Description

Subject-level DCM posterior summaries for 48 subjects, used to demonstrate the robust and sparse group-level model rsdcm. The summaries are the inputs a group analysis needs: a posterior mean and covariance per subject.

Usage

narps_dcm

Format

A list with:

eta

48 x 22 matrix of subject-level posterior means (subjects in rows, parameters in columns).

Cp

Length-48 list of 22 x 22 posterior covariance matrices.

parameter_names

Length-22 character vector, e.g. "A(1,1)", "B(2,1,1)".

subject

Length-48 character vector of de-identified subject ids.

covariates

Data frame of between-subject covariates: group, gender, age.

regions

The four region labels, in order.

Details

Each subject has 22 parameters: the 16 intrinsic connections of a four-region DCM (A) plus 6 task-modulatory connections (B), across the regions vmPFC, vStr (ventral striatum), amygdala, and anterior insula.

Source

Subject-level Dynamic Causal Modelling summaries computed by the package author from the openly shared NARPS dataset: Botvinik-Nezer, R., Holzmeister, F., Camerer, C. F., et al. (2020), "Variability in the analysis of a single neuroimaging dataset by many teams", Nature, 582(7810), 84-88, doi:10.1038/s41586-020-2314-9. The build script that prepared the shipped object is in data-raw/make-narps-dcm.R.

See Also

rsdcm, which this dataset is the example input for.

Examples

data(narps_dcm)
dim(narps_dcm$eta)               # 48 subjects x 22 parameters
length(narps_dcm$Cp)             # one posterior covariance per subject
table(narps_dcm$covariates$group)

Pade approximation of the matrix exponential

Description

Pure-R Pade approximation, used as a fallback. See also expm.

Usage

padm(A, p = 6)

Arguments

A

Square matrix.

p

Order.

Value

Matrix exponential of A.

Robust and sparse group-level DCM (Student-t + pMOM)

Description

Fits a group-level model to subject-level DCM parameter estimates that is robust to outlier subjects and performs Bayesian variable selection on the group effects. The model combines three ingredients: Student-t weighting of subjects (robustness), a nonlocal product-moment (pMOM) spike-and-slab prior on the group effects (sparsity / inclusion probabilities), and ReML-estimated between-subject variance components. Estimation is by EM.

Usage

rsdcm(
  eta_theta_y,
  C_theta_y_list,
  X_G,
  V_list,
  nu = 3,
  tau0 = 0.05,
  tau1 = 1,
  pi = 0.5,
  a0 = 2,
  b0 = 0.01,
  a1 = 2,
  b1 = 1,
  min_slab_spike_ratio = 1,
  max_iter = 500,
  tol = 1e-06,
  inner_sweeps = 3,
  verbose = TRUE,
  max_beta_step = 1,
  min_alpha = 1e-08,
  max_alpha = 100
)

Arguments

eta_theta_y

N x p matrix of subject-level posterior means (subjects in rows, parameters in columns).

C_theta_y_list

Length-N list of p x p subject-level posterior covariance matrices.

X_G

N x r group (between-subject) design matrix.

V_list

List of p x p basis matrices for the between-subject covariance Sigma_b = sum_k alpha_k V_k. A per-parameter diagonal basis is a common default.

nu

Student-t degrees of freedom (smaller = heavier-tailed, more robust).

tau0, tau1

Initial spike and slab standard deviations.

pi

Prior inclusion probability (slab weight).

a0, b0, a1, b1

Inverse-Gamma hyperprior parameters on tau0^2 and tau1^2.

min_slab_spike_ratio

Identifiability guard: enforce tau1 >= min_slab_spike_ratio * tau0.

max_iter

Maximum EM iterations.

tol

Convergence tolerance on the change in beta and alpha.

inner_sweeps

Coordinate-Newton sweeps per M-step for beta.

verbose

Logical. Report progress via message(); silence with suppressMessages() or verbose = FALSE.

max_beta_step

Per-coordinate cap on the beta Newton step.

min_alpha, max_alpha

Bounds on the variance components.

Details

This is an alternative to the Gaussian Parametric Empirical Bayes layer (dcm_peb_run); the numerics are original to the package, not a port of SPM. Most users will call the convenience wrapper rsdcm_fit, which assembles the arguments below from a list of fitted DCMs.

Value

A list with the group-level effects beta_mat (p x r) and beta_vec, posterior inclusion probabilities inclusion (p x r), Student-t weights, variance components alpha, learned scales tau0/tau1, residual variance sigma2, fitted means mu_hat (N x p), and the sizes N, p, r.

References

Arhin, G., Sanyal, N. (2026). Robust and sparse group dynamic causal modeling via Student-t parametric empirical Bayes and nonlocal priors. arXiv:2609.06379. doi:10.48550/arXiv.2609.06379

See Also

rsdcm_fit for the wrapper over fitted DCMs, dcm_peb_run for the Gaussian PEB alternative, narps_dcm for the example dataset.

Examples

# Real 48-subject NARPS group analysis. The inputs are precomputed
# subject-level DCM summaries (a posterior mean and covariance per subject).
data(narps_dcm)
N <- nrow(narps_dcm$eta)
p <- ncol(narps_dcm$eta)
# Per-parameter variance-component basis, and an intercept-only design.
V_list <- lapply(seq_len(p), function(k) { V <- matrix(0, p, p); V[k, k] <- 1; V })
X_G <- matrix(1, N, 1, dimnames = list(NULL, "intercept"))


fit <- rsdcm(narps_dcm$eta, narps_dcm$Cp, X_G, V_list, verbose = FALSE)
rownames(fit$beta_mat) <- rownames(fit$inclusion) <- narps_dcm$parameter_names

# Group-level connections selected by the spike-and-slab prior (PIP > 0.5)
sel <- fit$inclusion[, 1] > 0.5
round(cbind(estimate = fit$beta_mat[sel, 1],
            PIP      = fit$inclusion[sel, 1]), 3)

# Subjects the Student-t weighting down-weights most
w <- rowMeans(matrix(fit$weights, N, p, byrow = TRUE))
narps_dcm$subject[order(w)][1:5]

Robust sparse group DCM from a list of fitted DCMs

Description

Convenience wrapper around rsdcm that assembles its inputs from a list of estimated DCMs (as returned by dcm_estimate). For each subject it takes the posterior mean Ep and covariance Cp restricted to the requested parameter field, builds a between-subject design, and fits the robust + sparse group model.

Usage

rsdcm_fit(P, field = "A", covariates = NULL, X_G = NULL, V_list = NULL, ...)

Arguments

P

List of estimated DCMs; each must carry Ep, Cp, and M$pE/M$pC (or options) so the field indices can be resolved.

field

Parameter block(s) to model at the group level, named as in the DCM parameter structure Ep. For an fMRI DCM these are "A", "B", "C", "D", "transit", "decay" and "epsilon"; pass one (e.g. "A") or several (e.g. c("A", "B")). Case-sensitive: the connectivity blocks are uppercase ("A", not "a"). Passed to dcm_find_pC.

covariates

Optional data frame of between-subject covariates, one row per subject; NULL gives an intercept-only design (the group mean). Ignored if X_G is supplied.

X_G

Optional group design matrix (N x r); overrides covariates.

V_list

Optional list of p x p variance-component bases; defaults to a per-parameter diagonal basis.

...

Further arguments passed to rsdcm (e.g. nu, pi, verbose).

Details

All subjects must share the same parameterisation, so the selected field indices are required to match across P.

Value

The list returned by rsdcm, with the group effects and inclusion probabilities labelled by the selected parameters and design columns, plus the selected param_index.

See Also

rsdcm, dcm_estimate, dcm_peb_design.

Examples

# rsdcm_fit() consumes a list of fitted DCMs (each from dcm_estimate), then
# fits the group model on a chosen field:
#   fits <- lapply(dcm_files, function(f) dcm_estimate(readRDS(f)))
#   grp  <- rsdcm_fit(fits, field = "A", covariates = my_covariates)
# For a runnable group-level example on precomputed subject summaries, see
# ?rsdcm and the narps_dcm dataset.

Get or set rsDCM runtime options

Description

Read or modify package-internal numerical options. Currently the only option is GLOBAL_DX, the finite-difference step used by dcm_diff when computing numerical Jacobians.

Usage

rsdcm_options(...)

Arguments

...

Named arguments to set. Call with no arguments to retrieve the current option list.

Value

The current option list (invisibly when setting).

Examples

rsdcm_options()                 # show current options
rsdcm_options(GLOBAL_DX = 1e-4) # widen the FD step
rsdcm_options(GLOBAL_DX = exp(-8)) # restore default

Example DCM for fMRI inversion

Description

A three-region Dynamic Causal Model used in the documentation, vignette, and tests. It is an input specification (no posterior fields), ready to invert with dcm_estimate.

Usage

toy_dcm

Format

A list with the standard SPM-style DCM fields:

a

3 x 3 binary matrix of endogenous connections to estimate.

b

3 x 3 x 1 binary array of modulatory connections.

c

3 x 1 binary matrix of driving inputs; here the input enters region 3 only.

Y

Observed BOLD data. Y$y is the 482 x 3 time-series matrix (rows = scans, columns = regions); Y$dt is the sampling interval in seconds; Y$X0 is the confound design; Y$Q is the error-covariance component list; Y$name holds the region names.

U

Input design. U$u is the 482 x 1 input matrix at the microtime resolution; U$dt is the microtime step; U$name holds the input names.

n, v

Number of regions (3) and number of scans (482).

TE

Echo time in seconds.

name

Name of the model.

xY

Region (VOI) structure carried over from the source DCM.

options

Estimation options: nonlinear, two_state, stochastic, centre, mode, maxit, induced, maxnodes, hE, hC.

Details

The model has 3 regions and 482 scans, with a single driving input entering region 3 and a single modulatory input. All entries of the a matrix are enabled, so every directed connection between the three regions is estimated.

Source

Derived from the fMRI data of Valerio et al. (2025), "Neural and behavioral similarity-driven tuning curves for manipulable objects", Imaging Neuroscience, 3, doi:10.1162/imag_a_00482, via a Dynamic Causal Modelling exercise carried out with that data during package development. The build script that prepared the shipped object is in data-raw/make-toy-dcm.R. See the package CITATION for the full author list.

Examples

data(toy_dcm)
str(toy_dcm, max.level = 1)
# Invert it with dcm_estimate(toy_dcm); see ?dcm_estimate for a runnable
# (donttest) inversion example.