merge cran accepted version with main

Merge branch 'cran-submission'

# Conflicts:
#	man/outstandR-package.Rd

Author: n8thangreen

.Rbuildignore

@@ -28,4 +28,5 @@ scripts/create_package_data\.R
 ^R/conversion-functions\.R$
 ^cran-comments\.md$
 ^.*\.qmd$
-^.*\.Rmd$
+^.*\.Rmd$
+^check_results$

.gitignore

@@ -5,3 +5,7 @@
 inst/doc
 /Meta/
 /doc/
+inst/jss-section.html
+inst/jss-section.qmd
+inst/jss-section_files/
+check_results/

DESCRIPTION

@@ -11,9 +11,9 @@ Authors@R: c(
 Description:
     For the problem of indirect treatment comparison with limited subject-level data, this 
     package provides tools for model-based standardisation with several different computation approaches.
-    See Remiro‐Azócar A, Heath A, Baio G (2022. Parametric G‐computation for compatible
-    indirect treatment comparisons with limited individual patient data, 
-    Res. Synth. Methods, 1–31. ISSN 1759-2879, <doi:10.1002/jrsm.1565>.)
+        See Remiro‐Azócar A, Heath A, Baio G (2022) ``Parametric G‐computation for compatible
+        indirect treatment comparisons with limited individual patient data'',
+        Res. Synth. Methods, 1–31. ISSN 1759-2879, <doi:10.1002/jrsm.1565>.
 License: GPL (>= 3)
 Encoding: UTF-8
 Language: en-GB

NAMESPACE

@@ -1,5 +1,7 @@
 # Generated by roxygen2: do not edit by hand
 
+S3method(add_seed,default)
+S3method(add_seed,gcomp_bayes)
 S3method(calc_IPD_stats,default)
 S3method(calc_IPD_stats,gcomp_bayes)
 S3method(calc_IPD_stats,gcomp_ml)

R/calc_IPD_stats.R

@@ -95,16 +95,16 @@ IPD_stat_factory <- function(ipd_fun) {
       var_est <- var(hat.delta.AC, na.rm = TRUE)
     } else if (var_method == "rubin") {
 
-      if (is.null(out$model$hats.v) || is.null(out$model$M)) {
-        stop("Rubin's rules require 'hats.v' and 'M' in the model output.", call. = FALSE)
+      if (is.null(out$model$hats.v) || is.null(out$model$n_imp)) {
+        stop("Rubin's rules require 'hats.v' and 'n_imp' in the model output.", call. = FALSE)
       }
 
-      M <- out$model$M
+      n_imp <- out$model$n_imp
       bar.v <- mean(out$model$hats.v)
       b <- var(hat.delta.AC, na.rm = TRUE)
 
-      var_est <- var_by_pooling(M, bar.v, b)
-      out$model$nu <- wald_type_interval(M, bar.v, b)
+      var_est <- var_by_pooling(n_imp, bar.v, b)
+      out$model$nu <- wald_type_interval(n_imp, bar.v, b)
 
     } else {
       # This should be caught by get_var_method, but serves as a failsafe

R/cor2cov.R

@@ -1,12 +1,12 @@
 
 #' Compute covariance matrix
 #'
-#' @param R correlation matrix
+#' @param cormat correlation matrix
 #' @param S vector of standard deviations covariance matrix
 #' @return Required input for \code{mvrnorm}.
 #' @keywords internal
 #'
-cor2cov <- function(R, S) {
-  # R * S %*% t(S)  # alternative
-  sweep(sweep(R, 1, S, "*"), 2, S, "*")
+cor2cov <- function(cormat, S) {
+  # cormat * S %*% t(S)  # alternative
+  sweep(sweep(cormat, 1, S, "*"), 2, S, "*")
 }

R/gcomp_bayes.R

@@ -59,11 +59,19 @@
 calc_gcomp_bayes <- function(strategy,
                              analysis_params, ...) {
 
+  # extract seed, or default to a random integer if missing/NULL
+  bayes_seed <- analysis_params$seed
+  
+  if (is.null(bayes_seed)) {
+    bayes_seed <- sample.int(.Machine$integer.max, 1)
+  }
+  
   default_stan_args <- list(
     algorithm = "sampling",
     chains = 2,
     iter = 2000,
-    refresh = 0  # quiet
+    refresh = 0,  # quiet
+    seed = bayes_seed
   )
 
   # merge with user-provided dots
@@ -144,7 +152,7 @@ calc_gcomp_bayes <- function(strategy,
 #'   family = binomial(),
 #'   rho = NA,
 #'   N = 1000L,
-#'   R = 100L,
+#'   n_boot = 100L,
 #'   marginal_distns = NA,
 #'   marginal_params = NA,
 #'   trt_var = "trt")
@@ -188,7 +196,7 @@ calc_gcomp_ml <- function(strategy,
   args_boot <- c(
     common_args, list(
       data = analysis_params$ipd, 
-      R = strategy$R))
+      R = strategy$n_boot))
 
   gcomp_boot <- do.call(boot::boot, c(statistic = gcomp_ml.boot, args_boot))
 

R/gcomp_ml.R

@@ -23,7 +23,8 @@
 #' @keywords internal
 #' 
 gcomp_ml.boot <- function(data, indices,
-                          R, formula = NULL,
+                          R, 
+                          formula = NULL,
                           family, trt_var,
                           ref_trt = NA,
                           comp_trt = NA,

R/maic.R

@@ -189,7 +189,7 @@ maic.boot <- function(ipd, indices = 1:nrow(ipd),
 calc_maic <- function(strategy,
                       analysis_params) {
   args_list <- 
-    list(R = strategy$R,
+    list(R = strategy$n_boot,
          formula = strategy$formula,
          family = strategy$family,
          trt_var = strategy$trt_var,

R/mim.R

@@ -7,13 +7,13 @@
 #' @param ... Additional argument to pass to Stan model
 #' 
 #' @return A list containing:
-#' * `means`: A list containing vectors of posterior means (one per synthesis `M`):
+#' * `means`: A list containing vectors of posterior means (one per synthesis `n_imp`):
 #'     * `A`: Comparator means.
 #'     * `C`: Reference means.
 #' * `model`: A list containing:
 #'     * `fit`: The first-stage [rstanarm::stan_glm()] object.
 #'     * `hats.v`: Vector of variance point estimates for each synthesis.
-#'     * `M`: Number of posterior prediction draws (syntheses).
+#'     * `n_imp`: Number of posterior prediction draws (syntheses).
 #'     * `rho`, `N`, `stan_args`: Strategy and model parameters.
 #'
 #' @importFrom rstanarm posterior_predict stan_glm
@@ -72,10 +72,10 @@ calc_mim <- function(strategy,
 
   # ANALYSIS STAGE ---
 
-  M <- nrow(y_star)
+  n_imp <- nrow(y_star)
 
   # fit second-stage regression to each synthesis using maximum-likelihood estimation
-  reg2.fits <- lapply(1:M, function(m) {
+  reg2.fits <- lapply(1:n_imp, function(m) {
     data_m <- aug.target
     data_m$y <- y_star[m, ]
     glm(as.formula(paste("y ~", trt_var)), data = data_m, family = family)
@@ -103,7 +103,7 @@ calc_mim <- function(strategy,
     model = list(
       fit = outcome_model,
       hats.v = hats.v,
-      M = M,
+      n_imp = n_imp,
       rho = rho,
       N = N,
       stan_args = stan_args)
@@ -114,27 +114,27 @@ calc_mim <- function(strategy,
 #' 
 #' Constructed using t-distribution with nu degrees of freedom.
 #'
-#' @param M Number of syntheses used in analysis stage (high for low Monte Carlo error)
+#' @param n_imp Number of syntheses used in analysis stage (high for low Monte Carlo error)
 #' @param bar.v "within" variance (average of variance point estimates)
 #' @param b "between" variance (sample variance of point estimates)
 #' @return Numeric value of Wald-type interval estimates.
 #' @keywords internal
 #' 
-wald_type_interval <- function(M, bar.v, b) {
-  (M - 1) * (1 + bar.v / ((1 + 1 / M) * b)) ^ 2
+wald_type_interval <- function(n_imp, bar.v, b) {
+  (n_imp - 1) * (1 + bar.v / ((1 + 1 / n_imp) * b)) ^ 2
 }
 
 #' Variance estimate by pooling
 #' 
 #' Use combining rules to estimate.
 #' 
-#' @param M Number of syntheses used in analysis stage (high for low Monte Carlo error)
+#' @param n_imp Number of syntheses used in analysis stage (high for low Monte Carlo error)
 #' @param bar.v "within" variance (average of variance point estimates)
 #' @param b "between" variance (sample variance of point estimates)
 #' @return Numeric value of variance estimate using pooling.
 #' @keywords internal
 #' 
-var_by_pooling <- function(M, bar.v, b) {
-  (1 + (1 / M)) * b - bar.v
+var_by_pooling <- function(n_imp, bar.v, b) {
+  (1 + (1 / n_imp)) * b - bar.v
 }