Implement baseline risk meta-regression

NAMESPACE

@@ -59,6 +59,7 @@ export(Surv)
 export(add_integration)
 export(as.stanfit)
 export(as.tibble.nodesplit_summary)
+export(calculate_baseline_risk)
 export(cauchy)
 export(combine_network)
 export(dbern)
@@ -121,6 +122,7 @@ importFrom(grDevices,nclass.Sturges)
 importFrom(graphics,pairs)
 importFrom(graphics,plot)
 importFrom(igraph,as.igraph)
+importFrom(rlang,"%||%")
 importFrom(rlang,":=")
 importFrom(rlang,.data)
 importFrom(rlang,abort)

R/data.R

@@ -400,3 +400,8 @@
 #'    \item{male}{proportion of male participants}
 #'  }
 "ndmm_agd_covs"
+
+#' Certolizumab
+#'
+#' Data analysed in \insertCite{TSD3}{multinma}.
+"certolizumab"

R/multinma-package.R

@@ -53,7 +53,7 @@
 #' @import methods
 #' @import Rcpp
 #' @importFrom dplyr %>%
-#' @importFrom rlang abort warn inform enquo .data :=
+#' @importFrom rlang abort warn inform enquo .data := %||%
 #' @importFrom rstan sampling
 #' @importFrom Rdpack reprompt
 #' @importFrom graphics plot pairs

R/nma.R

@@ -12,6 +12,9 @@
 #'   effect-modifying terms for a regression model. Any references to treatment
 #'   should use the `.trt` special variable, for example specifying effect
 #'   modifier interactions as `variable:.trt` (see details).
+#' @param baseline_risk Set to `TRUE` to fit a meta-regression on the baseline risk.
+#'   The baseline risk for centering purposed will be calculated from the network data.
+#'   The baseline risk can also be set manually by passing a numeric scalar.
 #' @param class_interactions Character string specifying whether effect modifier
 #'   interactions are specified as `"common"`, `"exchangeable"`, or
 #'   `"independent"`.
@@ -37,6 +40,8 @@
 #'   precision \eqn{1/\tau^2} (`"prec"`).
 #' @param prior_reg Specification of prior distribution for the regression
 #'   coefficients (if `regression` formula specified)
+#' @param prior_br Specification of prior distribution for the baseline risk
+#'   meta-regression coefficient (if `baseline_risk` specified).
 #' @param prior_aux Specification of prior distribution for the auxiliary
 #'   parameter, if applicable (see details). For `likelihood = "gengamma"` this
 #'   should be a list of prior distributions with elements `sigma` and `k`.
