inst/doc/BayesSUR-RE.R

In BayesSUR: Bayesian Seemingly Unrelated Regression Models in High-Dimensional Settings

## ----setup, include=FALSE-----------------------------------------------------
knitr::opts_chunk$set(echo = TRUE, eval = FALSE)
options(rmarkdown.html_vignette.check_title = FALSE)

## ----eval=TRUE----------------------------------------------------------------
library("BayesSUR")

## -----------------------------------------------------------------------------
# library("gRbase")
# sim.ssur <- function(n, s, p, t0 = 0, seed = 123, mv = TRUE,
#                      t.df = Inf, random.intercept = 0, intercept = TRUE) {
#   # set seed to fix coefficients
#   set.seed(7193)
#   sd_b <- 1
#   mu_b <- 1
#   b <- matrix(rnorm((p + ifelse(t0 == 0, 1, 0)) * s, mu_b, sd_b), p + ifelse(t0 == 0, 1, 0), s)
# 
#   # design groups and pathways of Gamma matrix
#   gamma <- matrix(FALSE, p + ifelse(t0 == 0, 1, 0), s)
#   if (t0 == 0) gamma[1, ] <- TRUE
#   gamma[2:6 - ifelse(t0 == 0, 0, 1), 1:5] <- TRUE
#   gamma[11:21 - ifelse(t0 == 0, 0, 1), 6:12] <- TRUE
#   gamma[31:51 - ifelse(t0 == 0, 0, 1), 1:5] <- TRUE
#   gamma[31:51 - ifelse(t0 == 0, 0, 1), 13:15] <- TRUE
#   gamma[52:61 - ifelse(t0 == 0, 0, 1), 1:12] <- TRUE
#   gamma[71:91 - ifelse(t0 == 0, 0, 1), 6:15] <- TRUE
#   gamma[111:121 - ifelse(t0 == 0, 0, 1), 1:15] <- TRUE
#   gamma[122 - ifelse(t0 == 0, 0, 1), 16:18] <- TRUE
#   gamma[123 - ifelse(t0 == 0, 0, 1), 19] <- TRUE
#   gamma[124 - ifelse(t0 == 0, 0, 1), 20] <- TRUE
# 
#   G_kron <- matrix(0, s * p, s * p)
#   G_m <- bdiag(matrix(1, ncol = 5, nrow = 5),
#                matrix(1, ncol = 7, nrow = 7),
#                matrix(1, ncol = 8, nrow = 8))
#   G_p <- bdiag(matrix(1, ncol = 5, nrow = 5), diag(3),
#                matrix(1, ncol = 11, nrow = 11), diag(9),
#                matrix(1, ncol = 21, nrow = 21),
#                matrix(1, ncol = 10, nrow = 10), diag(9),
#                matrix(1, ncol = 21, nrow = 21), diag(19),
#                matrix(1, ncol = 11, nrow = 11), diag(181))
#   G_kron <- kronecker(G_m, G_p)
# 
#   combn11 <- combn(rep((1:5 - 1) * p, each = length(1:5)) +
#                      rep(1:5, times = length(1:5)), 2)
#   combn12 <- combn(rep((1:5 - 1) * p, each = length(30:60)) +
#                      rep(30:60, times = length(1:5)), 2)
#   combn13 <- combn(rep((1:5 - 1) * p, each = length(110:120)) +
#                      rep(110:120, times = length(1:5)), 2)
#   combn21 <- combn(rep((6:12 - 1) * p, each = length(10:20)) +
#                      rep(10:20, times = length(6:12)), 2)
#   combn22 <- combn(rep((6:12 - 1) * p, each = length(51:60)) +
#                      rep(51:60, times = length(6:12)), 2)
#   combn23 <- combn(rep((6:12 - 1) * p, each = length(70:90)) +
#                      rep(70:90, times = length(6:12)), 2)
#   combn24 <- combn(rep((6:12 - 1) * p, each = length(110:120)) +
#                      rep(110:120, times = length(6:12)), 2)
#   combn31 <- combn(rep((13:15 - 1) * p, each = length(30:50)) +
#                      rep(30:50, times = length(13:15)), 2)
#   combn32 <- combn(rep((13:15 - 1) * p, each = length(70:90)) +
#                      rep(70:90, times = length(13:15)), 2)
#   combn33 <- combn(rep((13:15 - 1) * p, each = length(110:120)) +
#                      rep(110:120, times = length(13:15)), 2)
#   combn4 <- combn(rep((16:18 - 1) * p, each = length(121)) +
#                     rep(121, times = length(16:18)), 2)
#   combn5 <- matrix(rep((19 - 1) * p, each = length(122)) +
#                      rep(122, times = length(19)), nrow = 1, ncol = 2)
#   combn6 <- matrix(rep((20 - 1) * p, each = length(123)) +
#                      rep(123, times = length(20)), nrow = 1, ncol = 2)
# 
#   combnAll <- rbind(t(combn11), t(combn12), t(combn13),
#                     t(combn21), t(combn22), t(combn23), t(combn24),
#                     t(combn31), t(combn32), t(combn33),
#                     t(combn4), combn5, combn6)
# 
#   set.seed(seed + 7284)
#   sd_x <- 1
#   x <- matrix(rnorm(n * p, 0, sd_x), n, p)
# 
#   if (t0 == 0 & intercept) x <- cbind(rep(1, n), x)
#   if (!intercept) {
#     gamma <- gamma[-1, ]
#     b <- b[-1, ]
#   }
#   xb <- matrix(NA, n, s)
#   if (mv) {
#     for (i in 1:s) {
#       if (sum(gamma[, i]) >= 1) {
#         if (sum(gamma[, i]) == 1) {
#           xb[, i] <- x[, gamma[, i]] * b[gamma[, i], i]
#         } else {
#           xb[, i] <- x[, gamma[, i]] %*% b[gamma[, i], i]
#         }
#       } else {
#         xb[, i] <- sapply(1:s, function(i) rep(1, n) * b[1, i])
#       }
#     }
#   } else {
#     if (sum(gamma) >= 1) {
#       xb <- x[, gamma] %*% b[gamma, ]
#     } else {
#       xb <- sapply(1:s, function(i) rep(1, n) * b[1, i])
#     }
#   }
# 
#   corr_param <- 0.9
#   M <- matrix(corr_param, s, s)
#   diag(M) <- rep(1, s)
# 
#   ## wanna make it decomposable
#   Prime <- list(c(1:(s * .4), (s * .8):s),
#                 c((s * .4):(s * .6)),
#                 c((s * .65):(s * .75)),
#                 c((s * .8):s))
#   G <- matrix(0, s, s)
#   for (i in 1:length(Prime)) {
#     G[Prime[[i]], Prime[[i]]] <- 1
#   }
# 
#   # check
#   dimnames(G) <- list(1:s, 1:s)
#   length(gRbase::mcsMAT(G - diag(s))) > 0
# 
#   var <- solve(BDgraph::rgwish(n = 1, adj = G, b = 3, D = M))
# 
#   # change seeds to add randomness on error
#   set.seed(seed + 8493)
#   sd_err <- 0.5
#   if (is.infinite(t.df)) {
#     err <- matrix(rnorm(n * s, 0, sd_err), n, s) %*% chol(as.matrix(var))
#   } else {
#     err <- matrix(rt(n * s, t.df), n, s) %*% chol(as.matrix(var))
#   }
# 
#   if (t0 == 0) {
#     b.re <- NA
#     z <- NA
#     y <- xb + err
#     if (random.intercept != 0) {
#       y <- y + matrix(rnorm(n * s, 0, sqrt(random.intercept)), n, s)
#     }
# 
#     z <- sample(1:4, n, replace = T, prob = rep(1 / 4, 4))
# 
#     return(list(y = y, x = x, b = b, gamma = gamma, z = model.matrix(~ factor(z) + 0)[, ],
#                 b.re = b.re, Gy = G, mrfG = combnAll))
#   } else {
#     # add random effects
#     z <- t(rmultinom(n, size = 1, prob = c(.1, .2, .3, .4)))
#     z <- sample(1:t0, n, replace = T, prob = rep(1 / t0, t0))
#     set.seed(1683)
#     b.re <- rnorm(t0, 0, 2)
#     y <- matrix(b.re[z], nrow = n, ncol = s) + xb + err
# 
#     return(list(
#       y = y, x = x, b = b, gamma = gamma, z = model.matrix(~ factor(z) + 0)[, ],
#       b.re = b.re, Gy = G, mrfG = combnAll
#     ))
#   }
# }

