multinma (version 0.6.1)

plot_integration_error: Plot numerical integration error

Description

For ML-NMR models, plot the estimated numerical integration error over the entire posterior distribution, as the number of integration points increases. See methods_paper,Phillippo_thesismultinma for details.

Usage

plot_integration_error(
  x,
  ...,
  stat = "violin",
  orientation = c("vertical", "horizontal", "x", "y"),
  show_expected_rate = TRUE
)

Value

A ggplot object.

Arguments

x

An object of type stan_mlnmr

...

Additional arguments passed to the ggdist plot stat.

stat

Character string specifying the ggdist plot stat used to summarise the integration error over the posterior. Default is "violin", which is equivalent to "eye" with some cosmetic tweaks.

orientation

Whether the ggdist geom is drawn horizontally ("horizontal") or vertically ("vertical"), default "vertical"

show_expected_rate

Logical, show typical convergence rate \(1/N\)? Default TRUE.

Note for survival models

This function is not supported for survival/time-to-event models. These do not save cumulative integration points for efficiency reasons (both time and memory).

Details

The total number of integration points is set by the n_int argument to add_integration(), and the intervals at which integration error is estimated are set by the int_thin argument to nma(). The typical convergence rate of Quasi-Monte Carlo integration (as used here) is \(1/N\), which by default is displayed on the plot output.

The integration error at each thinning interval \(N_\mathrm{thin}\) is estimated for each point in the posterior distribution by subtracting the final estimate (using all n_int points) from the estimate using only the first \(N_\mathrm{thin}\) points.

Examples

## Plaque psoriasis ML-NMR
# Set up plaque psoriasis network combining IPD and AgD
library(dplyr)
pso_ipd <- filter(plaque_psoriasis_ipd,
                  studyc %in% c("UNCOVER-1", "UNCOVER-2", "UNCOVER-3"))

pso_agd <- filter(plaque_psoriasis_agd,
                  studyc == "FIXTURE")

head(pso_ipd)
head(pso_agd)

pso_ipd <- pso_ipd %>%
  mutate(# Variable transformations
    bsa = bsa / 100,
    prevsys = as.numeric(prevsys),
    psa = as.numeric(psa),
    weight = weight / 10,
    durnpso = durnpso / 10,
    # Treatment classes
    trtclass = case_when(trtn == 1 ~ "Placebo",
                         trtn %in% c(2, 3, 5, 6) ~ "IL blocker",
                         trtn == 4 ~ "TNFa blocker"),
    # Check complete cases for covariates of interest
    complete = complete.cases(durnpso, prevsys, bsa, weight, psa)
  )

pso_agd <- pso_agd %>%
  mutate(
    # Variable transformations
    bsa_mean = bsa_mean / 100,
    bsa_sd = bsa_sd / 100,
    prevsys = prevsys / 100,
    psa = psa / 100,
    weight_mean = weight_mean / 10,
    weight_sd = weight_sd / 10,
    durnpso_mean = durnpso_mean / 10,
    durnpso_sd = durnpso_sd / 10,
    # Treatment classes
    trtclass = case_when(trtn == 1 ~ "Placebo",
                         trtn %in% c(2, 3, 5, 6) ~ "IL blocker",
                         trtn == 4 ~ "TNFa blocker")
  )

# Exclude small number of individuals with missing covariates
pso_ipd <- filter(pso_ipd, complete)

pso_net <- combine_network(
  set_ipd(pso_ipd,
          study = studyc,
          trt = trtc,
          r = pasi75,
          trt_class = trtclass),
  set_agd_arm(pso_agd,
              study = studyc,
              trt = trtc,
              r = pasi75_r,
              n = pasi75_n,
              trt_class = trtclass)
)

# Print network details
pso_net

# Add integration points to the network
pso_net <- add_integration(pso_net,
  durnpso = distr(qgamma, mean = durnpso_mean, sd = durnpso_sd),
  prevsys = distr(qbern, prob = prevsys),
  bsa = distr(qlogitnorm, mean = bsa_mean, sd = bsa_sd),
  weight = distr(qgamma, mean = weight_mean, sd = weight_sd),
  psa = distr(qbern, prob = psa),
  n_int = 64)

# \donttest{
# Fit the ML-NMR model
pso_fit <- nma(pso_net, refresh = if (interactive()) 200 else 0,
               trt_effects = "fixed",
               link = "probit",
               likelihood = "bernoulli2",
               regression = ~(durnpso + prevsys + bsa + weight + psa)*.trt,
               class_interactions = "common",
               prior_intercept = normal(scale = 10),
               prior_trt = normal(scale = 10),
               prior_reg = normal(scale = 10),
               init_r = 0.1,
               QR = TRUE,
               # Set the thinning factor for saving the cumulative results
               # (This also sets int_check = FALSE)
               int_thin = 8)
pso_fit

# Plot numerical integration error
plot_integration_error(pso_fit)
# }