Add npost to control number of post processing samples

haziqj · May 28, 2024 · 9304f32 · 9304f32
1 parent 5aa7a00
commit 9304f32
Show file tree

Hide file tree

Showing 2 changed files with 14 additions and 19 deletions.
diff --git a/R/30-inlavaan.R b/R/30-inlavaan.R
@@ -5,6 +5,7 @@
 #' @param dp TBC
 #' @param save.lvs TBC
 #' @param bcontrol TBC
+#' @param npost Number of post-processing samples to do, if any.
 #'
 #' @return TBC
 #' @export
@@ -15,6 +16,7 @@ inlavaan <- function(
     target = "INLA",
     dp = NULL,
     save.lvs = FALSE,
+    npost = 250,
     bcontrol = list(num.threads = 6)) {
 
   # To play nice with blavaan code
@@ -880,10 +882,12 @@ inlavaan <- function(
       #   pxpartable = jagtrans$pxpartable,
       #   res = res
       # ))
+      if(verbose) cat("Post processing... ")
       parests <- coeffun_inla(
         lavpartable = lavpartable,
         pxpartable = jagtrans$pxpartable,
-        res = res
+        res = res,
+        npost = npost
       )
       stansumm <- parests$stansumm
     # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
@@ -1141,7 +1145,7 @@ inlavaan <- function(
                             jagextra            = mcmcextra,
                             stansumm            = stansumm,
                             domll               = domll)
-    if(verbose) cat(" done.\n")
+    if(verbose) cat("done.\n")
   }
   timing$TEST <- (proc.time()[3] - start.time)
   start.time <- proc.time()[3]

diff --git a/R/40-lav_export_INLA.R b/R/40-lav_export_INLA.R
@@ -1330,7 +1330,7 @@ coeffun_inla <- function(
     pxpartable,
     res,
     fun = "mean",
-    nsamp = 100) {
+    npost = 2000) {
 
   # Get INLA estimates ---------------------------------------------------------
   idx_nu    <- pxpartable$free[pxpartable$mat == "nu" & pxpartable$free > 0]
@@ -1375,21 +1375,17 @@ coeffun_inla <- function(
 
   # Psi (may require sampling)
   if (length(idx_lvrho) > 0) {
-    # samps <-
-    #   purrr::map(res$internal.marginals.hyperpar[c(idx_psi, idx_lvrho)], \(m) {
-    #     INLA::inla.rmarginal(nsamp, m)
-    #   }) |>
-    #   as.data.frame()
     marginals <- res$internal.marginals.hyperpar[c(idx_psi, idx_lvrho)]
     samps <- lapply(marginals, \(m) {
-      INLA::inla.rmarginal(nsamp, m)
+      INLA::inla.rmarginal(npost, m)
     })
     samps <- as.data.frame(samps)
 
     samps <- apply(samps, 1, simplify = FALSE, \(x) {
       psi <- exp(x[seq_along(idx_psi)])
-      lvrho <- exp(x[length(idx_psi) + seq_along(idx_lvrho)])
-      lvrho <- lvrho / (1 + lvrho)
+      u <- exp(-x[length(idx_psi) + seq_along(idx_lvrho)])
+      u <- 1 / (1 + u)
+      lvrho <- 2 * u - 1
 
       SD <- Diagonal(x = sqrt(psi))
       Rho_df <- pxpartable[pxpartable$mat == "lvrho", ]
@@ -1453,21 +1449,16 @@ coeffun_inla <- function(
 
   # Theta may require sampling
   if (length(idx_rho) > 0) {
-    # samps <-
-    #   purrr::map(res$internal.marginals.hyperpar[c(idx_theta, idx_rho)], \(m) {
-    #     INLA::inla.rmarginal(nsamp, m)
-    #   }) |>
-    #   as.data.frame()
     marginals <- res$internal.marginals.hyperpar[c(idx_theta, idx_rho)]
     samps <- lapply(marginals, \(m) {
-      INLA::inla.rmarginal(nsamp, m)
+      INLA::inla.rmarginal(npost, m)
     })
     samps <- as.data.frame(samps)
 
     samps <- apply(samps, 1, simplify = FALSE, \(x) {
       theta_e <- exp(x[seq_along(idx_theta)])
-      u <- exp(x[length(idx_theta) + seq_along(idx_rho)])
-      u <- u / (1 + u)
+      u <- exp(-x[length(idx_theta) + seq_along(idx_rho)])
+      u <- 1 / (1 + u)
       rho <- 2 * u - 1
 
       SD <- Diagonal(x = sqrt(theta_e))