## -----------------------------------------------------------------------------
# library("Matrix")
# n <- 250
# s <- 20
# p <- 300
# sim1 <- sim.ssur(n, s, p, seed = 1)

## -----------------------------------------------------------------------------
# t0 <- 4
# sim2 <- sim.ssur(n, s, p, t0, seed = 1) # learning data
# sim2.val <- sim.ssur(n, s, p, t0, seed=101) # validation data

## -----------------------------------------------------------------------------
# hyperpar <- list(mrf_d = -2, mrf_e = 1.6, a_w0 = 100, b_w0 = 500, a_w = 15, b_w = 60)
# set.seed(1038)
# fit2 <- BayesSUR(
#   data = cbind(sim2$y, sim2$z, sim2$x),
#   Y = 1:s,
#   X_0 = s + 1:t0,
#   X = s + t0 + 1:p,
#   outFilePath = "sim2_mrf_re",
#   hyperpar = hyperpar,
#   gammaInit = "0",
#   betaPrior = "reGroup",
#   nIter = 300, burnin = 100,
#   covariancePrior = "HIW",
#   standardize = F,
#   standardize.response = F,
#   gammaPrior = "MRF",
#   mrfG = sim2$mrfG,
#   output_CPO = T
# )

## -----------------------------------------------------------------------------
# summary(fit2)

