Simulation 1: covariate-varying interactions with a categorical design • BCAIA

This vignette reproduces the design of Simulation 1 (paper §4.1): a categorical design with a binary and a ternary covariate giving six experimental conditions, J = 15 features and N = 30 samples. The true covariance is rank-2 and varies with the covariates.

library(BCAIA)
set.seed(6)
sim <- simulate_bcaia(n = 30, J = 15, K_true = 2, seed = 6)
dim(sim$Y)
#> [1] 30 15
round(mean(sim$Y == 0), 3)   # zero rate
#> [1] 0.133

Choosing K

K is set to a value large enough to capture the covariance; redundant factors are shrunk out by the Dirichlet–Horseshoe prior. A data-driven starting point is the clr-PCA scree rule:

choose_K(sim$Y)
#> [1] 5
#> attr(,"eigenvalues")
#>  [1] 2.175147e+01 1.264311e+01 9.413063e+00 5.432482e+00 1.437905e+00
#>  [6] 5.794548e-01 3.896298e-01 3.547540e-01 2.226568e-01 2.022696e-01
#> [11] 1.639294e-01 1.068431e-01 4.800173e-02 3.954617e-02 3.989187e-16

The paper uses K = 5 here (a small over-specification is harmless).

Fitting the model

For illustration we run a short chain. The paper uses niter = 160000; increase niter for real analyses (best run on a server).

fit <- bcaia(sim$Y, sim$Xmean, sim$Xcov, K = 5,
             niter = 160000, burnin = 80000, thin = 10, seed = 1)

#> BCAIA fit (simple model)
#>   data: n = 30, J = 15, Pmean = 5, Pcov = 4, K = 5
#>   MCMC: 1500 iters, burn-in 750, thin 5 -> 150 saved samples
#>   runtime: 0.2 min

Posterior covariance vs. truth

We compare the posterior median covariance to the truth at two conditions.

x1  <- sim$Xcov[1, ]    # condition 1
x26 <- sim$Xcov[26, ]   # condition 6
S1  <- posterior_Sigma(fit, x1)
S26 <- posterior_Sigma(fit, x26)

ut <- upper.tri(S1, diag = TRUE)
c(cond1 = cor(S1[ut],  sim$truth$Sigma[, , 1][ut]),
  cond6 = cor(S26[ut], sim$truth$Sigma[, , 26][ut]))
#>     cond1     cond6 
#> 0.7970499 0.7194578

op <- par(mfrow = c(1, 2), mar = c(2, 2, 2, 1))
image(S1,  main = "Sigma-hat (condition 1)",  axes = FALSE)
image(S26, main = "Sigma-hat (condition 6)",  axes = FALSE)

par(op)

With the full niter = 160000 chain, the posterior estimates closely match the true covariate-varying covariance across all six conditions (paper Fig. 1).