more updates to the release

Novartis · Nov 21, 2024 · dee0638 · dee0638
1 parent 5b83632
commit dee0638
Show file tree

Hide file tree

Showing 10 changed files with 40 additions and 25 deletions.
diff --git a/Makefile b/Makefile
@@ -79,7 +79,7 @@ src/package-binary: R/stanmodels.R
 	echo "importFrom(rstan, sampling)" >> NAMESPACE
 	echo "useDynLib($(RPKG), .registration = TRUE)" >> NAMESPACE
 	install -d src
-	"${R_HOME}/bin/Rscript" -e 'pkgbuild::compile_dll()'
+	"${R_HOME}/bin/Rscript" -e 'pkgbuild::compile_dll(debug=FALSE)'
 	touch src/package-binary
 
 man/package-doc: $(R_PKG_SRCS) $(BIN_OBJS)

diff --git a/NEWS.md b/NEWS.md
@@ -2,12 +2,17 @@
 
 ## Enhancements
 
-* Upgraded testhat edition to version 3
 * Added for `mixstanvar` automatic generation of distribution
   functions allowing truncated priors in `brms` with mixtures
+* Slight speed increase for Stan model by more efficient construction
+  of likelihood. Also reducing object size of `gMAP` objects by
+  avoiding to store redundant variables.
+* Upgraded testhat edition to version 3.
 
 ## Bugfixes
 
+* Fix issue #18 of ESS ELIR aborting whenever one mixture component
+  has zero weight.
 * Fixed a rare issue when estimating ESS ELIR for beta mixtures. The
   calculation aborted due to boundary issues which are now avoided.
 * Avoid over and underflow of beta-binomial distribution.

diff --git a/R/plot_gMAP.R b/R/plot_gMAP.R
@@ -11,7 +11,7 @@
 #'
 #' @template plot-help
 #' 
-#' @return The function returns a list of \code{\link{ggplot}}
+#' @return The function returns a list of \code{\link[ggplot2:ggplot]{ggplot}}
 #' objects.
 #'
 #' @method plot gMAP

diff --git a/inst/examples/ess.R b/inst/examples/ess.R
@@ -24,14 +24,12 @@ round(sum(ab_matched))
 ess(bmix, method="morita")
 
 # Predictive consistency of elir
-# \donttest{
 n_forward <- 1E2
 bmixPred <- preddist(bmix, n=n_forward)
 pred_samp <- rmix(bmixPred, 1E3)
 pred_ess <- sapply(pred_samp, function(r) ess(postmix(bmix, r=r, n=n_forward), "elir") )
 ess(bmix, "elir")
 mean(pred_ess) - n_forward
-# }
 
 # Normal mixture example
 nmix <- mixnorm(rob=c(0.5, 0, 2), inf=c(0.5, 3, 4), sigma=10)

diff --git a/inst/sbc/calibration.md5 b/inst/sbc/calibration.md5
@@ -1,3 +1,3 @@
-Created:  2024-11-20 13:57:46 UTC
-git hash: c17af53badaeb321d3de398b84d998dff5cf685f
-MD5:      95c9b0f2f899285fb53b68663ec1b6e3
+Created:  2024-11-21 10:30:15 UTC
+git hash: 5f5afb101f98533e4229f9dc267b3540b172c573
+MD5:      c8f566f591c59887653bac53ab03bb45
diff --git a/inst/sbc/calibration.rds b/inst/sbc/calibration.rds
diff --git a/inst/sbc/make_reference_rankhist.R b/inst/sbc/make_reference_rankhist.R
@@ -23,7 +23,7 @@ library(data.table)
 library(knitr)
 sbc_tools <- new.env()
 source(here("inst", "sbc", "sbc_tools.R"), local=sbc_tools)
-set.seed(453453)
+set.seed(45348346)
 
 sbc_tools$rbest_lib_dir <- rbest_lib_dir
 sbc_tools$rbest_source_dir <- rbest_source_dir
@@ -80,7 +80,7 @@ num_simulations <- nrow(scenarios)
 cat("Total number of jobs to dispatch:", num_simulations, "\n")
 
 RNGkind("L'Ecuyer-CMRG")
-set.seed(56969)
+set.seed(269698974)
 rng_seeds <- sbc_tools$setup_lecuyer_seeds(.Random.seed, num_simulations)
 
 sim_result <- Q_rows(scenarios, sbc_tools$run_sbc_case, const=list(base_scenarios=base_scenarios, seeds=rng_seeds), export=as.list(sbc_tools), n_jobs=n_jobs, pkgs=pkg)