## -----------------------------------------------------------------------------
# # show estimates
# plot(fit2, estimator=c("beta","gamma","Gy"), type="heatmap", name.predictors = "auto")

## -----------------------------------------------------------------------------
# # compute accuracy, sensitivity, specificity of variable selection
# gamma <- getEstimator(fit2)
# (accuracy <- sum(data.matrix(gamma > 0.5) == sim2$gamma) / prod(dim(gamma)))

## -----------------------------------------------------------------------------
# (sensitivity <- sum((data.matrix(gamma > 0.5) == 1) & (sim2$gamma == 1)) / sum(sim2$gamma == 1))

## -----------------------------------------------------------------------------
# (specificity <- sum((data.matrix(gamma > 0.5) == 0) & (sim2$gamma == 0)) / sum(sim2$gamma == 0))

## -----------------------------------------------------------------------------
# # compute RMSE and RMSPE for prediction performance
# beta <- getEstimator(fit2, estimator = "beta", Pmax = .5, beta.type = "conditional")
# (RMSE <- sqrt(sum((sim2$y - cbind(sim2$z, sim2$x) %*% beta)^2) / prod(dim(sim2$y))))

## -----------------------------------------------------------------------------
# (RMSPE <- sqrt(sum((sim2.val$y - cbind(sim2.val$z, sim2.val$x) %*% beta)^2) / prod(dim(sim2.val$y))))

## -----------------------------------------------------------------------------
# # compute bias of beta estimates
# b <- sim2$b
# b[sim2$gamma == 0] <- 0
# (beta.l2 <- sqrt(sum((beta[-c(1:4), ] - b)^2) / prod(dim(b))))

## -----------------------------------------------------------------------------
# g.re <- getEstimator(fit2, estimator = "Gy")
# (g.accuracy <- sum((g.re > 0.5) == sim2$Gy) / prod(dim(g.re)))

## -----------------------------------------------------------------------------
# (g.sensitivity <- sum(((g.re > 0.5) == sim2$Gy)[sim2$Gy == 1]) / sum(sim2$Gy == 1))

## -----------------------------------------------------------------------------
# (g.specificity <- sum(((g.re > 0.5) == sim2$Gy)[sim2$Gy == 0]) / sum(sim2$Gy == 0))