drafted adding in mim code. need to finish

R/IPD_stats.R

@@ -120,3 +120,18 @@ IPD_stats.gcomp_stan <- function(strategy,
        var = var(hat.delta.AC))
 } 
 
+
+#' @rdname IPD_stats
+#' @section Multiple imputation marginalisation:
+#'
+#' @export
+#'
+IPD_stats.mim <- function(strategy,
+                          ipd, ald) {
+
+  hat.delta.AC <- mim(ipd, ald)
+
+  list(mean = mean(hat.delta.AC),
+       var = var(hat.delta.AC))
+} 
+

R/mim.R

@@ -0,0 +1,71 @@
+
+#' Multiple imputation marginalization (MIM)
+#' 
+#' @param ipd.index,ipd.target Index RCT and target covariate datasets
+#' @param M Number of syntheses used in analysis stage (high for low Monte Carlo error)
+#'
+mim <- function(ipd.index,
+                ipd.target,
+                M = 1,
+                n.chains = 2,
+                warmup = 10,
+                iters = 1000) {
+
+  ## SYNTHESIS STAGE ##
+
+  # first-stage logistic regression model fitted to index RCT using MCMC (Stan)
+  outcome.model <- stan_glm(
+    y ~ (trt * x1) + (trt * x2),
+    data = ipd.index,
+    family = binomial,
+    algorithm = "sampling",
+    iter = iters,
+    warmup = warmup,
+    chains = n.chains,
+    # thin to use M independent samples in analysis stage
+    thin = (n.chains * (iters - warmup)) / M
+  )
+
+  # create augmented target dataset
+  target.t1 <- target.t0 <- ipd.target
+  target.t1$trt <- 1
+  target.t0$trt <- 0
+
+  aug.target <- rbind(target.t0, target.t1)
+
+  # complete syntheses by drawing binary outcomes from their posterior predictive distribution
+  y.star <- posterior_predict(outcome.model, newdata = aug.target)
+
+  ## ANALYSIS STAGE ##
+
+  # fit second-stage regression to each synthesis using maximum-likelihood estimation
+  reg2.fits <- lapply(1:M, function(m)
+    glm(y.star[m, ] ~ trt, data = aug.target, family = binomial))
+
+  # treatment effect point estimates in each synthesis
+  hats.delta <- unlist(lapply(reg2.fits, function(fit)
+    coef(fit)["trt"][[1]]))
+
+  # point estimates for the variance in each synthesis
+  hats.v <- unlist(lapply(reg2.fits, function(fit)
+    vcov(fit)["trt", "trt"]))
+
+  # quantities originally defined by Rubin (1987) for multiple imputation
+  bar.delta <- mean(hats.delta)  # average of treatment effect point estimates
+  bar.v <- mean(hats.v)          # "within" variance (average of variance point estimates)
+  b <- var(hats.delta)           # "between" variance (sample variance of point estimates)
+
+  # pooling: average of point estimates is marginal log odds ratio
+  hat.Delta <- bar.delta
+
+  # pooling: use combining rules to estimate the variance
+  hat.var.Delta <- (1 + (1 / M)) * b - bar.v
+
+  # Wald-type interval estimates constructed using t-distribution with nu degrees of freedom
+  nu <- (M - 1) * (1 + bar.v / ((1 + 1 / M) * b)) ^ 2
+
+  lci.Delta <- hat.Delta + qt(0.025, df = nu) * sqrt(hat.var.Delta)
+  uci.Delta <- hat.Delta + qt(0.975, df = nu) * sqrt(hat.var.Delta)
+
+  hat.Delta
+}

R/strategy_.R

@@ -153,6 +153,22 @@ strategy_gcomp_stan <- function(formula = NULL) {
   do.call(new_strategy, c(strategy = "gcomp_stan", args))
 }
 
+#' @rdname strategy
+#' 
+#' @section Multiple imputation marginalization (MIM):
+#' 
+#' @return `mim` class object
+#' @export
+# 
+strategy_mim <- function(formula = NULL) {
+  check_formula(formula)
+
+  default_args <- formals()
+  args <- c(formula = formula, as.list(match.call())[-c(1,2)])
+  args <- modifyList(default_args, args)
+  do.call(new_strategy, c(strategy = "mim", args))
+}
+
 #' @name strategy
 #' @title New strategy objects
 #'