Skip to contents

Simulation 3: arbitrary covariance and choosing K

Source: vignettes/simulation3.Rmd
simulation3.Rmd

Simulation 3 (paper §4.3) stress-tests the model when the true covariance is arbitrary (not low-rank) and varies with a binary covariate. It also illustrates choosing the working dimension K and checking robustness to it. The paper uses S=25S = 25, N=50N = 50, J=100J = 100 and K=25K = 25; a smaller design is used here.

library(BCAIA)
set.seed(2)
s <- 15; n <- 2 * s; J <- 40
subject <- rep(seq_len(s), 2)
xd <- rep(c(0, 1), each = s)
Xcov  <- cbind(1, xd)     # intercept + binary covariate
Xmean <- cbind(xd)

## two arbitrary sparse covariance matrices, one per level of xd
rand_cov <- function(J, sparsity = 0.85) {
  A <- matrix(rnorm(J * J), J, J); A[abs(A) < sparsity * 2] <- 0
  S <- tcrossprod(A) / J + diag(1, J); S
}
Sig0 <- rand_cov(J); Sig1 <- rand_cov(J)
alpha_s <- matrix(rnorm(s * J, 3, 1), s, J)
Y <- matrix(0, n, J)
for (i in seq_len(n)) {
  Sig <- if (xd[i] == 0) Sig0 else Sig1
  L <- chol(Sig + diag(1e-8, J))
  Y[i, ] <- floor(exp(alpha_s[subject[i], ] + as.numeric(crossprod(L, rnorm(J)))))
}

Choosing K from a clr-PCA scree rule 

choose_K(Y, plot = TRUE)

#> [1] 19
#> attr(,"eigenvalues")
#>  [1] 1.382640e+01 9.324728e+00 8.692113e+00 7.118524e+00 6.762344e+00
#>  [6] 6.343753e+00 5.857729e+00 4.942973e+00 3.528089e+00 3.357174e+00
#> [11] 3.130019e+00 2.753344e+00 2.438689e+00 2.316600e+00 2.226240e+00
#> [16] 1.841705e+00 1.489687e+00 1.275955e+00 1.110679e+00 1.073295e+00
#> [21] 8.547790e-01 6.638565e-01 5.032603e-01 4.920572e-01 3.355124e-01
#> [26] 2.869716e-01 1.690666e-01 1.002683e-01 1.541650e-02 2.133646e-15
#> [31] 1.435907e-15 9.652194e-16 6.969157e-16 3.531268e-16 7.975653e-17
#> [36] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00

Fitting 

fit <- bcaia(Y, Xmean, Xcov, subject = subject, K = 25,
             niter = 160000, burnin = 80000, thin = 10, seed = 1)
#> BCAIA fit (subject model)
#>   data: n = 30, J = 40, Pmean = 1, Pcov = 2, K = 10
#>   subjects: s = 15
#>   MCMC: 1000 iters, burn-in 500, thin 5 -> 100 saved samples
#>   runtime: 0.5 min

Group-specific covariance estimates 

S0 <- posterior_cor(fit, c(1, 0))
S1 <- posterior_cor(fit, c(1, 1))
op <- par(mfrow = c(1, 2), mar = c(2, 2, 2, 1))
image(S0, main = "rho-hat (xd = 0)", axes = FALSE)
image(S1, main = "rho-hat (xd = 1)", axes = FALSE)

par(op)

Even though the model assumes a low-rank covariance regression, with a sufficiently large K and the full chain it approximates the arbitrary group-specific covariance structures well, and the posterior estimates are robust to the choice of K (paper §4.3 and Supp. Fig. 7).