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.
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
)Numeric N by p matrix of subject-level
posterior means.
Length-N list of p by p
posterior covariance matrices.
Numeric N by r between-subject design matrix.
List of symmetric p by p bases defining
\(\Sigma_b=\sum_k\alpha_k V_k\). The matrices must be linearly
independent.
Positive Student-t degrees of freedom.
Initial spike and slab standard deviations.
Prior slab probability in (0, 1).
Inverse-Gamma hyperparameters for the squared spike and slab scales.
Minimum allowed ratio tau1/tau0.
Maximum generalized-EM iterations per start.
Relative parameter and objective convergence tolerance.
Number of coordinate-Newton sweeps in each conditional maximization of the group effects.
Logical. Print fitting and MCMC diagnostics.
Maximum absolute coordinate-Newton step.
Bounds for covariance-component amplitudes.
Pre-whitening factorization. One of
"precision_chol", "covariance_chol", or
"symmetric".
Logical. Estimate the global residual scale. The default fixes it to avoid confounding with covariance-component scale.
Positive residual scale used when
estimate_sigma2 = FALSE.
Optional initial coefficient and covariance- component vectors.
Number of empirical-Bayes starts.
Standard deviation multiplier for randomized starts.
Optional random seed for randomized starts.
Numerical tolerance for monotonicity safeguards.
Maximum coefficient-score norm required for convergence.
Maximum backtracking halvings for a coefficient update.
Logical. Compute integrated PIPs for every start rather than only the selected start.
Parallel-tempering MCMC iteration, burn-in, and thinning settings.
Number of independent temperature ladders.
Number of replicas in each ladder.
Inverse temperature of the hottest replica.
Number of iterations between replica swaps.
Probability of a slab-scale proposal.
Optional random seed for posterior PIP integration.
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.
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.
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.
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
})
# \donttest{
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])
#> intercept PIP
#> P1 0.764428748 1.000000000
#> P2 0.011277761 0.005746535
#> P3 -0.005717361 0.001908935
# }