Merge branch 'main' of https://github.com/StatisticsHealthEconomics/mimR

NAMESPACE

@@ -21,3 +21,8 @@ export(strategy_stc)
 export(trial_treatment_effect)
 export(trial_variance)
 importFrom(Rdpack,reprompt)
+importFrom(copula,mvdc)
+importFrom(copula,normalCopula)
+importFrom(copula,rMvdc)
+importFrom(rstanarm,posterior_predict)
+importFrom(rstanarm,stan_glm)

R/IPD_stats.R

@@ -34,9 +34,7 @@
 #' @export
 #'
 strategy_maic <- function(formula = as.formula("y ~ X3 + X4 + trt*X1 + trt*X2"),
-                          R = 1000,
-                          ald) {
-
+                          R = 1000) {
   if (class(formula) != "formula")
     stop("formula argument must be of formula class.")
 
@@ -247,7 +245,7 @@ new_strategy <- function(strategy, ...) {
 hat_Delta_stats <- function(AC.IPD, BC.ALD, strategy, CI = 0.95, ...) {
 
   if (CI <= 0 || CI >= 1) stop("CI argument must be between 0 and 1.")
-  
+
   ##TODO: as method instead?
   if (!inherits(strategy, "strategy"))
     stop("strategy argument must be of a class strategy.")
@@ -329,7 +327,7 @@ IPD_stats.maic <- function(strategy,
                           statistic = maic.boot,
                           R = strategy$R,
                           formula = strategy$formula,
-                          ald = strategy$dat_ALD)
+                          ald = ald)
 
   list(mean = mean(maic_boot$t),
        var = var(maic_boot$t))

R/gcomp_ml.R

@@ -30,6 +30,7 @@ gcomp_ml.boot <- function(data, indices,
 #'
 #' @return Difference of log-odds 
 #' @seealso [strategy_gcomp_ml()], [gcomp_ml.boot()]
+#' @importFrom copula normalCopula mvdc rMvdc
 #' @export
 #'
 gcomp_ml_log_odds_ratio <- function(formula, dat) {
@@ -56,9 +57,9 @@ gcomp_ml_log_odds_ratio <- function(formula, dat) {
   # AC IPD pairwise correlations
   t_rho <- t(rho)  # extract along rows
   cop <-
-    normalCopula(param = t_rho[lower.tri(t_rho)],
-                 dim = n_covariates,
-                 dispstr = "un")
+    copula::normalCopula(param = t_rho[lower.tri(t_rho)],
+                         dim = n_covariates,
+                         dispstr = "un")
 
   # BC covariate means & standard deviations
   mean_sd_margins <- list()
@@ -69,12 +70,12 @@ gcomp_ml_log_odds_ratio <- function(formula, dat) {
   }
 
   # sample covariates from approximate joint distribution using copula
-  mvd <- mvdc(copula = cop,
-              margins = rep("norm", n_covariates),  # Gaussian marginals
-              paramMargins = mean_sd_margins)
+  mvd <- copula::mvdc(copula = cop,
+                      margins = rep("norm", n_covariates),  # Gaussian marginals
+                      paramMargins = mean_sd_margins)
 
   # simulated BC pseudo-population
-  x_star <- as.data.frame(rMvdc(n = 1000, mvd))
+  x_star <- as.data.frame(copula::rMvdc(n = 1000, mvd))
   colnames(x_star) <- covariate_names
 
   # outcome logistic regression fitted to IPD using maximum likelihood

R/gcomp_stan.R

@@ -9,6 +9,8 @@
 #' @template args-ald
 #'
 #' @return A list of \eqn{y^*_A} and \eqn{y^*_C} posterior predictions
+#' @importFrom copula normalCopula mvdc rMvdc
+#' @importFrom rstanarm stan_glm posterior_predict
 #' @export
 #'
 gcomp_stan <- function(formula = as.formula("y ~ X3 + X4 + trt*X1 + trt*X2"),
@@ -27,11 +29,11 @@ gcomp_stan <- function(formula = as.formula("y ~ X3 + X4 + trt*X1 + trt*X2"),
 
   # covariate simulation for BC trial using copula package
   cop <-
-    normalCopula(param = c(rho[1,2],rho[1,3],rho[1,4],
-                                    rho[2,3],rho[2,4],
-                                             rho[3,4]),
-                 dim = n_covariates,
-                 dispstr = "un")  # AC IPD pairwise correlations
+    copula::normalCopula(param = c(rho[1,2],rho[1,3],rho[1,4],
+                                   rho[2,3],rho[2,4],
+                                   rho[3,4]),
+                         dim = n_covariates,
+                         dispstr = "un")  # AC IPD pairwise correlations
 
   # sample covariates from approximate joint distribution using copula
   mvd <- mvdc(copula = cop,
@@ -48,11 +50,11 @@ gcomp_stan <- function(formula = as.formula("y ~ X3 + X4 + trt*X1 + trt*X2"),
   colnames(x_star) <- cov_names
 
   # outcome logistic regression fitted to IPD using MCMC (Stan)
-  outcome.model <- stan_glm(formula,
-                            data = ipd,
-                            family = binomial,
-                            algorithm = "sampling",
-                            iter = 4000, warmup = 2000, chains = 2)
+  outcome.model <- rstanarm::stan_glm(formula,
+                                      data = ipd,
+                                      family = binomial,
+                                      algorithm = "sampling",
+                                      iter = 4000, warmup = 2000, chains = 2)
   # counterfactual datasets
   data.trtA <- data.trtC <- x_star
 
@@ -62,7 +64,7 @@ gcomp_stan <- function(formula = as.formula("y ~ X3 + X4 + trt*X1 + trt*X2"),
 
   # draw binary responses from posterior predictive distribution
   list(
-    y.star.A = posterior_predict(outcome.model, newdata=data.trtA),
-    y.star.C = posterior_predict(outcome.model, newdata=data.trtC))
+    y.star.A = rstanarm::posterior_predict(outcome.model, newdata=data.trtA),
+    y.star.C = rstanarm::posterior_predict(outcome.model, newdata=data.trtC))
 }
 

R/gen_data.R

@@ -1,20 +1,36 @@
 
-#' Generate simulated datasets for index and comparator trials
+#' Generate simulated datasets of IPD covariates and outcome for a trial
 #' 
-#' @param N Size
-#' @param b_trt b coefficient for treatment 
-#' @param b_X b coefficients for X 
-#' @param b_EM b coefficients effect modifiers
+#' @param N Total number of patients
+#' @param b_trt `b` coefficient for active treatment vs. common comparator 
+#' @param b_X `b` coefficients for each prognostic variable `X` 
+#' @param b_EM `b` coefficients effect modifiers
 #' @param b_0 Intercept coefficient
-#' @param meanX mean X 
-#' @param sdX standard deviation X 
+#' @param meanX Mean of each normally-distributed covariate `X` 
+#' @param sdX Standard deviation of each covariate `X` 
 #' @param event_rate Event rate 
-#' @param corX Correlation of X 
-#' @param allocation Allocation
+#' @param corX Covariate correlation coefficient of `X` 
+#' @param allocation Allocation to active treatment as proportion of total; 0 to 1
 #' @return Data frame of `X`, `trt` and `y`
 #' 
 #' @keywords internal
 #' 
+#' @examples
+#' 
+#' x <- gen_data(
+#'  N = 100,
+#'  b_trt = log(0.17),
+#'  b_X = -log(0.5),
+#'  b_EM = -log(0.67),
+#'  b_0 = -0.62,
+#'  meanX = 0.6,
+#'  sdX = 0.4,
+#'  event_rate = 0.35, 
+#'  corX = 0.2,
+#'  allocation = 2/3) 
+#' 
+#' head(x)
+#' 
 gen_data <- function(N, b_trt, b_X, b_EM, b_0,
                      meanX, sdX, event_rate, 
                      corX, allocation) {
@@ -26,13 +42,15 @@ gen_data <- function(N, b_trt, b_X, b_EM, b_0,
   sd.vec <- rep(sdX, 4)            # vector of standard deviations
 
   cov.mat <- cor2cov(rho, sd.vec)  # covariance matrix
+
   # simulate correlated continuous covariates using multivariate normal
   # patients under active treatment
   X_active <- 
     as.data.frame(
       MASS::mvrnorm(n = N_active,
                     mu = rep(meanX, 4),
                     Sigma = cov.mat))
+
   # patients under control treatment
   X_control <-
     as.data.frame(
@@ -41,7 +59,8 @@ gen_data <- function(N, b_trt, b_X, b_EM, b_0,
                     Sigma = cov.mat))  
   # all patients
   X <- rbind(X_active, X_control)
-  colnames(X) <- c("X1","X2","X3","X4") 
+  colnames(X) <- c("X1", "X2", "X3", "X4")
+
   # treatment assignment (1: active; 0: control)
   trt <- c(rep(1,N_active), rep(0,N_control))
 
@@ -52,7 +71,7 @@ gen_data <- function(N, b_trt, b_X, b_EM, b_0,
     b_X*X$X3 + b_X*X$X4 + b_trt*trt +
     b_EM*X$X1*trt + b_EM*X$X2*trt
 
-  yprob <- 1/(1 + exp(-LP))               # binary outcome probability
+  yprob <- 1/(1 + exp(-LP))             # binary outcome probability
   y <- rbinom(n=N, size=1, prob=yprob)  # binary outcome
 
   as.data.frame(cbind(X, trt, y))