Package {rsDCM}

rsDCM: Robust and Sparse Dynamic Causal Modelling for Functional MRI

Description

The package provides a robust and sparse group-level extension of Parametric Empirical Bayes for Dynamic Causal Modelling (DCM), together with an R port of deterministic single-state fMRI DCM routines from SPM25.

Details

The proposed group method combines a Student-t likelihood on pre-whitened first-level posterior summaries with a Gaussian-spike and nonlocal product-moment (pMOM)-slab prior. Mode fitting uses a multistart generalized EM algorithm with conditional maximization. Posterior inclusion probabilities are integrated conditionally on fitted empirical-Bayes hyperparameters by parallel-tempering MCMC. Conditional Laplace standard deviations are also returned for the group coefficients.

Main functions

rsdcm

Fit the robust and sparse group model.

rsdcm_fit

Assemble group inputs from fitted DCMs.

dcm_estimate

Invert a single-subject fMRI DCM.

dcm_peb_run

Fit the Gaussian SPM-style PEB model.

Other functions

rsdcm_options

Get or set package numerical options.

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.

dcm_peb_prepare, dcm_peb_run

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

dcm_peb_design, dcm_peb_files, dcm_peb_load

Design-matrix and subject-PEB loading helpers.

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.

Datasets

toy_dcm

Three-region single-subject DCM specification.

narps_dcm

Corrected 48-subject NARPS DCM summaries.

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:

• Godfred Arhin arhin0122@gmail.com [translator]

• Nilotpal Sanyal nsanyal@utep.edu

Other contributors:

• SPM25 Authors (Original MATLAB SPM25 (v25.01.02) implementation; see LICENSE.note) [contributor, copyright holder]

References

Arhin, G. and 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

Friston, K. J., Harrison, L., and 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

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

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:

• https://github.com/Kay202/rsDCM

• Report bugs at https://github.com/Kay202/rsDCM/issues

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

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

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

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

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

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

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

Value

Numeric scalar.

Polynomial detrend (column-wise)

Description

Polynomial detrend (column-wise)

Usage

dcm_detrend(x, p = 0)

Arguments

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

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

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

Value

Update dx matching the structure of f.

Euclidean column normalization

Description

Euclidean column normalization

Usage

dcm_en(X, p = NULL)

Arguments

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Value

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

Corrected subject-level DCM summaries for NARPS

Description

Corrected within-subject PEB summaries for 48 participants in the openly shared mixed-gambles NARPS data. These summaries replace the preliminary package data used before the first-level input correction.

Usage

narps_dcm

Format

A list with:

eta

A 48 by 22 matrix of posterior means.

Cp

A named length-48 list of 22 by 22 posterior covariances.

parameter_names

Names of the 22 parameters.

subject

De-identified participant identifiers.

covariates

Task version, gender, and age for each participant.

regions

The four region labels in model order.

Details

Each participant has 22 parameters: 16 intrinsic connections from a four-region DCM and six task-modulatory connections. The regions are vmPFC, ventral striatum, amygdala, and anterior insula.

Source

Corrected DCM summaries derived from the NARPS data described by Botvinik-Nezer et al. (2020), Nature 582, 84–88. doi:10.1038/s41586-020-2314-9

Examples

data(narps_dcm)
dim(narps_dcm$eta)
length(narps_dcm$Cp)
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

Value

Matrix exponential of A.

Robust and sparse group-level DCM

Description

Fits a group model to subject-level DCM posterior means and covariances. The likelihood is Student-t after block-covariance pre-whitening. Group coefficients receive a Gaussian-spike and nonlocal pMOM-slab prior.

Usage

rsdcm(
  eta_theta_y,
  C_theta_y_list,
  X_G,
  V_list,
  nu = 3,
  tau0 = 0.05,
  tau1 = 1,
  pi = 0.4,
  a0 = 2,
  b0 = 0.01,
  a1 = 2,
  b1 = 1,
  min_slab_spike_ratio = 1,
  max_iter = 200,
  tol = 1e-06,
  inner_sweeps = 3,
  verbose = TRUE,
  max_beta_step = 1,
  min_alpha = 1e-08,
  max_alpha = 100,
  whitening_method = c("precision_chol", "covariance_chol", "symmetric"),
  estimate_sigma2 = FALSE,
  sigma2_fixed = 1,
  beta_init = NULL,
  alpha_init = NULL,
  n_starts = 3L,
  start_jitter = 0.25,
  start_seed = NULL,
  objective_tol = 1e-10,
  gradient_tol = 1e-04,
  max_backtrack = 20L,
  pip_all_starts = FALSE,
  pip_mcmc_iter = 10000L,
  pip_mcmc_burn = 3000L,
  pip_mcmc_thin = 2L,
  pip_mcmc_chains = 4L,
  pip_mcmc_temperatures = 8L,
  pip_mcmc_min_inverse_temp = 0.05,
  pip_mcmc_swap_every = 2L,
  pip_mcmc_large_jump_prob = 0.1,
  pip_mcmc_seed = NULL
)

Arguments

Details

Mode fitting uses multistart generalized expectation-maximization with conditional maximization. Integrated posterior inclusion probabilities (PIPs) are then computed by parallel-tempering MCMC while holding the empirical-Bayes hyperparameters fixed. The returned coefficient standard deviations use a conditional Laplace approximation at the selected mode.

The Student-t weights apply to residual coordinates defined by the selected whitening factor and the supplied parameter ordering. They are not invariant to arbitrary rotations of the whitened coordinates.

Value

A list containing empirical-Bayes modes beta_mat, conditional posterior means beta_posterior_mean, conditional Laplace standard deviations beta_sd_approx, integrated PIPs inclusion, PIP Monte Carlo errors and chain diagnostics, Student-t weights, fitted scales and covariance components, convergence traces, multistart diagnostics, whitening information, and fitted subject-level means.

References

Johnson, V. E. and Rossell, D. (2010). On the use of non-local prior densities in Bayesian hypothesis tests. Journal of the Royal Statistical Society: Series B, 72(2), 143–170.

See Also

rsdcm_fit, dcm_peb_run

Examples

set.seed(1)
N <- 12L
p <- 3L
X_G <- matrix(1, N, 1, dimnames = list(NULL, "intercept"))
eta <- matrix(rnorm(N * p, sd = 0.25), N, p)
eta[, 1] <- eta[, 1] + 0.7
colnames(eta) <- paste0("P", seq_len(p))
Cp <- replicate(N, diag(0.1, p), simplify = FALSE)
V <- lapply(seq_len(p), function(j) {
  z <- matrix(0, p, p)
  z[j, j] <- 1
  z
})

fit <- rsdcm(
  eta, Cp, X_G, V, n_starts = 1L, max_iter = 100L,
  pip_mcmc_iter = 1000L, pip_mcmc_burn = 400L,
  pip_mcmc_thin = 5L, pip_mcmc_chains = 2L,
  pip_mcmc_temperatures = 4L, verbose = FALSE
)
cbind(fit$beta_mat, PIP = fit$inclusion[, 1])

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

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, conditional Laplace standard deviations, integrated PIPs, and diagnostics labelled by parameter and design columns, plus 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)
# See ?rsdcm for a self-contained group-level example.

Get or set rsDCM runtime options

Description

Read or modify package-internal numerical options. The option GLOBAL_DX controls the finite-difference step used by dcm_diff.

Usage

rsdcm_options(...)

Arguments

Value

The current option list, invisibly when setting values.

Examples

old <- rsdcm_options()$GLOBAL_DX
rsdcm_options(GLOBAL_DX = 1e-4)
rsdcm_options(GLOBAL_DX = old)

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.