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

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.

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.

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.

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.

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

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

# \donttest{
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)
#>        estimate   PIP
#> A(1,1)   -0.292 0.991
#> A(2,2)   -0.298 0.993
#> A(3,3)   -0.208 0.670

# Subjects the Student-t weighting down-weights most
w <- rowMeans(matrix(fit$weights, N, p, byrow = TRUE))
narps_dcm$subject[order(w)][1:5]
#> [1] "sub-002" "sub-003" "sub-013" "sub-026" "sub-040"
# }