diff --git a/inst/sbc/sbc_tools.R b/inst/sbc/sbc_tools.R
@@ -87,12 +87,16 @@ fit_rbest <- function(fake, draw, draw_theta, family, prior_mean_mu, prior_sd_mu
                     gaussian=cbind(y, y_se) ~ 1 | group)
 
     options(RBesT.MC.warmup=2000, RBesT.MC.iter=4000, RBesT.MC.thin=1, RBesT.MC.init=0.1,
-            RBesT.MC.control=list(adapt_delta=0.999, stepsize=0.01, max_treedepth=10,
+            RBesT.MC.control=list(##adapt_delta=0.999, ## 2024-11-21 lowered to 0.95 for the sake of performance
+                                  adapt_delta=0.95,
+                                  stepsize=0.01, max_treedepth=10,
                                   adapt_init_buffer=100, adapt_term_buffer=300))
 
     if(Ng > 5) {
         ## for the dense case we can be a bit less aggressive with the sampling tuning parameters
-        options(RBesT.MC.control=list(adapt_delta=0.95, stepsize=0.01, max_treedepth=10,
+        options(RBesT.MC.control=list(##adapt_delta=0.95,
+                                      adapt_delta=0.90,
+                                      stepsize=0.01, max_treedepth=10,
                                       adapt_init_buffer=100, adapt_term_buffer=300))
     }
     fit <- gMAP(model, data=fake, family=family,

diff --git a/inst/stan/gMAP.stan b/inst/stan/gMAP.stan
@@ -134,10 +134,8 @@ parameters {
 }
 transformed parameters {
   vector[H] theta;
-  vector[n_groups] eta;
   vector[mX] beta;
   vector[n_tau_strata] tau;
-  vector[n_groups] tau_group;
 
   beta = beta_raw_guess[1] + beta_raw_guess[2] .* beta_raw;
 
@@ -147,14 +145,24 @@ transformed parameters {
   else
     tau = exp(tau_raw_guess[1] + tau_raw_guess[2] * tau_raw);
 
-  tau_group = tau[tau_strata_gindex];
-
-  if (ncp)  // NCP
-    eta = xi_eta .* tau_group;
-  else // CP places overall intercept into random effect
-    eta = beta_raw_guess[1,1] + beta_raw_guess[2,1] * xi_eta;
-
-  theta = X_param * beta + eta[group_index];
+  // expand random effect to groups in loop for performance reasons
+  if(ncp) {
+    if(n_tau_strata == 1) {
+      // most common case of just one stratum which simplifies things
+      // and in ncp mode
+      for(h in 1:H) {
+        theta[h] = X_param[h] * beta + xi_eta[ group_index[h] ] * tau[1];
+      }
+    } else {
+      for(h in 1:H) {
+        theta[h] = X_param[h] * beta + xi_eta[ group_index[h] ] * tau[ tau_strata_gindex[ group_index[h] ] ];
+      }
+    }
+  } else {
+    for(h in 1:H) {
+      theta[h] = X_param[h] * beta + beta_raw_guess[1,1] + beta_raw_guess[2,1] * xi_eta[ group_index[h] ];
+    }
+  }
 }
 model {
   if (ncp) {
@@ -163,8 +171,8 @@ model {
     if(re_dist == 1) xi_eta ~ student_t(re_dist_t_df, 0, 1);
   } else {
     // random effect distribution
-    if(re_dist == 0) xi_eta ~ normal( (beta[1] - beta_raw_guess[1,1])/beta_raw_guess[2,1], tau_group / beta_raw_guess[2,1]);
-    if(re_dist == 1) xi_eta ~ student_t(re_dist_t_df, (beta[1] - beta_raw_guess[1,1])/beta_raw_guess[2,1], tau_group / beta_raw_guess[2,1]);
+    if(re_dist == 0) xi_eta ~ normal( (beta[1] - beta_raw_guess[1,1])/beta_raw_guess[2,1], tau[tau_strata_gindex] / beta_raw_guess[2,1]);
+    if(re_dist == 1) xi_eta ~ student_t(re_dist_t_df, (beta[1] - beta_raw_guess[1,1])/beta_raw_guess[2,1], tau[tau_strata_gindex] / beta_raw_guess[2,1]);
   }
 
   // assign priors to coefficients

diff --git a/tools/make-ds.R b/tools/make-ds.R
@@ -43,7 +43,7 @@ make_internal_ds <- function() {
     calibration_meta["MD5"]  <- vals["MD5"]
 
     pkg_create_date  <- Sys.time()
-    pkg_sha <- "1f6edf8"
+    pkg_sha <- "a0c573a"
 
     if (gsub("\\$", "", pkg_sha) == "Format:%h") {
         pkg_sha <- system("git rev-parse --short HEAD", intern=TRUE)