@@ -262,6 +267,7 @@ nma <- function(network,
                 consistency = c("consistency", "ume", "nodesplit"),
                 trt_effects = c("fixed", "random"),
                 regression = NULL,
+                baseline_risk = NULL,
                 class_interactions = c("common", "exchangeable", "independent"),
                 likelihood = NULL,
                 link = NULL,
@@ -272,6 +278,7 @@ nma <- function(network,
                 prior_het = .default(half_normal(scale = 5)),
                 prior_het_type = c("sd", "var", "prec"),
                 prior_reg = .default(normal(scale = 10)),
+                prior_br = .default(normal(scale = 10)),
                 prior_aux = .default(),
                 prior_aux_reg = .default(),
                 aux_by = NULL,
@@ -297,6 +304,7 @@ nma <- function(network,
 
   # Check model arguments
   consistency <- rlang::arg_match(consistency)
+  # TODO: unreachable, arg_match already does this
   if (length(consistency) > 1) abort("`consistency` must be a single string.")
   trt_effects <- rlang::arg_match(trt_effects)
   if (length(trt_effects) > 1) abort("`trt_effects` must be a single string.")
@@ -388,6 +396,7 @@ nma <- function(network,
                          prior_het = prior_het,
                          prior_het_type = prior_het_type,
                          prior_reg = prior_reg,
+                         prior_br = prior_br,
                          prior_aux = prior_aux,
                          prior_aux_reg = prior_aux_reg,
                          aux_by = aux_by,
@@ -515,6 +524,8 @@ nma <- function(network,
     prior_defaults$prior_het <- get_prior_call(prior_het)
   if (!is.null(regression) && !is_only_offset(regression) && .is_default(prior_reg))
     prior_defaults$prior_reg <- get_prior_call(prior_reg)
+  if (!is.null(baseline_risk) && .is_default(prior_br))
+    prior_defaults$prior_br <- get_prior_call(prior_br)
   if (has_aux && .is_default(prior_aux)) {
     if (likelihood == "normal" && has_ipd(network)) {
       prior_aux <- .default(half_normal(scale = 5))
@@ -891,6 +902,22 @@ nma <- function(network,
     .which_RE <- NULL
   }
 
+  if (!is.null(baseline_risk) && !isFALSE(baseline_risk)) {
+    if (!(likelihood %in% c("bernoulli", "binomial", "ordered", "normal"))) {
+      abort(glue::glue("Baseline risk meta-regression not yet supported with likelihood '{likelihood}'."))
+    }
+    if (is.numeric(baseline_risk)) {
+      baseline_risk <- link_fun(baseline_risk, link)
+    } else {
+      if (!isTRUE(baseline_risk)) {
+        abort("`baseline_risk` should be a single logical or numeric value.")
+      }
+      baseline_risk <- calculate_baseline_risk(network, link)
+    }
+  } else {
+    baseline_risk <- FALSE
+  }
+
   # Set up spline basis for mspline and pexp models
   if (likelihood %in% c("mspline", "pexp") && (has_ipd(network) || has_agd_arm(network))) {
     require_pkg("splines2")
@@ -1054,6 +1081,7 @@ nma <- function(network,
     agd_arm_offset = offset_agd_arm,
     agd_contrast_offset = offset_agd_contrast,
     trt_effects = trt_effects,
+    baseline_risk = baseline_risk,
     RE_cor = .RE_cor,
     which_RE = .which_RE,
     likelihood = likelihood,
@@ -1065,6 +1093,7 @@ nma <- function(network,
     prior_het = prior_het,
     prior_het_type = prior_het_type,
     prior_reg = prior_reg,
+    prior_br = prior_br,
     prior_aux = prior_aux,
     prior_aux_reg = prior_aux_reg,
     aux_id = aux_id,
@@ -1219,6 +1248,7 @@ nma <- function(network,
               consistency = consistency,
               regression = regression,
               aux_regression = aux_regression,
+              baseline_risk = baseline_risk,
               class_interactions = if (!is.null(regression) && !is.null(network$classes)) class_interactions else NULL,
               xbar = xbar,
               likelihood = likelihood,
@@ -1229,6 +1259,7 @@ nma <- function(network,
                             prior_het = if (trt_effects == "random") prior_het else NULL,
                             prior_het_type = if (trt_effects == "random") prior_het_type else NULL,
                             prior_reg = if (!is.null(regression) && !is_only_offset(regression)) prior_reg else NULL,
+                            prior_br = if (!is.null(baseline_risk)) prior_br else NULL,
                             prior_aux = if (has_aux) prior_aux else NULL,
                             prior_aux_reg = if (has_aux_regression) prior_aux_reg else NULL))
 
@@ -1270,6 +1301,7 @@ nma.fit <- function(ipd_x, ipd_y,
                     n_int,
                     ipd_offset = NULL, agd_arm_offset = NULL, agd_contrast_offset = NULL,
                     trt_effects = c("fixed", "random"),
+                    baseline_risk = NULL,
                     RE_cor = NULL,
                     which_RE = NULL,
                     likelihood = NULL,
@@ -1281,6 +1313,7 @@ nma.fit <- function(ipd_x, ipd_y,
                     prior_het,
                     prior_het_type = c("sd", "var", "prec"),
                     prior_reg,
+                    prior_br,
                     prior_aux,
                     prior_aux_reg,
                     aux_id = integer(),
@@ -1382,6 +1415,7 @@ nma.fit <- function(ipd_x, ipd_y,
 
   # Check model arguments
   trt_effects <- rlang::arg_match(trt_effects)
+  # TODO: unreachable
   if (length(trt_effects) > 1) abort("`trt_effects` must be a single string.")
 
   likelihood <- check_likelihood(likelihood)
@@ -1397,6 +1431,7 @@ nma.fit <- function(ipd_x, ipd_y,
   check_prior(prior_trt)
   if (trt_effects == "random") check_prior(prior_het)
   check_prior(prior_reg)
+  check_prior(prior_br)
   if (has_aux) {
     if (likelihood == "gengamma") check_prior(prior_aux, c("sigma", "k"))
     else check_prior(prior_aux)
@@ -1604,7 +1639,10 @@ nma.fit <- function(ipd_x, ipd_y,
     R_inv = if (QR) X_all_R_inv else matrix(0, 0, 0),
     # Offsets
     has_offset = has_offsets,
-    offsets = if (has_offsets) as.array(c(ipd_offset, agd_arm_offset, agd_contrast_offset)) else numeric()
+    offsets = if (has_offsets) as.array(c(ipd_offset, agd_arm_offset, agd_contrast_offset)) else numeric(),
+    baseline_risk = 0,
+    baseline_risk_flag = 0,
+    baseline_risk_which = numeric()
     )
 
   # Add priors
@@ -1615,6 +1653,8 @@ nma.fit <- function(ipd_x, ipd_y,
                       valid = c("Normal", "Cauchy", "Student t", "flat (implicit)")),
     !!! prior_standat(prior_reg, "prior_reg",
                       valid = c("Normal", "Cauchy", "Student t", "flat (implicit)")),
+    !!! prior_standat(prior_br, "prior_br",
+                      valid = c("Normal", "Cauchy", "Student t", "flat (implicit)")),
     !!! prior_standat(prior_het, "prior_het",
                       valid = c("Normal", "half-Normal", "log-Normal",
                                 "Cauchy",  "half-Cauchy",
@@ -1659,6 +1699,22 @@ nma.fit <- function(ipd_x, ipd_y,
   # Set chain_id to make CHAIN_ID available in data block
   stanargs$chain_id <- 1L
 
+  # Baseline risk meta-regression
+  if (!isFALSE(baseline_risk)) {
+    pars <- c(pars, "beta_baseline_risk")
+
+    studies <- X_all[, seq_len(standat$ns_agd_arm)]
+    is_reference <- unname(agd_arm_trt == 1L)
+    includes_reference <- as.logical(studies %*% (t(studies) %*% as.matrix(is_reference)))
+
+    baseline_risk_which <- (!is_reference) * includes_reference
+
+    standat <-
+      standat %>%
+      purrr::list_modify(baseline_risk = baseline_risk, baseline_risk_flag = 1L,
+                         baseline_risk_which = baseline_risk_which)
+  }
+
   # Call Stan model for given likelihood
 
   # -- Normal likelihood

R/nma_data.R

@@ -1810,3 +1810,54 @@ nfactor <- function(x, ..., numeric = TRUE, resort = FALSE) {
     return(factor(x, levels = stringr::str_sort(unique(x), numeric = numeric), ...))
   }
 }
+
+#' Calculate baseline risk
+#'
+#' @param network An `nma_data` object.
+#' @param link The name of the link function to use.
+#' @param treatment The reference treatment. If `NULL` (the default), use the `network`
+#'   reference treatment.
+#'
+#' @return The baseline risk as a single numeric value.
+#' @export
+calculate_baseline_risk <- function(network, link, treatment = NULL) {
+  if (!inherits(network, "nma_data")) abort("Not nma_data object.")
+
+  if (!has_agd_arm(network) || has_agd_contrast(network) || has_ipd(network)) {
+    abort("Should only have agd_arm.")
+  }
+
+  treatments <- levels(network$treatments)
+  treatment <- treatment %||% treatments[1L]
+  treatment <- rlang::arg_match(treatment, treatments)
+
+  agd_arm <- network$agd_arm[network$agd_arm$.trt == treatment, ]
+
+  if (network$outcome$agd_arm == "continuous") {
+    out <- agd_arm$.y
+
+  } else if (network$outcome$agd_arm %in% c("count", "ordered")) {
+    r <- agd_arm$.r
+    if (inherits(r, "matrix")) {
+      n <- rowSums(r, na.rm = TRUE)
+      r <- rowSums(r[, -1L], na.rm = TRUE)
+    } else {
+      n <- agd_arm$.n
+    }
+
+    out <- r / n
+
+    # Assume odds of 0.01 if r is 0
+    out[r == 0] <- 0.01 / (1 + 0.01)
+
+    stopifnot(all(out >= 0 & out <= 1))
+
+  } else {
+    abort(paste(
+      "Calculation of baseline risk not yet implemented for",
+      network$outcome$agd_arm, "outcomes."
+    ))
+  }
+
+  mean(link_fun(out, link))
+}

inst/stan/binomial_1par.stan

@@ -34,6 +34,10 @@ transformed parameters {
       X_agd_arm * beta_tilde + offset_agd_arm :
       X_agd_arm * beta_tilde;
 
+    if (baseline_risk_flag) {
+      eta_agd_arm_noRE += beta_baseline_risk[1] * baseline_risk_which .* (X_agd_arm[, 1:ns_agd_arm] * beta_tilde[1:ns_agd_arm] - baseline_risk);
+    }
+
     if (nint_max > 1) { // -- If integration points are used --
       if (RE) {
         if (link == 1) { // logit link

inst/stan/binomial_2par.stan

@@ -37,6 +37,10 @@ transformed parameters {
       X_agd_arm * beta_tilde + offset_agd_arm :
       X_agd_arm * beta_tilde;
 
+    if (baseline_risk_flag) {
+      eta_agd_arm_noRE += beta_baseline_risk[1] * baseline_risk_which .* (X_agd_arm[, 1:ns_agd_arm] * beta_tilde[1:ns_agd_arm] - baseline_risk);
+    }
+
     if (nint_max > 1) { // -- If integration points are used --
       if (RE) {
 

inst/stan/include/data_common.stan

@@ -76,3 +76,12 @@ int<lower=0,upper=3> prior_reg_dist;
 real prior_reg_location;
 real<lower=0> prior_reg_scale;
 real<lower=0> prior_reg_df;
+
+int<lower=0,upper=3> prior_br_dist;
+real prior_br_location;
+real<lower=0> prior_br_scale;
+real<lower=0> prior_br_df;
+
+real baseline_risk;
+int<lower=0, upper=1> baseline_risk_flag;
+vector[baseline_risk_flag ? ni_agd_arm : 0] baseline_risk_which;

inst/stan/include/model_common.stan

@@ -7,6 +7,8 @@ prior_select_lp(mu, prior_intercept_dist, prior_intercept_location, prior_interc
 prior_select_lp(d, prior_trt_dist, prior_trt_location, prior_trt_scale, prior_trt_df);
 // Regression parameters
 prior_select_lp(beta, prior_reg_dist, prior_reg_location, prior_reg_scale, prior_reg_df);
+// Baseline risk meta-regression
+prior_select_lp(beta_baseline_risk, prior_br_dist, prior_br_location, prior_br_scale, prior_br_df);
 
 // Node-splitting - just use prior for d
 prior_select_lp(omega, prior_trt_dist, prior_trt_location, prior_trt_scale, prior_trt_df);

inst/stan/include/parameters_common.stan

@@ -5,3 +5,5 @@ vector[nX] beta_tilde;
 // (Uncorrelated) random effects and heterogeneity SD
 vector[n_delta] u_delta;
 vector<lower = 0>[RE ? 1 : 0] tau;
+
+vector[baseline_risk_flag ? 1 : 0] beta_baseline_risk;

inst/stan/normal.stan

@@ -39,6 +39,10 @@ transformed parameters {
       X_agd_arm * beta_tilde + offset_agd_arm :
       X_agd_arm * beta_tilde;
 
+    if (baseline_risk_flag) {
+      eta_agd_arm_noRE += beta_baseline_risk[1] * baseline_risk_which .* (X_agd_arm[, 1:ns_agd_arm] * beta_tilde[1:ns_agd_arm] - baseline_risk);
+    }
+
     if (nint_max > 1) { // -- If integration points are used --
       if (RE) {
 

inst/stan/ordered_multinomial.stan

@@ -87,6 +87,10 @@ transformed parameters {
           X_agd_arm * beta_tilde + offset_agd_arm :
           X_agd_arm * beta_tilde;
 
+    if (baseline_risk_flag) {
+      eta_agd_arm_noRE += beta_baseline_risk[1] * baseline_risk_which .* (X_agd_arm[, 1:ns_agd_arm] * beta_tilde[1:ns_agd_arm] - baseline_risk);
+    }
+
     if (nint_max > 1) { // -- If integration points are used --
 
       if (RE) {