plot_integration_error {multinma} | R Documentation |
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 (Phillippo et al. 2020; Phillippo 2019) for details.
Usage
plot_integration_error(
x,
...,
stat = "violin",
orientation = c("vertical", "horizontal", "x", "y"),
show_expected_rate = TRUE
)
Arguments
x |
An object of type |
... |
Additional arguments passed to the |
stat |
Character string specifying the |
orientation |
Whether the |
show_expected_rate |
Logical, show typical convergence rate |
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
, which by default is displayed on the plot output.
The integration error at each thinning interval 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 points.
Value
A ggplot
object.
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).
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)
# Fit the ML-NMR model
pso_fit <- nma(pso_net,
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)