changed from hardcoded A,B,C to something more generic. Its still a bit hacky but the binomial data examples seems to work ok. Need to check the other data types n ext.

closes #35

Author: n8thangreen

NAMESPACE

@@ -31,7 +31,6 @@ export(strategy_gcomp_stan)
 export(strategy_maic)
 export(strategy_mim)
 export(strategy_stc)
-import(plyr)
 import(stats)
 import(stringr)
 import(tidyr)

R/calc_IPD_stats.R

@@ -74,8 +74,14 @@ calc_IPD_stats.mim <- function(strategy,
   lci.Delta <- coef_est + qt(0.025, df = nu) * sqrt(var_est)
   uci.Delta <- coef_est + qt(0.975, df = nu) * sqrt(var_est)
 
-  list(mean = coef_est,
-       var = var_est)
+  list(
+    contrasts = list(
+      mean = coef_est,
+      var = var_est),
+    absolute = list(
+      mean = NA,  #p_est,  ##TODO:
+      var = NA)   #p_var)
+  )
 } 
 
 #' Factory function for creating calc_IPD_stats methods

R/calc_stc.R

@@ -5,7 +5,7 @@
 #' 
 calc_stc <- function(strategy, ipd, ...) {
 
-  treat_nm <- strategy$trt_var
+  trt_var <- strategy$trt_var
 
   # centre covariates
   centre_vars <- get_eff_mod_names(strategy$formula)
@@ -19,8 +19,8 @@ calc_stc <- function(strategy, ipd, ...) {
   # extract model coefficients
   coef_fit <- coef(fit)
 
-  # safer than treat_nm in case of factor level append
-  treat_coef_name <- grep(paste0("^", treat_nm), names(coef_fit), value = TRUE)
+  # safer than trt_var in case of factor level append
+  treat_coef_name <- grep(paste0("^", trt_var, "[^:]*$"), names(coef_fit), value = TRUE)
 
   # probability for control and treatment group
   # estimating treatment effect at means because of centring

R/gcomp_ml.R

@@ -23,10 +23,12 @@
 gcomp_ml.boot <- function(data, indices,
                           R, formula = NULL,
                           family, trt_var,
+                          ref_trt = "C", comp_trt = "A",
                           rho = NA,
                           N = 1000, ald) {
   dat <- data[indices, ]
-  gcomp_ml_means(formula, family, dat, ald, trt_var, rho, N) 
+  gcomp_ml_means(formula, family, dat, ald, trt_var, rho, N,
+                 ref_trt, comp_trt) 
 }
 
 
@@ -62,7 +64,9 @@ gcomp_ml_means <- function(formula,
                            ipd, ald,
                            trt_var,
                            rho = NA,
-                           N = 1000) {
+                           N = 1000,
+                           ref_trt = "C",
+                           comp_trt = "A") {
 
   x_star <- simulate_ALD_pseudo_pop(formula = formula,
                                     ipd = ipd, ald = ald,
@@ -75,19 +79,19 @@ gcomp_ml_means <- function(formula,
              data = ipd)
 
   # counterfactual datasets
-  data.trtA <- data.trtC <- x_star
+  data.comp <- data.ref <- x_star
 
   # intervene on treatment while keeping set covariates fixed
-  data.trtA[[trt_var]] <- 0  # all receive A
-  data.trtC[[trt_var]] <- 1  # all receive C
+  data.comp[[trt_var]] <- comp_trt  # all receive A
+  data.ref[[trt_var]] <- ref_trt    # all receive C
 
   # predict counterfactual event probs, conditional on treatment/covariates
-  hat.mu.A <- predict(fit, type="response", newdata=data.trtA)
-  hat.mu.C <- predict(fit, type="response", newdata=data.trtC)
+  hat.mu.comp <- predict(fit, type="response", newdata=data.comp)
+  hat.mu.ref <- predict(fit, type="response", newdata=data.ref)
 
   # (marginal) mean probability prediction under A and C
-  c(`0` = mean(hat.mu.A),
-    `1` = mean(hat.mu.C))
+  c(`0` = mean(hat.mu.ref),
+    `1` = mean(hat.mu.comp))
 }
 
 

R/gcomp_stan.R

@@ -37,7 +37,9 @@
 #' @export
 #'
 calc_gcomp_stan <- function(strategy,
-                            ipd, ald, ...) {
+                            ipd, ald, 
+                            ref_trt = "C",
+                            comp_trt = "A", ...) {
 
   formula <- strategy$formula
   family <- strategy$family
@@ -56,21 +58,21 @@ calc_gcomp_stan <- function(strategy,
                        ...)
 
   # counterfactual datasets
-  data.trtA <- data.trtC <- x_star
+  data.comp <- data.ref <- x_star
 
   # intervene on treatment while keeping set covariates fixed
-  data.trtA[[trt_var]] <- 1  # everyone receives treatment A
-  data.trtC[[trt_var]] <- 0  # receive treatment C
+  data.comp[[trt_var]] <- comp_trt  # all receive treatment A
+  data.ref[[trt_var]] <- ref_trt    # all receive treatment C
 
   ##TODO: is this going to work for all of the different data types?
   # draw responses from posterior predictive distribution
-  y.star.A <- rstanarm::posterior_predict(outcome.model, newdata = data.trtA)
-  y.star.C <- rstanarm::posterior_predict(outcome.model, newdata = data.trtC)
+  y.star.comp <- rstanarm::posterior_predict(outcome.model, newdata = data.comp)
+  y.star.ref <- rstanarm::posterior_predict(outcome.model, newdata = data.ref)
 
   # posterior means for each treatment group
   list(
-    mean_A = rowMeans(y.star.A),
-    mean_C = rowMeans(y.star.C))
+    mean_A = rowMeans(y.star.comp),
+    mean_C = rowMeans(y.star.ref))
 }
 
 
@@ -116,6 +118,8 @@ calc_gcomp_ml <- function(strategy,
          formula = strategy$formula,
          family = strategy$family,
          trt_var = strategy$trt_var,
+         ref_trt = "C",
+         comp_trt = "A",
          rho = strategy$rho,
          N = strategy$N,
          data = ipd,

R/get_study_treatments.R

@@ -9,7 +9,9 @@
 #' @keywords internal
 #' 
 calc_mim <- function(strategy,
-                     ipd, ald, ...) {
+                     ipd, ald, 
+                     ref_trt = "C",
+                     comp_trt = "A", ...) {
 
   formula <- strategy$formula
   family <- strategy$family
@@ -29,11 +31,15 @@ calc_mim <- function(strategy,
     algorithm = "sampling", ...)
 
   # create augmented target dataset
-  target.t1 <- target.t0 <- x_star
-  target.t1$trt_var <- 1
-  target.t0$trt_var <- 0
+  target.comp <- target.ref <- x_star
+  target.comp[[trt_var]] <- comp_trt
+  target.ref[[trt_var]] <- ref_trt
 
-  aug.target <- rbind(target.t0, target.t1)
+  aug.target <- rbind(target.ref, target.comp)
+  
+  # set reference treatment as base level
+  aug.target[[trt_var]] <- factor(aug.target[[trt_var]],
+                                  levels = c(ref_trt, comp_trt))
 
   # complete syntheses by drawing binary outcomes
   # from their posterior predictive distribution
@@ -52,18 +58,21 @@ calc_mim <- function(strategy,
     glm(as.formula(paste("y ~", trt_var)), data = data_m, family = family)
   })
 
-  
   # treatment effect point estimates in each synthesis
   coef_fit <- do.call(rbind, lapply(reg2.fits, function(fit) coef(fit)))
 
+  # safer than trt_var in case of factor level append
+  coef_names <- names(coef(reg2.fits[[1]]))
+  treat_coef_name <- grep(pattern = paste0("^", trt_var, "[^:]*$"), coef_names, value = TRUE)
+  
   ##TODO: how to transform this to the prob scale?
   # point estimates for the variance in each synthesis
   hats.v <- unlist(lapply(reg2.fits,
                           function(fit)
-                            vcov(fit)[trt_var, trt_var]))
+                            vcov(fit)[treat_coef_name, treat_coef_name]))
 
-  mean_C <- family$linkinv(coef_fit[, 1])                  # probability for control
-  mean_A <- family$linkinv(coef_fit[, 1] + coef_fit[, 2])  # probability for treatment
+  mean_C <- family$linkinv(coef_fit[, 1])                  # probability for reference
+  mean_A <- family$linkinv(coef_fit[, 1] + coef_fit[, treat_coef_name])  # probability for comparator
 
   tibble::lst(mean_A, mean_C,
               hats.v, M)

R/mim.R

@@ -6,10 +6,18 @@
 #' Methods taken from
 #' \insertCite{RemiroAzocar2022}{outstandR}.
 #' 
-#' @param ipd_trial Individual-level patient data. Suppose between studies _A_ and _C_.
-#' @param ald_trial Aggregate-level data. Suppose between studies _B_ and _C_. 
+#' @param ipd_trial Individual-level patient data. For example, suppose between studies _A_ and _C_.
+#'   In a long format and must contain a treatment column and outcome column consistent with the formula object.
+#'   The labels in the treatment are used internally so there must be a common treatment with the aggregate-level data trial.
+#' @param ald_trial Aggregate-level data. For example, suppose between studies _B_ and _C_. The column names take the form
+#'  - `mean.X*`: mean patient measurement
+#'  - `sd.X*`: standard deviation of patient measurement
+#'  - `y.*.sum`: total number of events
+#'  - `y.*.bar`: proportion of events
+#'  - `N.*`: total number of individuals
 #' @param strategy Computation strategy function. These can be
 #'    `strategy_maic()`, `strategy_stc()`, `strategy_gcomp_ml()` and  `strategy_gcomp_stan()`
+#' @param ref_trt Reference / common / anchoring treatment name; default "C"
 #' @param CI Confidence interval; between 0,1
 #' @param scale Relative treatment effect scale. If `NULL`, the scale is automatically determined from the model.
 #' @param ... Additional arguments
@@ -27,44 +35,63 @@
 #' data(AC_IPD)  # AC patient-level data
 #' data(BC_ALD)  # BC aggregate-level data
 #' 
+#' # linear formula
 #' lin_form <- as.formula("y ~ X3 + X4 + trt*X1 + trt*X2")
 #' 
 #' # matching-adjusted indirect comparison
-#' outstandR_maic <- outstandR(AC_IPD, BC_ALD, strategy = strategy_maic(formula = lin_form))
+#' outstandR_maic <- outstandR(AC_IPD, BC_ALD,
+#'                             strategy = strategy_maic(formula = lin_form))
 #' 
 #' # simulated treatment comparison
-#' outstandR_stc <- outstandR(AC_IPD, BC_ALD, strategy = strategy_stc(lin_form))
+#' outstandR_stc <- outstandR(AC_IPD, BC_ALD,
+#'                            strategy = strategy_stc(lin_form))
 #' 
 #' # G-computation with maximum likelihood
-#' # outstandR_gcomp_ml <- outstandR(AC_IPD, BC_ALD, strategy = strategy_gcomp_ml(lin_form))
+#' # outstandR_gcomp_ml <- outstandR(AC_IPD, BC_ALD,
+#'                                   strategy = strategy_gcomp_ml(lin_form))
 #' 
 #' # G-computation with Bayesian inference
-#' outstandR_gcomp_stan <- outstandR(AC_IPD, BC_ALD, strategy = strategy_gcomp_stan(lin_form))
+#' outstandR_gcomp_stan <- outstandR(AC_IPD, BC_ALD,
+#'                                   strategy = strategy_gcomp_stan(lin_form))
 #' 
 #' # Multiple imputation marginalization
-#' outstandR_gcomp_stan <- outstandR(AC_IPD, BC_ALD, strategy = strategy_mim(lin_form))
+#' outstandR_mim <- outstandR(AC_IPD, BC_ALD,
+#'                            strategy = strategy_mim(lin_form))
 #' 
-outstandR <- function(ipd_trial, ald_trial, ref_trt = "C",
-                      strategy, CI = 0.95, scale = NULL, ...) {
+outstandR <- function(ipd_trial, ald_trial, strategy,
+                      ref_trt = "C",
+                      CI = 0.95, scale = NULL, ...) {
 
   validate_outstandr(ipd_trial, ald_trial, strategy, CI, scale)
 
   ipd <- prep_ipd(strategy$formula, ipd_trial)
   ald <- prep_ald(strategy$formula, ald_trial, trt_var = strategy$trt_var)
 
-  ipd_trts <- get_ipd_trts(ipd, ref_trt)
-  ald_trts <- get_ald_trts(ald, ref_trt)
+  # treatment names for each study
+  
+  ipd_comp <- get_ipd_comparator(ipd, ref_trt, strategy$trt_var)
+  ald_comp <- get_ald_comparator(ald, ref_trt)
+  
+  ipd_trts <- list(ipd_comp, ref_trt)
+  ald_trts <- list(ald_comp, ref_trt)
 
   if (is.null(scale)) scale <- get_treatment_effect(strategy$family$link)
 
-  ipd_stats <- calc_IPD_stats(strategy, ipd = ipd, ald = ald, scale, ...) 
-  ald_stats <- calc_ALD_stats(strategy, ald = ald, treatments = ald_trts, scale = scale) 
+  ipd_stats <- calc_IPD_stats(strategy, 
+                              ipd = ipd, ald = ald,
+                              scale, ...) 
+  
+  ald_stats <- calc_ALD_stats(strategy, ald = ald,
+                              treatments = ald_trts,
+                              scale = scale) 
 
   stats <- result_stats(ipd_stats, ald_stats, CI)
 
-  structure(stats,
-            CI = CI,
-            scale = scale,
-            model = strategy$family$family,
-            class = c("outstandR", class(stats)))
+  structure(
+    stats,
+    CI = CI,
+    ref_trt = ref_trt,
+    scale = scale,
+    model = strategy$family$family,
+    class = c("outstandR", class(stats)))
 }

R/outstandR.R

@@ -39,7 +39,6 @@ prep_ald <- function(form, data, trt_var = "trt") {
 
 #' Convert from long to wide format
 #'
-#' @import plyr
 #' @import tidyr
 #' @import stringr
 #' 
@@ -75,3 +74,25 @@ reshape_ald_to_wide <- function(df) {
 
     bind_cols(variable_df, y_df)
 }
+
+# Get study comparator treatment names
+
+#
+get_ald_comparator <- function(ald, ref_trt = "C") {
+  
+  pattern <- paste0("^y\\.(?!", ref_trt, ")")
+  
+  # filter for names starting with "y." but not with "y.C"
+  y_non_ref <- grep(pattern, colnames(ald),
+                    value = TRUE, perl = TRUE)
+  
+  # extract treatment of first match
+  sub("^y\\.([A-Z]).*", "\\1", y_non_ref[1])
+}
+
+#
+get_ipd_comparator <- function(ipd, ref_trt = "C", trt_var = "trt") {
+  all_trt <- levels(ipd[[trt_var]])
+  all_trt[all_trt != ref_trt]
+}
+

R/prep_data.R

@@ -14,10 +14,12 @@ print.outstandR <- function(x, ...) {
     stop("Please install the 'tibble' and 'pillar' packages for colored tibble output.")
   }
 
+  ref_trt <- attr(x, "ref_trt")
+  
   cat(pillar::style_bold("Object of class 'outstandR'"), "\n")
   cat("Model:", pillar::style_subtle(attr(x, "model")), "\n")
   cat("Scale:", pillar::style_subtle(attr(x, "scale")), "\n")
-  cat("Common treatment:", pillar::style_subtle("C"), "\n")
+  cat("Common treatment:", pillar::style_subtle(ref_trt), "\n")
   # cat(pillar::style_subtle("Common treatment:"), attr(x, "reference"), "\n")  ##TODO:
   cat("Individual patient data study:", pillar::style_subtle("AC"), "\n")
   cat("Aggregate level data study:", pillar::style_subtle("BC"), "\n")
@@ -46,8 +48,8 @@ print.outstandR <- function(x, ...) {
     Treatments = names(absolute$means),
     Estimate = unlist(absolute$means),
     Std.Error = unlist(absolute$variances),
-    lower.0.95 = sapply(absolute$CI, \(x) x[1]),
-    upper.0.95 = sapply(absolute$CI, \(x) x[2])
+    lower.0.95 = c(NA,NA),  #sapply(absolute$CI, \(x) x[1]),  ##TODO:
+    upper.0.95 = c(NA,NA)   #sapply(absolute$CI, \(x) x[2])
   )
 
   cat("Contrasts:\n\n")