In dmphillippo/multinma: Bayesian Network Meta-Analysis of Individual and Aggregate Data
set.seed(1042020)
library(multinma)
options(mc.cores = parallel::detectCores())
library(multinma)
nc <- switch(tolower(Sys.getenv("_R_CHECK_LIMIT_CORES_")),
"true" =, "warn" = 2,
parallel::detectCores())
options(mc.cores = nc)
This vignette describes the analysis of data on the mean off-time reduction in patients given dopamine agonists as adjunct therapy in Parkinson's disease, in a network of 7 trials of 4 active drugs plus placebo [@TSD2].
The data are available in this package as parkinsons:
head(parkinsons)
We consider analysing these data in three separate ways:
- Using arm-based data (means
y and corresponding standard errors se);
- Using contrast-based data (mean differences
diff and corresponding standard errors se_diff);
- A combination of the two, where some studies contribute arm-based data, and other contribute contrast-based data.
Note: In this case, with Normal likelihoods for both arms and contrasts, we will see that the three analyses give identical results.
In general, unless the arm-based likelihood is Normal, results from a model using a contrast-based likelihood will not exactly match those from a model using an arm-based likelihood, since the contrast-based Normal likelihood is only an approximation.
Similarity of results depends on the suitability of the Normal approximation, which may not always be appropriate - e.g. with a small number of events or small sample size for a binary outcome.
The use of an arm-based likelihood (sometimes called an "exact" likelihood) is therefore preferable where possible in general.
Analysis of arm-based data
We begin with an analysis of the arm-based data - means and standard errors.
Setting up the network
We have arm-level continuous data giving the mean off-time reduction (y) and standard error (se) in each arm.
We use the function set_agd_arm() to set up the network.
arm_net <- set_agd_arm(parkinsons,
study = studyn,
trt = trtn,
y = y,
se = se,
sample_size = n)
arm_net
We let treatment 4 be set by default as the network reference treatment, since this results in considerably improved sampling efficiency over choosing treatment 1 as the network reference.
The sample_size argument is optional, but enables the nodes to be weighted by sample size in the network plot.
Plot the network structure.
plot(arm_net, weight_edges = TRUE, weight_nodes = TRUE)
Meta-analysis models
We fit both fixed effect (FE) and random effects (RE) models.
Fixed effect meta-analysis
First, we fit a fixed effect model using the nma() function with trt_effects = "fixed".
We use $\mathrm{N}(0, 100^2)$ prior distributions for the treatment effects $d_k$ and study-specific intercepts $\mu_j$.
We can examine the range of parameter values implied by these prior distributions with the summary() method:
summary(normal(scale = 100))
The model is fitted using the nma() function.
arm_fit_FE <- nma(arm_net,
trt_effects = "fixed",
prior_intercept = normal(scale = 100),
prior_trt = normal(scale = 10))
arm_fit_FE <- nma(arm_net,
trt_effects = "fixed",
prior_intercept = normal(scale = 100),
prior_trt = normal(scale = 10),
iter = 10000)
Basic parameter summaries are given by the print() method:
arm_fit_FE
By default, summaries of the study-specific intercepts $\mu_j$ are hidden, but could be examined by changing the pars argument:
# Not run
print(arm_fit_FE, pars = c("d", "mu"))
The prior and posterior distributions can be compared visually using the plot_prior_posterior() function:
plot_prior_posterior(arm_fit_FE)
Random effects meta-analysis
We now fit a random effects model using the nma() function with trt_effects = "random".
Again, we use $\mathrm{N}(0, 100^2)$ prior distributions for the treatment effects $d_k$ and study-specific intercepts $\mu_j$, and we additionally use a $\textrm{half-N}(5^2)$ prior for the heterogeneity standard deviation $\tau$.
We can examine the range of parameter values implied by these prior distributions with the summary() method:
summary(normal(scale = 100))
summary(half_normal(scale = 5))
Fitting the RE model
arm_fit_RE <- nma(arm_net,
seed = 379394727,
trt_effects = "random",
prior_intercept = normal(scale = 100),
prior_trt = normal(scale = 100),
prior_het = half_normal(scale = 5),
adapt_delta = 0.99)
arm_fit_RE <- nowarn_on_ci(nma(arm_net,
seed = 379394727,
trt_effects = "random",
prior_intercept = normal(scale = 100),
prior_trt = normal(scale = 100),
prior_het = half_normal(scale = 5),
adapt_delta = 0.99,
iter = 10000))
We do see a small number of divergent transition errors, which cannot simply be removed by increasing the value of the adapt_delta argument (by default set to 0.95 for RE models).
To diagnose, we use the pairs() method to investigate where in the posterior distribution these divergences are happening (indicated by red crosses):
pairs(arm_fit_RE, pars = c("mu[4]", "d[3]", "delta[4: 3]", "tau"))
The divergent transitions occur in the upper tail of the heterogeneity standard deviation.
In this case, with only a small number of studies, there is not very much information to estimate the heterogeneity standard deviation and the prior distribution may be too heavy-tailed.
We could consider a more informative prior distribution for the heterogeneity variance to aid estimation.
Basic parameter summaries are given by the print() method:
arm_fit_RE
By default, summaries of the study-specific intercepts $\mu_j$ and study-specific relative effects $\delta_{jk}$ are hidden, but could be examined by changing the pars argument:
# Not run
print(arm_fit_RE, pars = c("d", "mu", "delta"))
The prior and posterior distributions can be compared visually using the plot_prior_posterior() function:
plot_prior_posterior(arm_fit_RE)
Model comparison
Model fit can be checked using the dic() function:
(arm_dic_FE <- dic(arm_fit_FE))
(arm_dic_RE <- dic(arm_fit_RE))
Both models fit the data well, having posterior mean residual deviance close to the number of data points.
The DIC is similar between models, so we choose the FE model based on parsimony.
We can also examine the residual deviance contributions with the corresponding plot() method.
plot(arm_dic_FE)
plot(arm_dic_RE)
Further results
For comparison with @TSD2, we can produce relative effects against placebo using the relative_effects() function with trt_ref = 1:
(arm_releff_FE <- relative_effects(arm_fit_FE, trt_ref = 1))
plot(arm_releff_FE, ref_line = 0)
(arm_releff_RE <- relative_effects(arm_fit_RE, trt_ref = 1))
plot(arm_releff_RE, ref_line = 0)
Following @TSD2, we produce absolute predictions of the mean off-time reduction on each treatment assuming a Normal distribution for the outcomes on treatment 1 (placebo) with mean $-0.73$ and precision $21$.
We use the predict() method, where the baseline argument takes a distr() distribution object with which we specify the corresponding Normal distribution, and we specify baseline_trt = 1 to indicate that the baseline distribution corresponds to treatment 1.
(Strictly speaking, type = "response" is unnecessary here, since the identity link function was used.)
arm_pred_FE <- predict(arm_fit_FE,
baseline = distr(qnorm, mean = -0.73, sd = 21^-0.5),
type = "response",
baseline_trt = 1)
arm_pred_FE
plot(arm_pred_FE)
arm_pred_RE <- predict(arm_fit_RE,
baseline = distr(qnorm, mean = -0.73, sd = 21^-0.5),
type = "response",
baseline_trt = 1)
arm_pred_RE
plot(arm_pred_RE)
If the baseline argument is omitted, predictions of mean off-time reduction will be produced for every study in the network based on their estimated baseline response $\mu_j$:
arm_pred_FE_studies <- predict(arm_fit_FE, type = "response")
arm_pred_FE_studies
plot(arm_pred_FE_studies)
We can also produce treatment rankings, rank probabilities, and cumulative rank probabilities.
(arm_ranks <- posterior_ranks(arm_fit_FE))
plot(arm_ranks)
(arm_rankprobs <- posterior_rank_probs(arm_fit_FE))
plot(arm_rankprobs)
(arm_cumrankprobs <- posterior_rank_probs(arm_fit_FE, cumulative = TRUE))
plot(arm_cumrankprobs)
Analysis of contrast-based data
We now perform an analysis using the contrast-based data (mean differences and standard errors).
Setting up the network
With contrast-level data giving the mean difference in off-time reduction (diff) and standard error (se_diff), we use the function set_agd_contrast() to set up the network.
contr_net <- set_agd_contrast(parkinsons,
study = studyn,
trt = trtn,
y = diff,
se = se_diff,
sample_size = n)
contr_net
The sample_size argument is optional, but enables the nodes to be weighted by sample size in the network plot.
Plot the network structure.
plot(contr_net, weight_edges = TRUE, weight_nodes = TRUE)
Meta-analysis models
We fit both fixed effect (FE) and random effects (RE) models.
Fixed effect meta-analysis
First, we fit a fixed effect model using the nma() function with trt_effects = "fixed".
We use $\mathrm{N}(0, 100^2)$ prior distributions for the treatment effects $d_k$.
We can examine the range of parameter values implied by these prior distributions with the summary() method:
summary(normal(scale = 100))
The model is fitted using the nma() function.
contr_fit_FE <- nma(contr_net,
trt_effects = "fixed",
prior_trt = normal(scale = 100))
contr_fit_FE <- nma(contr_net,
trt_effects = "fixed",
prior_trt = normal(scale = 100),
iter = 10000)
Basic parameter summaries are given by the print() method:
contr_fit_FE
The prior and posterior distributions can be compared visually using the plot_prior_posterior() function:
plot_prior_posterior(contr_fit_FE)
Random effects meta-analysis
We now fit a random effects model using the nma() function with trt_effects = "random".
Again, we use $\mathrm{N}(0, 100^2)$ prior distributions for the treatment effects $d_k$, and we additionally use a $\textrm{half-N}(5^2)$ prior for the heterogeneity standard deviation $\tau$.
We can examine the range of parameter values implied by these prior distributions with the summary() method:
summary(normal(scale = 100))
summary(half_normal(scale = 5))
Fitting the RE model
contr_fit_RE <- nma(contr_net,
seed = 1150676438,
trt_effects = "random",
prior_trt = normal(scale = 100),
prior_het = half_normal(scale = 5),
adapt_delta = 0.99)
contr_fit_RE <- nowarn_on_ci(nma(contr_net,
seed = 1150676438,
trt_effects = "random",
prior_trt = normal(scale = 100),
prior_het = half_normal(scale = 5),
adapt_delta = 0.99,
iter = 10000))
We do see a small number of divergent transition errors, which cannot simply be removed by increasing the value of the adapt_delta argument (by default set to 0.95 for RE models).
To diagnose, we use the pairs() method to investigate where in the posterior distribution these divergences are happening (indicated by red crosses):
pairs(contr_fit_RE, pars = c("d[3]", "delta[4: 4 vs. 3]", "tau"))
The divergent transitions occur in the upper tail of the heterogeneity standard deviation.
In this case, with only a small number of studies, there is not very much information to estimate the heterogeneity standard deviation and the prior distribution may be too heavy-tailed.
We could consider a more informative prior distribution for the heterogeneity variance to aid estimation.
Basic parameter summaries are given by the print() method:
contr_fit_RE
By default, summaries of the study-specific relative effects $\delta_{jk}$ are hidden, but could be examined by changing the pars argument:
# Not run
print(contr_fit_RE, pars = c("d", "delta"))
The prior and posterior distributions can be compared visually using the plot_prior_posterior() function:
plot_prior_posterior(contr_fit_RE)
Model comparison
Model fit can be checked using the dic() function:
(contr_dic_FE <- dic(contr_fit_FE))
(contr_dic_RE <- dic(contr_fit_RE))
Both models fit the data well, having posterior mean residual deviance close to the number of data points.
The DIC is similar between models, so we choose the FE model based on parsimony.
We can also examine the residual deviance contributions with the corresponding plot() method.
plot(contr_dic_FE)
plot(contr_dic_RE)
Further results
For comparison with @TSD2, we can produce relative effects against placebo using the relative_effects() function with trt_ref = 1:
(contr_releff_FE <- relative_effects(contr_fit_FE, trt_ref = 1))
plot(contr_releff_FE, ref_line = 0)
(contr_releff_RE <- relative_effects(contr_fit_RE, trt_ref = 1))
plot(contr_releff_RE, ref_line = 0)
Following @TSD2, we produce absolute predictions of the mean off-time reduction on each treatment assuming a Normal distribution for the outcomes on treatment 1 (placebo) with mean $-0.73$ and precision $21$.
We use the predict() method, where the baseline argument takes a distr() distribution object with which we specify the corresponding Normal distribution, and we specify baseline_trt = 1 to indicate that the baseline distribution corresponds to treatment 1.
(Strictly speaking, type = "response" is unnecessary here, since the identity link function was used.)
contr_pred_FE <- predict(contr_fit_FE,
baseline = distr(qnorm, mean = -0.73, sd = 21^-0.5),
type = "response",
baseline_trt = 1)
contr_pred_FE
plot(contr_pred_FE)
contr_pred_RE <- predict(contr_fit_RE,
baseline = distr(qnorm, mean = -0.73, sd = 21^-0.5),
type = "response",
baseline_trt = 1)
contr_pred_RE
plot(contr_pred_RE)
If the baseline argument is omitted an error will be raised, as there are no study baselines estimated in this network.
# Not run
predict(contr_fit_FE, type = "response")
We can also produce treatment rankings, rank probabilities, and cumulative rank probabilities.
(contr_ranks <- posterior_ranks(contr_fit_FE))
plot(contr_ranks)
(contr_rankprobs <- posterior_rank_probs(contr_fit_FE))
plot(contr_rankprobs)
(contr_cumrankprobs <- posterior_rank_probs(contr_fit_FE, cumulative = TRUE))
plot(contr_cumrankprobs)
Analysis of mixed arm-based and contrast-based data
We now perform an analysis where some studies contribute arm-based data, and other contribute contrast-based data.
Replicating @TSD2, we consider arm-based data from studies 1-3, and contrast-based data from studies 4-7.
studies <- parkinsons$studyn
(parkinsons_arm <- parkinsons[studies %in% 1:3, ])
(parkinsons_contr <- parkinsons[studies %in% 4:7, ])
Setting up the network
We use the functions set_agd_arm() and set_agd_contrast() to set up the respective data sources within the network, and then combine together with combine_network().
mix_arm_net <- set_agd_arm(parkinsons_arm,
study = studyn,
trt = trtn,
y = y,
se = se,
sample_size = n)
mix_contr_net <- set_agd_contrast(parkinsons_contr,
study = studyn,
trt = trtn,
y = diff,
se = se_diff,
sample_size = n)
mix_net <- combine_network(mix_arm_net, mix_contr_net)
mix_net
The sample_size argument is optional, but enables the nodes to be weighted by sample size in the network plot.
Plot the network structure.
plot(mix_net, weight_edges = TRUE, weight_nodes = TRUE)
Meta-analysis models
We fit both fixed effect (FE) and random effects (RE) models.
Fixed effect meta-analysis
First, we fit a fixed effect model using the nma() function with trt_effects = "fixed".
We use $\mathrm{N}(0, 100^2)$ prior distributions for the treatment effects $d_k$ and study-specific intercepts $\mu_j$.
We can examine the range of parameter values implied by these prior distributions with the summary() method:
summary(normal(scale = 100))
The model is fitted using the nma() function.
mix_fit_FE <- nma(mix_net,
trt_effects = "fixed",
prior_intercept = normal(scale = 100),
prior_trt = normal(scale = 100))
mix_fit_FE <- nma(mix_net,
trt_effects = "fixed",
prior_intercept = normal(scale = 100),
prior_trt = normal(scale = 100),
iter = 10000)
Basic parameter summaries are given by the print() method:
mix_fit_FE
By default, summaries of the study-specific intercepts $\mu_j$ are hidden, but could be examined by changing the pars argument:
# Not run
print(mix_fit_FE, pars = c("d", "mu"))
The prior and posterior distributions can be compared visually using the plot_prior_posterior() function:
plot_prior_posterior(mix_fit_FE)
Random effects meta-analysis
We now fit a random effects model using the nma() function with trt_effects = "random".
Again, we use $\mathrm{N}(0, 100^2)$ prior distributions for the treatment effects $d_k$ and study-specific intercepts $\mu_j$, and we additionally use a $\textrm{half-N}(5^2)$ prior for the heterogeneity standard deviation $\tau$.
We can examine the range of parameter values implied by these prior distributions with the summary() method:
summary(normal(scale = 100))
summary(half_normal(scale = 5))
Fitting the RE model
mix_fit_RE <- nma(mix_net,
seed = 437219664,
trt_effects = "random",
prior_intercept = normal(scale = 100),
prior_trt = normal(scale = 100),
prior_het = half_normal(scale = 5),
adapt_delta = 0.99)
mix_fit_RE <- nowarn_on_ci(nma(mix_net,
seed = 437219664,
trt_effects = "random",
prior_intercept = normal(scale = 100),
prior_trt = normal(scale = 100),
prior_het = half_normal(scale = 5),
adapt_delta = 0.99,
iter=10000))
We do see a small number of divergent transition errors, which cannot simply be removed by increasing the value of the adapt_delta argument (by default set to 0.95 for RE models).
To diagnose, we use the pairs() method to investigate where in the posterior distribution these divergences are happening (indicated by red crosses):
pairs(mix_fit_RE, pars = c("d[3]", "delta[4: 4 vs. 3]", "tau"))
The divergent transitions occur in the upper tail of the heterogeneity standard deviation.
In this case, with only a small number of studies, there is not very much information to estimate the heterogeneity standard deviation and the prior distribution may be too heavy-tailed.
We could consider a more informative prior distribution for the heterogeneity variance to aid estimation.
Basic parameter summaries are given by the print() method:
mix_fit_RE
By default, summaries of the study-specific intercepts $\mu_j$ and study-specific relative effects $\delta_{jk}$ are hidden, but could be examined by changing the pars argument:
# Not run
print(mix_fit_RE, pars = c("d", "mu", "delta"))
The prior and posterior distributions can be compared visually using the plot_prior_posterior() function:
plot_prior_posterior(mix_fit_RE)
Model comparison
Model fit can be checked using the dic() function:
(mix_dic_FE <- dic(mix_fit_FE))
(mix_dic_RE <- dic(mix_fit_RE))
Both models fit the data well, having posterior mean residual deviance close to the number of data points.
The DIC is similar between models, so we choose the FE model based on parsimony.
We can also examine the residual deviance contributions with the corresponding plot() method.
plot(mix_dic_FE)
plot(mix_dic_RE)
Further results
For comparison with @TSD2, we can produce relative effects against placebo using the relative_effects() function with trt_ref = 1:
(mix_releff_FE <- relative_effects(mix_fit_FE, trt_ref = 1))
plot(mix_releff_FE, ref_line = 0)
(mix_releff_RE <- relative_effects(mix_fit_RE, trt_ref = 1))
plot(mix_releff_RE, ref_line = 0)
Following @TSD2, we produce absolute predictions of the mean off-time reduction on each treatment assuming a Normal distribution for the outcomes on treatment 1 (placebo) with mean $-0.73$ and precision $21$.
We use the predict() method, where the baseline argument takes a distr() distribution object with which we specify the corresponding Normal distribution, and we specify baseline_trt = 1 to indicate that the baseline distribution corresponds to treatment 1.
(Strictly speaking, type = "response" is unnecessary here, since the identity link function was used.)
mix_pred_FE <- predict(mix_fit_FE,
baseline = distr(qnorm, mean = -0.73, sd = 21^-0.5),
type = "response",
baseline_trt = 1)
mix_pred_FE
plot(mix_pred_FE)
mix_pred_RE <- predict(mix_fit_RE,
baseline = distr(qnorm, mean = -0.73, sd = 21^-0.5),
type = "response",
baseline_trt = 1)
mix_pred_RE
plot(mix_pred_RE)
If the baseline argument is omitted, predictions of mean off-time reduction will be produced for every arm-based study in the network based on their estimated baseline response $\mu_j$:
mix_pred_FE_studies <- predict(mix_fit_FE, type = "response")
mix_pred_FE_studies
plot(mix_pred_FE_studies)
We can also produce treatment rankings, rank probabilities, and cumulative rank probabilities.
(mix_ranks <- posterior_ranks(mix_fit_FE))
plot(mix_ranks)
(mix_rankprobs <- posterior_rank_probs(mix_fit_FE))
plot(mix_rankprobs)
(mix_cumrankprobs <- posterior_rank_probs(mix_fit_FE, cumulative = TRUE))
plot(mix_cumrankprobs)
References
#--- Test against TSD 2 results ---
library(testthat)
library(dplyr)
tol <- 0.05
tol_dic <- 0.1
#--- Arm-level data ---
# FE Relative effects
tsd_FE <- tribble(
~trt, ~mean, ~sd , ~median, ~lower, ~upper,
2 , -1.81, 0.33, -1.81 , -2.46 ,-1.16 ,
3 , -0.47, 0.49, -0.47 , -1.43 ,0.49 ,
4 , -0.52, 0.48, -0.52 , -1.46 ,0.43 ,
5 , -0.82, 0.52, -0.82 , -1.84 ,0.22 ,
) %>%
mutate(trt = ordered(trt, levels = levels(arm_net$treatments))) %>%
arrange(trt)
arm_releff_FE <- as.data.frame(arm_releff_FE)
test_that("FE relative effects", {
expect_equivalent(arm_releff_FE$mean, tsd_FE$mean, tolerance = tol)
expect_equivalent(arm_releff_FE$sd, tsd_FE$sd, tolerance = tol)
expect_equivalent(arm_releff_FE$`50%`, tsd_FE$median, tolerance = tol)
expect_equivalent(arm_releff_FE$`2.5%`, tsd_FE$lower, tolerance = tol)
expect_equivalent(arm_releff_FE$`97.5%`, tsd_FE$upper, tolerance = tol)
})
# RE Relative effects
tsd_RE <- tribble(
~trt, ~mean, ~sd , ~median, ~lower, ~upper,
2 , -1.85, 0.54, -1.84 , -2.91 ,-0.85 ,
3 , -0.50, 0.66, -0.50 , -1.78 ,0.75 ,
4 , -0.53, 0.65, -0.53 , -1.77 ,0.71 ,
5 , -0.83, 0.80, -0.83 , -2.35 ,0.69 ,
) %>%
mutate(trt = ordered(trt, levels = levels(arm_net$treatments))) %>%
arrange(trt)
arm_releff_RE <- as.data.frame(arm_releff_RE)
test_that("RE relative effects", {
expect_equivalent(arm_releff_RE$mean, tsd_RE$mean, tolerance = tol)
expect_equivalent(arm_releff_RE$`50%`, tsd_RE$median, tolerance = tol)
skip_on_ci()
expect_equivalent(arm_releff_RE$sd, tsd_RE$sd, tolerance = tol)
expect_equivalent(arm_releff_RE$`2.5%`, tsd_RE$lower, tolerance = tol)
expect_equivalent(arm_releff_RE$`97.5%`, tsd_RE$upper, tolerance = tol)
})
# Heterogeneity SD
arm_tau <- as.data.frame(summary(arm_fit_RE, pars = "tau"))
test_that("RE heterogeneity SD", {
expect_equivalent(arm_tau$mean, 0.4, tolerance = tol)
expect_equivalent(arm_tau$`50%`, 0.28, tolerance = tol)
skip_on_ci()
expect_equivalent(arm_tau$sd, 0.43, tolerance = tol)
expect_equivalent(arm_tau$`2.5%`, 0.01, tolerance = tol)
expect_equivalent(arm_tau$`97.5%`, 1.55, tolerance = tol)
})
# FE probabilities
tsd_pred_FE <- tribble(
~trt, ~mean, ~sd , ~median, ~lower, ~upper,
1 , -0.73, 0.22, -0.73 , -1.16 ,-0.30 ,
2 , -2.54, 0.40, -2.54 , -3.32 ,-1.76 ,
3 , -1.21, 0.53, -1.20 , -2.25 ,-0.15 ,
4 , -1.25, 0.53, -1.25 , -2.28 ,-0.21 ,
5 , -1.55, 0.57, -1.55 , -2.66 ,-0.43 ,
) %>%
mutate(trt = ordered(trt, levels = levels(arm_net$treatments))) %>%
arrange(trt)
arm_pred_FE <- as.data.frame(arm_pred_FE)
test_that("FE predicted probabilities at 3 years", {
expect_equivalent(arm_pred_FE$mean, tsd_pred_FE$mean, tolerance = tol)
expect_equivalent(arm_pred_FE$sd, tsd_pred_FE$sd, tolerance = tol)
expect_equivalent(arm_pred_FE$`50%`, tsd_pred_FE$median, tolerance = tol)
expect_equivalent(arm_pred_FE$`2.5%`, tsd_pred_FE$lower, tolerance = tol)
expect_equivalent(arm_pred_FE$`97.5%`, tsd_pred_FE$upper, tolerance = tol)
})
# RE probabilities
tsd_pred_RE <- tribble(
~trt, ~mean, ~sd , ~median, ~lower, ~upper,
1 , -0.73, 0.22, -0.73 , -1.16 ,-0.30 ,
2 , -2.58, 0.58, -2.57 , -3.72 ,-1.50 ,
3 , -1.23, 0.70, -1.23 , -2.57 ,0.10 ,
4 , -1.26, 0.69, -1.26 , -2.57 ,0.05 ,
5 , -1.57, 0.83, -1.56 , -3.14 ,0.02 ,
) %>%
mutate(trt = ordered(trt, levels = levels(arm_net$treatments))) %>%
arrange(trt)
arm_pred_RE <- as.data.frame(arm_pred_RE)
test_that("RE predicted probabilities at 3 years", {
expect_equivalent(arm_pred_RE$mean, tsd_pred_RE$mean, tolerance = tol)
expect_equivalent(arm_pred_RE$`50%`, tsd_pred_RE$median, tolerance = tol)
skip_on_ci()
expect_equivalent(arm_pred_RE$sd, tsd_pred_RE$sd, tolerance = tol)
expect_equivalent(arm_pred_RE$`2.5%`, tsd_pred_RE$lower, tolerance = tol)
expect_equivalent(arm_pred_RE$`97.5%`, tsd_pred_RE$upper, tolerance = tol)
})
# FE DIC
test_that("FE DIC", {
expect_equivalent(arm_dic_FE$resdev, 13.3, tolerance = tol_dic)
expect_equivalent(arm_dic_FE$pd, 11.0, tolerance = tol_dic)
expect_equivalent(arm_dic_FE$dic, 24.3, tolerance = tol_dic)
})
# RE DIC
test_that("RE DIC", {
expect_equivalent(arm_dic_RE$resdev, 13.6, tolerance = tol_dic)
expect_equivalent(arm_dic_RE$pd, 12.4, tolerance = tol_dic)
expect_equivalent(arm_dic_RE$dic, 26.0, tolerance = tol_dic)
})
#--- Contrast-level data ---
# Relative effects, heterogeneity, absolute effects are all equal across data scenarios
contr_releff_FE <- as.data.frame(contr_releff_FE)
test_that("FE relative effects", {
expect_equivalent(contr_releff_FE$mean, tsd_FE$mean, tolerance = tol)
expect_equivalent(contr_releff_FE$sd, tsd_FE$sd, tolerance = tol)
expect_equivalent(contr_releff_FE$`50%`, tsd_FE$median, tolerance = tol)
expect_equivalent(contr_releff_FE$`2.5%`, tsd_FE$lower, tolerance = tol)
expect_equivalent(contr_releff_FE$`97.5%`, tsd_FE$upper, tolerance = tol)
})
contr_releff_RE <- as.data.frame(contr_releff_RE)
test_that("RE relative effects", {
expect_equivalent(contr_releff_RE$mean, tsd_RE$mean, tolerance = tol)
expect_equivalent(contr_releff_RE$`50%`, tsd_RE$median, tolerance = tol)
skip_on_ci()
expect_equivalent(contr_releff_RE$sd, tsd_RE$sd, tolerance = tol)
expect_equivalent(contr_releff_RE$`2.5%`, tsd_RE$lower, tolerance = tol)
expect_equivalent(contr_releff_RE$`97.5%`, tsd_RE$upper, tolerance = tol)
})
# Heterogeneity SD
contr_tau <- as.data.frame(summary(contr_fit_RE, pars = "tau"))
test_that("RE heterogeneity SD", {
expect_equivalent(contr_tau$mean, 0.4, tolerance = tol)
expect_equivalent(contr_tau$`50%`, 0.28, tolerance = tol)
skip_on_ci()
expect_equivalent(contr_tau$sd, 0.43, tolerance = tol)
expect_equivalent(contr_tau$`2.5%`, 0.01, tolerance = tol)
expect_equivalent(contr_tau$`97.5%`, 1.55, tolerance = tol)
})
contr_pred_FE <- as.data.frame(contr_pred_FE)
test_that("FE predicted probabilities at 3 years", {
expect_equivalent(contr_pred_FE$mean, tsd_pred_FE$mean, tolerance = tol)
expect_equivalent(contr_pred_FE$sd, tsd_pred_FE$sd, tolerance = tol)
expect_equivalent(contr_pred_FE$`50%`, tsd_pred_FE$median, tolerance = tol)
expect_equivalent(contr_pred_FE$`2.5%`, tsd_pred_FE$lower, tolerance = tol)
expect_equivalent(contr_pred_FE$`97.5%`, tsd_pred_FE$upper, tolerance = tol)
})
contr_pred_RE <- as.data.frame(contr_pred_RE)
test_that("RE predicted probabilities at 3 years", {
expect_equivalent(contr_pred_RE$mean, tsd_pred_RE$mean, tolerance = tol)
expect_equivalent(contr_pred_RE$`50%`, tsd_pred_RE$median, tolerance = tol)
skip_on_ci()
expect_equivalent(contr_pred_RE$sd, tsd_pred_RE$sd, tolerance = tol)
expect_equivalent(contr_pred_RE$`2.5%`, tsd_pred_RE$lower, tolerance = tol)
expect_equivalent(contr_pred_RE$`97.5%`, tsd_pred_RE$upper, tolerance = tol)
})
test_that("Error message when using predict() on only contrast data", {
expect_error(predict(contr_fit_FE),
"No arm-based data \\(IPD or AgD\\) in network")
})
# FE DIC
test_that("FE DIC", {
expect_equivalent(contr_dic_FE$resdev, 6.3, tolerance = tol_dic)
expect_equivalent(contr_dic_FE$pd, 4.0, tolerance = tol_dic)
expect_equivalent(contr_dic_FE$dic, 10.3, tolerance = tol_dic)
})
# RE DIC
test_that("RE DIC", {
expect_equivalent(contr_dic_RE$resdev, 6.6, tolerance = tol_dic)
expect_equivalent(contr_dic_RE$pd, 5.5, tolerance = tol_dic)
expect_equivalent(contr_dic_RE$dic, 12.1, tolerance = tol_dic)
})
# Test contrast-data rows in incorrect order
park_reorder <- sample_frac(parkinsons, size = 1, replace = FALSE)
reorder_net <- set_agd_contrast(park_reorder,
study = studyn,
trt = trtn,
y = diff,
se = se_diff,
sample_size = n)
reorder_fit_FE <- nma(reorder_net,
seed = 1878819627,
trt_effects = "fixed",
prior_trt = normal(scale = 100),
iter = 10000)
reorder_fit_RE <- nowarn_on_ci(nma(reorder_net,
seed = 1147315,
trt_effects = "random",
prior_trt = normal(scale = 100),
prior_het = half_normal(scale = 5),
adapt_delta = 0.99,
iter = 10000))
expect_equivalent_nma_summary <- function(object, expected, ...) {
act <- quasi_label(enquo(object), arg = "object")
exp <- quasi_label(enquo(expected), arg = "expected")
comp <- compare(act$val %>% as.data.frame() %>% select(-c(Bulk_ESS, Tail_ESS, Rhat)),
exp$val %>% as.data.frame() %>% select(-c(Bulk_ESS, Tail_ESS, Rhat)),
...)
expect(comp$equal, sprintf("%s not equal to %s.\n%s",
act$lab, exp$lab, comp$message))
invisible(act$val)
}
test_that("Analysis of reordered contrast data gives same results", {
expect_equivalent_nma_summary(relative_effects(reorder_fit_FE, trt_ref = 1),
contr_releff_FE,
tolerance = tol)
expect_equivalent_nma_summary(predict(reorder_fit_FE,
baseline = distr(qnorm, mean = -0.73, sd = 21^-0.5),
type = "response",
baseline_trt = 1),
contr_pred_FE,
tolerance = tol)
expect_equivalent(dic(reorder_fit_FE)$resdev, contr_dic_FE$resdev, tolerance = tol_dic)
expect_equivalent(dic(reorder_fit_FE)$pd, contr_dic_FE$pd, tolerance = tol_dic)
expect_equivalent(dic(reorder_fit_FE)$dic, contr_dic_FE$dic, tolerance = tol_dic)
expect_equivalent(dic(reorder_fit_RE)$resdev, contr_dic_RE$resdev, tolerance = tol_dic)
expect_equivalent(dic(reorder_fit_RE)$pd, contr_dic_RE$pd, tolerance = tol_dic)
expect_equivalent(dic(reorder_fit_RE)$dic, contr_dic_RE$dic, tolerance = tol_dic)
skip_on_ci()
expect_equivalent_nma_summary(relative_effects(reorder_fit_RE, trt_ref = 1),
contr_releff_RE,
tolerance = tol)
expect_equivalent_nma_summary(predict(reorder_fit_RE,
baseline = distr(qnorm, mean = -0.73, sd = 21^-0.5),
type = "response",
baseline_trt = 1),
contr_pred_RE,
tolerance = tol)
expect_equivalent_nma_summary(summary(reorder_fit_RE, pars = "tau"),
contr_tau,
tolerance = tol)
})
#--- Mixed data ---
# Relative effects, heterogeneity, absolute effects are all equal across data scenarios
mix_releff_FE <- as.data.frame(mix_releff_FE)
test_that("FE relative effects", {
expect_equivalent(mix_releff_FE$mean, tsd_FE$mean, tolerance = tol)
expect_equivalent(mix_releff_FE$sd, tsd_FE$sd, tolerance = tol)
expect_equivalent(mix_releff_FE$`50%`, tsd_FE$median, tolerance = tol)
expect_equivalent(mix_releff_FE$`2.5%`, tsd_FE$lower, tolerance = tol)
expect_equivalent(mix_releff_FE$`97.5%`, tsd_FE$upper, tolerance = tol)
})
mix_releff_RE <- as.data.frame(mix_releff_RE)
test_that("RE relative effects", {
expect_equivalent(mix_releff_RE$mean, tsd_RE$mean, tolerance = tol)
expect_equivalent(mix_releff_RE$`50%`, tsd_RE$median, tolerance = tol)
skip_on_ci()
expect_equivalent(mix_releff_RE$sd, tsd_RE$sd, tolerance = tol)
expect_equivalent(mix_releff_RE$`2.5%`, tsd_RE$lower, tolerance = tol)
expect_equivalent(mix_releff_RE$`97.5%`, tsd_RE$upper, tolerance = tol)
})
# Heterogeneity SD
mix_tau <- as.data.frame(summary(mix_fit_RE, pars = "tau"))
test_that("RE heterogeneity SD", {
expect_equivalent(mix_tau$mean, 0.4, tolerance = tol)
expect_equivalent(mix_tau$`50%`, 0.28, tolerance = tol)
skip_on_ci()
expect_equivalent(mix_tau$sd, 0.43, tolerance = tol)
expect_equivalent(mix_tau$`2.5%`, 0.01, tolerance = tol)
expect_equivalent(mix_tau$`97.5%`, 1.55, tolerance = tol)
})
mix_pred_FE <- as.data.frame(mix_pred_FE)
test_that("FE predicted probabilities at 3 years", {
expect_equivalent(mix_pred_FE$mean, tsd_pred_FE$mean, tolerance = tol)
expect_equivalent(mix_pred_FE$sd, tsd_pred_FE$sd, tolerance = tol)
expect_equivalent(mix_pred_FE$`50%`, tsd_pred_FE$median, tolerance = tol)
expect_equivalent(mix_pred_FE$`2.5%`, tsd_pred_FE$lower, tolerance = tol)
expect_equivalent(mix_pred_FE$`97.5%`, tsd_pred_FE$upper, tolerance = tol)
})
mix_pred_RE <- as.data.frame(mix_pred_RE)
test_that("RE predicted probabilities at 3 years", {
expect_equivalent(mix_pred_RE$mean, tsd_pred_RE$mean, tolerance = tol)
expect_equivalent(mix_pred_RE$`50%`, tsd_pred_RE$median, tolerance = tol)
skip_on_ci()
expect_equivalent(mix_pred_RE$sd, tsd_pred_RE$sd, tolerance = tol)
expect_equivalent(mix_pred_RE$`2.5%`, tsd_pred_RE$lower, tolerance = tol)
expect_equivalent(mix_pred_RE$`97.5%`, tsd_pred_RE$upper, tolerance = tol)
})
# FE DIC
test_that("FE DIC", {
expect_equivalent(mix_dic_FE$resdev, 9.3, tolerance = tol_dic)
expect_equivalent(mix_dic_FE$pd, 7.0, tolerance = tol_dic)
expect_equivalent(mix_dic_FE$dic, 16.3, tolerance = tol_dic)
})
# RE DIC
test_that("RE DIC", {
expect_equivalent(mix_dic_RE$resdev, 9.6, tolerance = tol_dic)
expect_equivalent(mix_dic_RE$pd, 8.5, tolerance = tol_dic)
expect_equivalent(mix_dic_RE$dic, 18.1, tolerance = tol_dic)
})
# --- Check UME models ---
arm_fit_FE_UME <- nma(arm_net,
trt_effects = "fixed",
consistency = "ume",
prior_intercept = normal(scale = 100),
prior_trt = normal(scale = 10),
iter = 10000)
# NOTE: RE UME models are not really tractable - lots of divergences - not
# enough data to estimate tau
mix_fit_FE_UME <- nma(mix_net,
trt_effects = "fixed",
consistency = "ume",
prior_intercept = normal(scale = 100),
prior_trt = normal(scale = 100),
iter = 10000)
test_that("Arm/Mixed UME results are the same", {
expect_equivalent_nma_summary(summary(arm_fit_FE_UME, pars = "d"),
summary(mix_fit_FE_UME, pars = "d"),
tolerance = tol)
})
# NOTE: Contrast (and reordered contrast) results are different, since the
# choice of independent contrasts is based on the "observed" baseline arms in
# the data rather than the treatment ordering
contr_fit_FE_UME <- nma(contr_net,
trt_effects = "fixed",
consistency = "ume",
prior_trt = normal(scale = 100),
iter = 10000)
reorder_fit_FE_UME <- nma(reorder_net,
seed = 1878819627,
trt_effects = "fixed",
consistency = "ume",
prior_trt = normal(scale = 100),
iter = 10000)
test_that("Contr/Reordered UME results are the same", {
expect_equivalent_nma_summary(summary(contr_fit_FE_UME, pars = "d"),
summary(reorder_fit_FE_UME, pars = "d"),
tolerance = tol)
})
# --- Test regression models with contrast data ---
park_dummy <- dplyr::mutate(parkinsons, x1 = rnorm(dplyr::n(), 0, 1))
arm_net_reg <- set_agd_arm(park_dummy,
study = studyn,
trt = trtn,
y = y,
se = se,
sample_size = n)
arm_fit_FE_reg <- nma(arm_net_reg,
trt_effects = "fixed",
regression = ~x1:.trt,
prior_intercept = normal(scale = 100),
prior_trt = normal(scale = 10),
prior_reg = normal(scale = 1),
iter = 10000)
contr_net_reg <- set_agd_contrast(park_dummy,
study = studyn,
trt = trtn,
y = diff,
se = se_diff,
sample_size = n)
contr_fit_FE_reg <- nma(contr_net_reg,
trt_effects = "fixed",
regression = ~x1:.trt,
prior_trt = normal(scale = 100),
prior_reg = normal(scale = 1),
iter = 10000)
mix_net_reg <- combine_network(set_agd_arm(park_dummy[studies %in% 1:3, ],
study = studyn,
trt = trtn,
y = y,
se = se,
sample_size = n),
set_agd_contrast(park_dummy[studies %in% 4:7, ],
study = studyn,
trt = trtn,
y = diff,
se = se_diff,
sample_size = n))
mix_fit_FE_reg <- nma(mix_net_reg,
trt_effects = "fixed",
regression = ~x1:.trt,
prior_intercept = normal(scale = 100),
prior_trt = normal(scale = 100),
prior_reg = normal(scale = 1),
iter = 10000)
reorder_net_reg <- set_agd_contrast(sample_frac(park_dummy, size = 1, replace = FALSE),
study = studyn,
trt = trtn,
y = diff,
se = se_diff,
sample_size = n)
reorder_fit_FE_reg <- nma(reorder_net_reg,
seed = 1878819627,
trt_effects = "fixed",
regression = ~x1:.trt,
prior_trt = normal(scale = 100),
prior_reg = normal(scale = 1),
iter = 10000)
test_that("Regression models work with contrast data", {
expect_equivalent_nma_summary(summary(arm_fit_FE_reg, pars = c("d", "beta")),
summary(contr_fit_FE_reg, pars = c("d", "beta")),
tolerance = tol)
expect_equivalent_nma_summary(summary(arm_fit_FE_reg, pars = c("d", "beta")),
summary(mix_fit_FE_reg, pars = c("d", "beta")),
tolerance = tol)
expect_equivalent_nma_summary(summary(contr_fit_FE_reg, pars = c("d", "beta")),
summary(reorder_fit_FE_reg, pars = c("d", "beta")),
tolerance = tol)
})
test_that("Error message when using predict() on only contrast data", {
expect_error(predict(contr_fit_FE_reg),
"No arm-based data \\(IPD or AgD\\) in network")
})
# ML-NMR regression models too
park_dummy_mlnmr <- dplyr::mutate(parkinsons,
x1_mean = rnorm(dplyr::n(), 0, 1),
x1_sd = runif(dplyr::n(), 0.1, 1))
arm_net_mlnmr <- set_agd_arm(park_dummy_mlnmr,
study = studyn,
trt = trtn,
y = y,
se = se,
sample_size = n) %>%
add_integration(x1 = distr(qnorm, x1_mean, x1_sd), n_int = 10)
arm_fit_FE_mlnmr <- nma(arm_net_mlnmr,
trt_effects = "fixed",
regression = ~x1:.trt,
prior_intercept = normal(scale = 100),
prior_trt = normal(scale = 10),
prior_reg = normal(scale = 1),
iter = 20000, save_warmup = FALSE)
contr_net_mlnmr <- set_agd_contrast(park_dummy_mlnmr,
study = studyn,
trt = trtn,
y = diff,
se = se_diff,
sample_size = n) %>%
add_integration(x1 = distr(qnorm, x1_mean, x1_sd), n_int = 10)
contr_fit_FE_mlnmr <- nma(contr_net_mlnmr,
trt_effects = "fixed",
regression = ~x1:.trt,
prior_trt = normal(scale = 100),
prior_reg = normal(scale = 1),
iter = 20000, save_warmup = FALSE)
mix_net_mlnmr <- combine_network(set_agd_arm(park_dummy_mlnmr[studies %in% 1:3, ],
study = studyn,
trt = trtn,
y = y,
se = se,
sample_size = n),
set_agd_contrast(park_dummy_mlnmr[studies %in% 4:7, ],
study = studyn,
trt = trtn,
y = diff,
se = se_diff,
sample_size = n)) %>%
add_integration(x1 = distr(qnorm, x1_mean, x1_sd), n_int = 10)
mix_fit_FE_mlnmr <- nma(mix_net_mlnmr,
trt_effects = "fixed",
regression = ~x1:.trt,
prior_intercept = normal(scale = 100),
prior_trt = normal(scale = 100),
prior_reg = normal(scale = 1),
iter = 20000, save_warmup = FALSE)
reorder_net_mlnmr <- set_agd_contrast(sample_frac(park_dummy_mlnmr, size = 1, replace = FALSE),
study = studyn,
trt = trtn,
y = diff,
se = se_diff,
sample_size = n) %>%
add_integration(x1 = distr(qnorm, x1_mean, x1_sd), n_int = 10)
reorder_fit_FE_mlnmr <- nma(reorder_net_mlnmr,
seed = 1878819627,
trt_effects = "fixed",
regression = ~x1:.trt,
prior_trt = normal(scale = 100),
prior_reg = normal(scale = 1),
iter = 20000, save_warmup = FALSE)
test_that("ML-NMR models work with contrast data", {
expect_equivalent_nma_summary(summary(arm_fit_FE_mlnmr, pars = c("d", "beta")),
summary(contr_fit_FE_mlnmr, pars = c("d", "beta")),
tolerance = tol)
expect_equivalent_nma_summary(summary(arm_fit_FE_mlnmr, pars = c("d", "beta")),
summary(mix_fit_FE_mlnmr, pars = c("d", "beta")),
tolerance = tol)
expect_equivalent_nma_summary(summary(contr_fit_FE_mlnmr, pars = c("d", "beta")),
summary(reorder_fit_FE_mlnmr, pars = c("d", "beta")),
tolerance = tol)
})
test_that("Error message when using predict() on only contrast data", {
expect_error(predict(contr_fit_FE_mlnmr),
"No arm-based data \\(IPD or AgD\\) in network")
})
test_that("Robust to custom options(contrasts) settings", {
arm_fit_FE_SAS <- withr::with_options(list(contrasts = c(ordered = "contr.SAS",
unordered = "contr.SAS")),
nma(arm_net,
trt_effects = "fixed",
prior_intercept = normal(scale = 100),
prior_trt = normal(scale = 10),
iter = 10000))
expect_equivalent_nma_summary(summary(arm_fit_FE_SAS), summary(arm_fit_FE),
tolerance = tol)
expect_equivalent_nma_summary(relative_effects(arm_fit_FE_SAS), relative_effects(arm_fit_FE),
tolerance = tol)
contr_fit_FE_SAS <- withr::with_options(list(contrasts = c(ordered = "contr.SAS",
unordered = "contr.SAS")),
nma(contr_net,
trt_effects = "fixed",
prior_intercept = normal(scale = 100),
prior_trt = normal(scale = 10),
iter = 10000))
expect_equivalent_nma_summary(summary(contr_fit_FE_SAS), summary(contr_fit_FE),
tolerance = tol)
expect_equivalent_nma_summary(relative_effects(contr_fit_FE_SAS), relative_effects(contr_fit_FE),
tolerance = tol)
mix_fit_FE_SAS <- withr::with_options(list(contrasts = c(ordered = "contr.SAS",
unordered = "contr.SAS")),
nma(mix_net,
trt_effects = "fixed",
prior_intercept = normal(scale = 100),
prior_trt = normal(scale = 10),
iter = 10000))
expect_equivalent_nma_summary(summary(mix_fit_FE_SAS), summary(mix_fit_FE),
tolerance = tol)
expect_equivalent_nma_summary(relative_effects(mix_fit_FE_SAS), relative_effects(mix_fit_FE),
tolerance = tol)
})
dmphillippo/multinma documentation built on April 12, 2025, 11:41 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
set.seed(1042020)
library(multinma) options(mc.cores = parallel::detectCores())
library(multinma) nc <- switch(tolower(Sys.getenv("_R_CHECK_LIMIT_CORES_")), "true" =, "warn" = 2, parallel::detectCores()) options(mc.cores = nc)
This vignette describes the analysis of data on the mean off-time reduction in patients given dopamine agonists as adjunct therapy in Parkinson's disease, in a network of 7 trials of 4 active drugs plus placebo [@TSD2].
The data are available in this package as parkinsons:
head(parkinsons)
We consider analysing these data in three separate ways:
- Using arm-based data (means
yand corresponding standard errorsse); - Using contrast-based data (mean differences
diffand corresponding standard errorsse_diff); - A combination of the two, where some studies contribute arm-based data, and other contribute contrast-based data.
Note: In this case, with Normal likelihoods for both arms and contrasts, we will see that the three analyses give identical results. In general, unless the arm-based likelihood is Normal, results from a model using a contrast-based likelihood will not exactly match those from a model using an arm-based likelihood, since the contrast-based Normal likelihood is only an approximation. Similarity of results depends on the suitability of the Normal approximation, which may not always be appropriate - e.g. with a small number of events or small sample size for a binary outcome. The use of an arm-based likelihood (sometimes called an "exact" likelihood) is therefore preferable where possible in general.
Analysis of arm-based data
We begin with an analysis of the arm-based data - means and standard errors.
Setting up the network
We have arm-level continuous data giving the mean off-time reduction (y) and standard error (se) in each arm.
We use the function set_agd_arm() to set up the network.
arm_net <- set_agd_arm(parkinsons, study = studyn, trt = trtn, y = y, se = se, sample_size = n) arm_net
We let treatment 4 be set by default as the network reference treatment, since this results in considerably improved sampling efficiency over choosing treatment 1 as the network reference.
The sample_size argument is optional, but enables the nodes to be weighted by sample size in the network plot.
Plot the network structure.
plot(arm_net, weight_edges = TRUE, weight_nodes = TRUE)
Meta-analysis models
We fit both fixed effect (FE) and random effects (RE) models.
Fixed effect meta-analysis
First, we fit a fixed effect model using the nma() function with trt_effects = "fixed".
We use $\mathrm{N}(0, 100^2)$ prior distributions for the treatment effects $d_k$ and study-specific intercepts $\mu_j$.
We can examine the range of parameter values implied by these prior distributions with the summary() method:
summary(normal(scale = 100))
The model is fitted using the nma() function.
arm_fit_FE <- nma(arm_net, trt_effects = "fixed", prior_intercept = normal(scale = 100), prior_trt = normal(scale = 10))
arm_fit_FE <- nma(arm_net, trt_effects = "fixed", prior_intercept = normal(scale = 100), prior_trt = normal(scale = 10), iter = 10000)
Basic parameter summaries are given by the print() method:
arm_fit_FE
By default, summaries of the study-specific intercepts $\mu_j$ are hidden, but could be examined by changing the pars argument:
# Not run print(arm_fit_FE, pars = c("d", "mu"))
The prior and posterior distributions can be compared visually using the plot_prior_posterior() function:
plot_prior_posterior(arm_fit_FE)
Random effects meta-analysis
We now fit a random effects model using the nma() function with trt_effects = "random".
Again, we use $\mathrm{N}(0, 100^2)$ prior distributions for the treatment effects $d_k$ and study-specific intercepts $\mu_j$, and we additionally use a $\textrm{half-N}(5^2)$ prior for the heterogeneity standard deviation $\tau$.
We can examine the range of parameter values implied by these prior distributions with the summary() method:
summary(normal(scale = 100)) summary(half_normal(scale = 5))
Fitting the RE model
arm_fit_RE <- nma(arm_net, seed = 379394727, trt_effects = "random", prior_intercept = normal(scale = 100), prior_trt = normal(scale = 100), prior_het = half_normal(scale = 5), adapt_delta = 0.99)
arm_fit_RE <- nowarn_on_ci(nma(arm_net, seed = 379394727, trt_effects = "random", prior_intercept = normal(scale = 100), prior_trt = normal(scale = 100), prior_het = half_normal(scale = 5), adapt_delta = 0.99, iter = 10000))
We do see a small number of divergent transition errors, which cannot simply be removed by increasing the value of the adapt_delta argument (by default set to 0.95 for RE models).
To diagnose, we use the pairs() method to investigate where in the posterior distribution these divergences are happening (indicated by red crosses):
pairs(arm_fit_RE, pars = c("mu[4]", "d[3]", "delta[4: 3]", "tau"))
The divergent transitions occur in the upper tail of the heterogeneity standard deviation. In this case, with only a small number of studies, there is not very much information to estimate the heterogeneity standard deviation and the prior distribution may be too heavy-tailed. We could consider a more informative prior distribution for the heterogeneity variance to aid estimation.
Basic parameter summaries are given by the print() method:
arm_fit_RE
By default, summaries of the study-specific intercepts $\mu_j$ and study-specific relative effects $\delta_{jk}$ are hidden, but could be examined by changing the pars argument:
# Not run print(arm_fit_RE, pars = c("d", "mu", "delta"))
The prior and posterior distributions can be compared visually using the plot_prior_posterior() function:
plot_prior_posterior(arm_fit_RE)
Model comparison
Model fit can be checked using the dic() function:
(arm_dic_FE <- dic(arm_fit_FE))
(arm_dic_RE <- dic(arm_fit_RE))
Both models fit the data well, having posterior mean residual deviance close to the number of data points. The DIC is similar between models, so we choose the FE model based on parsimony.
We can also examine the residual deviance contributions with the corresponding plot() method.
plot(arm_dic_FE)
plot(arm_dic_RE)
Further results
For comparison with @TSD2, we can produce relative effects against placebo using the relative_effects() function with trt_ref = 1:
(arm_releff_FE <- relative_effects(arm_fit_FE, trt_ref = 1)) plot(arm_releff_FE, ref_line = 0)
(arm_releff_RE <- relative_effects(arm_fit_RE, trt_ref = 1)) plot(arm_releff_RE, ref_line = 0)
Following @TSD2, we produce absolute predictions of the mean off-time reduction on each treatment assuming a Normal distribution for the outcomes on treatment 1 (placebo) with mean $-0.73$ and precision $21$.
We use the predict() method, where the baseline argument takes a distr() distribution object with which we specify the corresponding Normal distribution, and we specify baseline_trt = 1 to indicate that the baseline distribution corresponds to treatment 1.
(Strictly speaking, type = "response" is unnecessary here, since the identity link function was used.)
arm_pred_FE <- predict(arm_fit_FE, baseline = distr(qnorm, mean = -0.73, sd = 21^-0.5), type = "response", baseline_trt = 1) arm_pred_FE plot(arm_pred_FE)
arm_pred_RE <- predict(arm_fit_RE, baseline = distr(qnorm, mean = -0.73, sd = 21^-0.5), type = "response", baseline_trt = 1) arm_pred_RE plot(arm_pred_RE)
If the baseline argument is omitted, predictions of mean off-time reduction will be produced for every study in the network based on their estimated baseline response $\mu_j$:
arm_pred_FE_studies <- predict(arm_fit_FE, type = "response") arm_pred_FE_studies plot(arm_pred_FE_studies)
We can also produce treatment rankings, rank probabilities, and cumulative rank probabilities.
(arm_ranks <- posterior_ranks(arm_fit_FE)) plot(arm_ranks)
(arm_rankprobs <- posterior_rank_probs(arm_fit_FE)) plot(arm_rankprobs)
(arm_cumrankprobs <- posterior_rank_probs(arm_fit_FE, cumulative = TRUE)) plot(arm_cumrankprobs)
Analysis of contrast-based data
We now perform an analysis using the contrast-based data (mean differences and standard errors).
Setting up the network
With contrast-level data giving the mean difference in off-time reduction (diff) and standard error (se_diff), we use the function set_agd_contrast() to set up the network.
contr_net <- set_agd_contrast(parkinsons, study = studyn, trt = trtn, y = diff, se = se_diff, sample_size = n) contr_net
The sample_size argument is optional, but enables the nodes to be weighted by sample size in the network plot.
Plot the network structure.
plot(contr_net, weight_edges = TRUE, weight_nodes = TRUE)
Meta-analysis models
We fit both fixed effect (FE) and random effects (RE) models.
Fixed effect meta-analysis
First, we fit a fixed effect model using the nma() function with trt_effects = "fixed".
We use $\mathrm{N}(0, 100^2)$ prior distributions for the treatment effects $d_k$.
We can examine the range of parameter values implied by these prior distributions with the summary() method:
summary(normal(scale = 100))
The model is fitted using the nma() function.
contr_fit_FE <- nma(contr_net, trt_effects = "fixed", prior_trt = normal(scale = 100))
contr_fit_FE <- nma(contr_net, trt_effects = "fixed", prior_trt = normal(scale = 100), iter = 10000)
Basic parameter summaries are given by the print() method:
contr_fit_FE
The prior and posterior distributions can be compared visually using the plot_prior_posterior() function:
plot_prior_posterior(contr_fit_FE)
Random effects meta-analysis
We now fit a random effects model using the nma() function with trt_effects = "random".
Again, we use $\mathrm{N}(0, 100^2)$ prior distributions for the treatment effects $d_k$, and we additionally use a $\textrm{half-N}(5^2)$ prior for the heterogeneity standard deviation $\tau$.
We can examine the range of parameter values implied by these prior distributions with the summary() method:
summary(normal(scale = 100)) summary(half_normal(scale = 5))
Fitting the RE model
contr_fit_RE <- nma(contr_net, seed = 1150676438, trt_effects = "random", prior_trt = normal(scale = 100), prior_het = half_normal(scale = 5), adapt_delta = 0.99)
contr_fit_RE <- nowarn_on_ci(nma(contr_net, seed = 1150676438, trt_effects = "random", prior_trt = normal(scale = 100), prior_het = half_normal(scale = 5), adapt_delta = 0.99, iter = 10000))
We do see a small number of divergent transition errors, which cannot simply be removed by increasing the value of the adapt_delta argument (by default set to 0.95 for RE models).
To diagnose, we use the pairs() method to investigate where in the posterior distribution these divergences are happening (indicated by red crosses):
pairs(contr_fit_RE, pars = c("d[3]", "delta[4: 4 vs. 3]", "tau"))
The divergent transitions occur in the upper tail of the heterogeneity standard deviation. In this case, with only a small number of studies, there is not very much information to estimate the heterogeneity standard deviation and the prior distribution may be too heavy-tailed. We could consider a more informative prior distribution for the heterogeneity variance to aid estimation.
Basic parameter summaries are given by the print() method:
contr_fit_RE
By default, summaries of the study-specific relative effects $\delta_{jk}$ are hidden, but could be examined by changing the pars argument:
# Not run print(contr_fit_RE, pars = c("d", "delta"))
The prior and posterior distributions can be compared visually using the plot_prior_posterior() function:
plot_prior_posterior(contr_fit_RE)
Model comparison
Model fit can be checked using the dic() function:
(contr_dic_FE <- dic(contr_fit_FE))
(contr_dic_RE <- dic(contr_fit_RE))
Both models fit the data well, having posterior mean residual deviance close to the number of data points. The DIC is similar between models, so we choose the FE model based on parsimony.
We can also examine the residual deviance contributions with the corresponding plot() method.
plot(contr_dic_FE)
plot(contr_dic_RE)
Further results
For comparison with @TSD2, we can produce relative effects against placebo using the relative_effects() function with trt_ref = 1:
(contr_releff_FE <- relative_effects(contr_fit_FE, trt_ref = 1)) plot(contr_releff_FE, ref_line = 0)
(contr_releff_RE <- relative_effects(contr_fit_RE, trt_ref = 1)) plot(contr_releff_RE, ref_line = 0)
Following @TSD2, we produce absolute predictions of the mean off-time reduction on each treatment assuming a Normal distribution for the outcomes on treatment 1 (placebo) with mean $-0.73$ and precision $21$.
We use the predict() method, where the baseline argument takes a distr() distribution object with which we specify the corresponding Normal distribution, and we specify baseline_trt = 1 to indicate that the baseline distribution corresponds to treatment 1.
(Strictly speaking, type = "response" is unnecessary here, since the identity link function was used.)
contr_pred_FE <- predict(contr_fit_FE, baseline = distr(qnorm, mean = -0.73, sd = 21^-0.5), type = "response", baseline_trt = 1) contr_pred_FE plot(contr_pred_FE)
contr_pred_RE <- predict(contr_fit_RE, baseline = distr(qnorm, mean = -0.73, sd = 21^-0.5), type = "response", baseline_trt = 1) contr_pred_RE plot(contr_pred_RE)
If the baseline argument is omitted an error will be raised, as there are no study baselines estimated in this network.
# Not run predict(contr_fit_FE, type = "response")
We can also produce treatment rankings, rank probabilities, and cumulative rank probabilities.
(contr_ranks <- posterior_ranks(contr_fit_FE)) plot(contr_ranks)
(contr_rankprobs <- posterior_rank_probs(contr_fit_FE)) plot(contr_rankprobs)
(contr_cumrankprobs <- posterior_rank_probs(contr_fit_FE, cumulative = TRUE)) plot(contr_cumrankprobs)
Analysis of mixed arm-based and contrast-based data
We now perform an analysis where some studies contribute arm-based data, and other contribute contrast-based data. Replicating @TSD2, we consider arm-based data from studies 1-3, and contrast-based data from studies 4-7.
studies <- parkinsons$studyn (parkinsons_arm <- parkinsons[studies %in% 1:3, ]) (parkinsons_contr <- parkinsons[studies %in% 4:7, ])
Setting up the network
We use the functions set_agd_arm() and set_agd_contrast() to set up the respective data sources within the network, and then combine together with combine_network().
mix_arm_net <- set_agd_arm(parkinsons_arm, study = studyn, trt = trtn, y = y, se = se, sample_size = n) mix_contr_net <- set_agd_contrast(parkinsons_contr, study = studyn, trt = trtn, y = diff, se = se_diff, sample_size = n) mix_net <- combine_network(mix_arm_net, mix_contr_net) mix_net
The sample_size argument is optional, but enables the nodes to be weighted by sample size in the network plot.
Plot the network structure.
plot(mix_net, weight_edges = TRUE, weight_nodes = TRUE)
Meta-analysis models
We fit both fixed effect (FE) and random effects (RE) models.
Fixed effect meta-analysis
First, we fit a fixed effect model using the nma() function with trt_effects = "fixed".
We use $\mathrm{N}(0, 100^2)$ prior distributions for the treatment effects $d_k$ and study-specific intercepts $\mu_j$.
We can examine the range of parameter values implied by these prior distributions with the summary() method:
summary(normal(scale = 100))
The model is fitted using the nma() function.
mix_fit_FE <- nma(mix_net, trt_effects = "fixed", prior_intercept = normal(scale = 100), prior_trt = normal(scale = 100))
mix_fit_FE <- nma(mix_net, trt_effects = "fixed", prior_intercept = normal(scale = 100), prior_trt = normal(scale = 100), iter = 10000)
Basic parameter summaries are given by the print() method:
mix_fit_FE
By default, summaries of the study-specific intercepts $\mu_j$ are hidden, but could be examined by changing the pars argument:
# Not run print(mix_fit_FE, pars = c("d", "mu"))
The prior and posterior distributions can be compared visually using the plot_prior_posterior() function:
plot_prior_posterior(mix_fit_FE)
Random effects meta-analysis
We now fit a random effects model using the nma() function with trt_effects = "random".
Again, we use $\mathrm{N}(0, 100^2)$ prior distributions for the treatment effects $d_k$ and study-specific intercepts $\mu_j$, and we additionally use a $\textrm{half-N}(5^2)$ prior for the heterogeneity standard deviation $\tau$.
We can examine the range of parameter values implied by these prior distributions with the summary() method:
summary(normal(scale = 100)) summary(half_normal(scale = 5))
Fitting the RE model
mix_fit_RE <- nma(mix_net, seed = 437219664, trt_effects = "random", prior_intercept = normal(scale = 100), prior_trt = normal(scale = 100), prior_het = half_normal(scale = 5), adapt_delta = 0.99)
mix_fit_RE <- nowarn_on_ci(nma(mix_net, seed = 437219664, trt_effects = "random", prior_intercept = normal(scale = 100), prior_trt = normal(scale = 100), prior_het = half_normal(scale = 5), adapt_delta = 0.99, iter=10000))
We do see a small number of divergent transition errors, which cannot simply be removed by increasing the value of the adapt_delta argument (by default set to 0.95 for RE models).
To diagnose, we use the pairs() method to investigate where in the posterior distribution these divergences are happening (indicated by red crosses):
pairs(mix_fit_RE, pars = c("d[3]", "delta[4: 4 vs. 3]", "tau"))
The divergent transitions occur in the upper tail of the heterogeneity standard deviation. In this case, with only a small number of studies, there is not very much information to estimate the heterogeneity standard deviation and the prior distribution may be too heavy-tailed. We could consider a more informative prior distribution for the heterogeneity variance to aid estimation.
Basic parameter summaries are given by the print() method:
mix_fit_RE
By default, summaries of the study-specific intercepts $\mu_j$ and study-specific relative effects $\delta_{jk}$ are hidden, but could be examined by changing the pars argument:
# Not run print(mix_fit_RE, pars = c("d", "mu", "delta"))
The prior and posterior distributions can be compared visually using the plot_prior_posterior() function:
plot_prior_posterior(mix_fit_RE)
Model comparison
Model fit can be checked using the dic() function:
(mix_dic_FE <- dic(mix_fit_FE))
(mix_dic_RE <- dic(mix_fit_RE))
Both models fit the data well, having posterior mean residual deviance close to the number of data points. The DIC is similar between models, so we choose the FE model based on parsimony.
We can also examine the residual deviance contributions with the corresponding plot() method.
plot(mix_dic_FE)
plot(mix_dic_RE)
Further results
For comparison with @TSD2, we can produce relative effects against placebo using the relative_effects() function with trt_ref = 1:
(mix_releff_FE <- relative_effects(mix_fit_FE, trt_ref = 1)) plot(mix_releff_FE, ref_line = 0)
(mix_releff_RE <- relative_effects(mix_fit_RE, trt_ref = 1)) plot(mix_releff_RE, ref_line = 0)
Following @TSD2, we produce absolute predictions of the mean off-time reduction on each treatment assuming a Normal distribution for the outcomes on treatment 1 (placebo) with mean $-0.73$ and precision $21$.
We use the predict() method, where the baseline argument takes a distr() distribution object with which we specify the corresponding Normal distribution, and we specify baseline_trt = 1 to indicate that the baseline distribution corresponds to treatment 1.
(Strictly speaking, type = "response" is unnecessary here, since the identity link function was used.)
mix_pred_FE <- predict(mix_fit_FE, baseline = distr(qnorm, mean = -0.73, sd = 21^-0.5), type = "response", baseline_trt = 1) mix_pred_FE plot(mix_pred_FE)
mix_pred_RE <- predict(mix_fit_RE, baseline = distr(qnorm, mean = -0.73, sd = 21^-0.5), type = "response", baseline_trt = 1) mix_pred_RE plot(mix_pred_RE)
If the baseline argument is omitted, predictions of mean off-time reduction will be produced for every arm-based study in the network based on their estimated baseline response $\mu_j$:
mix_pred_FE_studies <- predict(mix_fit_FE, type = "response") mix_pred_FE_studies plot(mix_pred_FE_studies)
We can also produce treatment rankings, rank probabilities, and cumulative rank probabilities.
(mix_ranks <- posterior_ranks(mix_fit_FE)) plot(mix_ranks)
(mix_rankprobs <- posterior_rank_probs(mix_fit_FE)) plot(mix_rankprobs)
(mix_cumrankprobs <- posterior_rank_probs(mix_fit_FE, cumulative = TRUE)) plot(mix_cumrankprobs)
References
#--- Test against TSD 2 results --- library(testthat) library(dplyr) tol <- 0.05 tol_dic <- 0.1 #--- Arm-level data --- # FE Relative effects tsd_FE <- tribble( ~trt, ~mean, ~sd , ~median, ~lower, ~upper, 2 , -1.81, 0.33, -1.81 , -2.46 ,-1.16 , 3 , -0.47, 0.49, -0.47 , -1.43 ,0.49 , 4 , -0.52, 0.48, -0.52 , -1.46 ,0.43 , 5 , -0.82, 0.52, -0.82 , -1.84 ,0.22 , ) %>% mutate(trt = ordered(trt, levels = levels(arm_net$treatments))) %>% arrange(trt) arm_releff_FE <- as.data.frame(arm_releff_FE) test_that("FE relative effects", { expect_equivalent(arm_releff_FE$mean, tsd_FE$mean, tolerance = tol) expect_equivalent(arm_releff_FE$sd, tsd_FE$sd, tolerance = tol) expect_equivalent(arm_releff_FE$`50%`, tsd_FE$median, tolerance = tol) expect_equivalent(arm_releff_FE$`2.5%`, tsd_FE$lower, tolerance = tol) expect_equivalent(arm_releff_FE$`97.5%`, tsd_FE$upper, tolerance = tol) }) # RE Relative effects tsd_RE <- tribble( ~trt, ~mean, ~sd , ~median, ~lower, ~upper, 2 , -1.85, 0.54, -1.84 , -2.91 ,-0.85 , 3 , -0.50, 0.66, -0.50 , -1.78 ,0.75 , 4 , -0.53, 0.65, -0.53 , -1.77 ,0.71 , 5 , -0.83, 0.80, -0.83 , -2.35 ,0.69 , ) %>% mutate(trt = ordered(trt, levels = levels(arm_net$treatments))) %>% arrange(trt) arm_releff_RE <- as.data.frame(arm_releff_RE) test_that("RE relative effects", { expect_equivalent(arm_releff_RE$mean, tsd_RE$mean, tolerance = tol) expect_equivalent(arm_releff_RE$`50%`, tsd_RE$median, tolerance = tol) skip_on_ci() expect_equivalent(arm_releff_RE$sd, tsd_RE$sd, tolerance = tol) expect_equivalent(arm_releff_RE$`2.5%`, tsd_RE$lower, tolerance = tol) expect_equivalent(arm_releff_RE$`97.5%`, tsd_RE$upper, tolerance = tol) }) # Heterogeneity SD arm_tau <- as.data.frame(summary(arm_fit_RE, pars = "tau")) test_that("RE heterogeneity SD", { expect_equivalent(arm_tau$mean, 0.4, tolerance = tol) expect_equivalent(arm_tau$`50%`, 0.28, tolerance = tol) skip_on_ci() expect_equivalent(arm_tau$sd, 0.43, tolerance = tol) expect_equivalent(arm_tau$`2.5%`, 0.01, tolerance = tol) expect_equivalent(arm_tau$`97.5%`, 1.55, tolerance = tol) }) # FE probabilities tsd_pred_FE <- tribble( ~trt, ~mean, ~sd , ~median, ~lower, ~upper, 1 , -0.73, 0.22, -0.73 , -1.16 ,-0.30 , 2 , -2.54, 0.40, -2.54 , -3.32 ,-1.76 , 3 , -1.21, 0.53, -1.20 , -2.25 ,-0.15 , 4 , -1.25, 0.53, -1.25 , -2.28 ,-0.21 , 5 , -1.55, 0.57, -1.55 , -2.66 ,-0.43 , ) %>% mutate(trt = ordered(trt, levels = levels(arm_net$treatments))) %>% arrange(trt) arm_pred_FE <- as.data.frame(arm_pred_FE) test_that("FE predicted probabilities at 3 years", { expect_equivalent(arm_pred_FE$mean, tsd_pred_FE$mean, tolerance = tol) expect_equivalent(arm_pred_FE$sd, tsd_pred_FE$sd, tolerance = tol) expect_equivalent(arm_pred_FE$`50%`, tsd_pred_FE$median, tolerance = tol) expect_equivalent(arm_pred_FE$`2.5%`, tsd_pred_FE$lower, tolerance = tol) expect_equivalent(arm_pred_FE$`97.5%`, tsd_pred_FE$upper, tolerance = tol) }) # RE probabilities tsd_pred_RE <- tribble( ~trt, ~mean, ~sd , ~median, ~lower, ~upper, 1 , -0.73, 0.22, -0.73 , -1.16 ,-0.30 , 2 , -2.58, 0.58, -2.57 , -3.72 ,-1.50 , 3 , -1.23, 0.70, -1.23 , -2.57 ,0.10 , 4 , -1.26, 0.69, -1.26 , -2.57 ,0.05 , 5 , -1.57, 0.83, -1.56 , -3.14 ,0.02 , ) %>% mutate(trt = ordered(trt, levels = levels(arm_net$treatments))) %>% arrange(trt) arm_pred_RE <- as.data.frame(arm_pred_RE) test_that("RE predicted probabilities at 3 years", { expect_equivalent(arm_pred_RE$mean, tsd_pred_RE$mean, tolerance = tol) expect_equivalent(arm_pred_RE$`50%`, tsd_pred_RE$median, tolerance = tol) skip_on_ci() expect_equivalent(arm_pred_RE$sd, tsd_pred_RE$sd, tolerance = tol) expect_equivalent(arm_pred_RE$`2.5%`, tsd_pred_RE$lower, tolerance = tol) expect_equivalent(arm_pred_RE$`97.5%`, tsd_pred_RE$upper, tolerance = tol) }) # FE DIC test_that("FE DIC", { expect_equivalent(arm_dic_FE$resdev, 13.3, tolerance = tol_dic) expect_equivalent(arm_dic_FE$pd, 11.0, tolerance = tol_dic) expect_equivalent(arm_dic_FE$dic, 24.3, tolerance = tol_dic) }) # RE DIC test_that("RE DIC", { expect_equivalent(arm_dic_RE$resdev, 13.6, tolerance = tol_dic) expect_equivalent(arm_dic_RE$pd, 12.4, tolerance = tol_dic) expect_equivalent(arm_dic_RE$dic, 26.0, tolerance = tol_dic) }) #--- Contrast-level data --- # Relative effects, heterogeneity, absolute effects are all equal across data scenarios contr_releff_FE <- as.data.frame(contr_releff_FE) test_that("FE relative effects", { expect_equivalent(contr_releff_FE$mean, tsd_FE$mean, tolerance = tol) expect_equivalent(contr_releff_FE$sd, tsd_FE$sd, tolerance = tol) expect_equivalent(contr_releff_FE$`50%`, tsd_FE$median, tolerance = tol) expect_equivalent(contr_releff_FE$`2.5%`, tsd_FE$lower, tolerance = tol) expect_equivalent(contr_releff_FE$`97.5%`, tsd_FE$upper, tolerance = tol) }) contr_releff_RE <- as.data.frame(contr_releff_RE) test_that("RE relative effects", { expect_equivalent(contr_releff_RE$mean, tsd_RE$mean, tolerance = tol) expect_equivalent(contr_releff_RE$`50%`, tsd_RE$median, tolerance = tol) skip_on_ci() expect_equivalent(contr_releff_RE$sd, tsd_RE$sd, tolerance = tol) expect_equivalent(contr_releff_RE$`2.5%`, tsd_RE$lower, tolerance = tol) expect_equivalent(contr_releff_RE$`97.5%`, tsd_RE$upper, tolerance = tol) }) # Heterogeneity SD contr_tau <- as.data.frame(summary(contr_fit_RE, pars = "tau")) test_that("RE heterogeneity SD", { expect_equivalent(contr_tau$mean, 0.4, tolerance = tol) expect_equivalent(contr_tau$`50%`, 0.28, tolerance = tol) skip_on_ci() expect_equivalent(contr_tau$sd, 0.43, tolerance = tol) expect_equivalent(contr_tau$`2.5%`, 0.01, tolerance = tol) expect_equivalent(contr_tau$`97.5%`, 1.55, tolerance = tol) }) contr_pred_FE <- as.data.frame(contr_pred_FE) test_that("FE predicted probabilities at 3 years", { expect_equivalent(contr_pred_FE$mean, tsd_pred_FE$mean, tolerance = tol) expect_equivalent(contr_pred_FE$sd, tsd_pred_FE$sd, tolerance = tol) expect_equivalent(contr_pred_FE$`50%`, tsd_pred_FE$median, tolerance = tol) expect_equivalent(contr_pred_FE$`2.5%`, tsd_pred_FE$lower, tolerance = tol) expect_equivalent(contr_pred_FE$`97.5%`, tsd_pred_FE$upper, tolerance = tol) }) contr_pred_RE <- as.data.frame(contr_pred_RE) test_that("RE predicted probabilities at 3 years", { expect_equivalent(contr_pred_RE$mean, tsd_pred_RE$mean, tolerance = tol) expect_equivalent(contr_pred_RE$`50%`, tsd_pred_RE$median, tolerance = tol) skip_on_ci() expect_equivalent(contr_pred_RE$sd, tsd_pred_RE$sd, tolerance = tol) expect_equivalent(contr_pred_RE$`2.5%`, tsd_pred_RE$lower, tolerance = tol) expect_equivalent(contr_pred_RE$`97.5%`, tsd_pred_RE$upper, tolerance = tol) }) test_that("Error message when using predict() on only contrast data", { expect_error(predict(contr_fit_FE), "No arm-based data \\(IPD or AgD\\) in network") }) # FE DIC test_that("FE DIC", { expect_equivalent(contr_dic_FE$resdev, 6.3, tolerance = tol_dic) expect_equivalent(contr_dic_FE$pd, 4.0, tolerance = tol_dic) expect_equivalent(contr_dic_FE$dic, 10.3, tolerance = tol_dic) }) # RE DIC test_that("RE DIC", { expect_equivalent(contr_dic_RE$resdev, 6.6, tolerance = tol_dic) expect_equivalent(contr_dic_RE$pd, 5.5, tolerance = tol_dic) expect_equivalent(contr_dic_RE$dic, 12.1, tolerance = tol_dic) }) # Test contrast-data rows in incorrect order park_reorder <- sample_frac(parkinsons, size = 1, replace = FALSE) reorder_net <- set_agd_contrast(park_reorder, study = studyn, trt = trtn, y = diff, se = se_diff, sample_size = n) reorder_fit_FE <- nma(reorder_net, seed = 1878819627, trt_effects = "fixed", prior_trt = normal(scale = 100), iter = 10000) reorder_fit_RE <- nowarn_on_ci(nma(reorder_net, seed = 1147315, trt_effects = "random", prior_trt = normal(scale = 100), prior_het = half_normal(scale = 5), adapt_delta = 0.99, iter = 10000)) expect_equivalent_nma_summary <- function(object, expected, ...) { act <- quasi_label(enquo(object), arg = "object") exp <- quasi_label(enquo(expected), arg = "expected") comp <- compare(act$val %>% as.data.frame() %>% select(-c(Bulk_ESS, Tail_ESS, Rhat)), exp$val %>% as.data.frame() %>% select(-c(Bulk_ESS, Tail_ESS, Rhat)), ...) expect(comp$equal, sprintf("%s not equal to %s.\n%s", act$lab, exp$lab, comp$message)) invisible(act$val) } test_that("Analysis of reordered contrast data gives same results", { expect_equivalent_nma_summary(relative_effects(reorder_fit_FE, trt_ref = 1), contr_releff_FE, tolerance = tol) expect_equivalent_nma_summary(predict(reorder_fit_FE, baseline = distr(qnorm, mean = -0.73, sd = 21^-0.5), type = "response", baseline_trt = 1), contr_pred_FE, tolerance = tol) expect_equivalent(dic(reorder_fit_FE)$resdev, contr_dic_FE$resdev, tolerance = tol_dic) expect_equivalent(dic(reorder_fit_FE)$pd, contr_dic_FE$pd, tolerance = tol_dic) expect_equivalent(dic(reorder_fit_FE)$dic, contr_dic_FE$dic, tolerance = tol_dic) expect_equivalent(dic(reorder_fit_RE)$resdev, contr_dic_RE$resdev, tolerance = tol_dic) expect_equivalent(dic(reorder_fit_RE)$pd, contr_dic_RE$pd, tolerance = tol_dic) expect_equivalent(dic(reorder_fit_RE)$dic, contr_dic_RE$dic, tolerance = tol_dic) skip_on_ci() expect_equivalent_nma_summary(relative_effects(reorder_fit_RE, trt_ref = 1), contr_releff_RE, tolerance = tol) expect_equivalent_nma_summary(predict(reorder_fit_RE, baseline = distr(qnorm, mean = -0.73, sd = 21^-0.5), type = "response", baseline_trt = 1), contr_pred_RE, tolerance = tol) expect_equivalent_nma_summary(summary(reorder_fit_RE, pars = "tau"), contr_tau, tolerance = tol) }) #--- Mixed data --- # Relative effects, heterogeneity, absolute effects are all equal across data scenarios mix_releff_FE <- as.data.frame(mix_releff_FE) test_that("FE relative effects", { expect_equivalent(mix_releff_FE$mean, tsd_FE$mean, tolerance = tol) expect_equivalent(mix_releff_FE$sd, tsd_FE$sd, tolerance = tol) expect_equivalent(mix_releff_FE$`50%`, tsd_FE$median, tolerance = tol) expect_equivalent(mix_releff_FE$`2.5%`, tsd_FE$lower, tolerance = tol) expect_equivalent(mix_releff_FE$`97.5%`, tsd_FE$upper, tolerance = tol) }) mix_releff_RE <- as.data.frame(mix_releff_RE) test_that("RE relative effects", { expect_equivalent(mix_releff_RE$mean, tsd_RE$mean, tolerance = tol) expect_equivalent(mix_releff_RE$`50%`, tsd_RE$median, tolerance = tol) skip_on_ci() expect_equivalent(mix_releff_RE$sd, tsd_RE$sd, tolerance = tol) expect_equivalent(mix_releff_RE$`2.5%`, tsd_RE$lower, tolerance = tol) expect_equivalent(mix_releff_RE$`97.5%`, tsd_RE$upper, tolerance = tol) }) # Heterogeneity SD mix_tau <- as.data.frame(summary(mix_fit_RE, pars = "tau")) test_that("RE heterogeneity SD", { expect_equivalent(mix_tau$mean, 0.4, tolerance = tol) expect_equivalent(mix_tau$`50%`, 0.28, tolerance = tol) skip_on_ci() expect_equivalent(mix_tau$sd, 0.43, tolerance = tol) expect_equivalent(mix_tau$`2.5%`, 0.01, tolerance = tol) expect_equivalent(mix_tau$`97.5%`, 1.55, tolerance = tol) }) mix_pred_FE <- as.data.frame(mix_pred_FE) test_that("FE predicted probabilities at 3 years", { expect_equivalent(mix_pred_FE$mean, tsd_pred_FE$mean, tolerance = tol) expect_equivalent(mix_pred_FE$sd, tsd_pred_FE$sd, tolerance = tol) expect_equivalent(mix_pred_FE$`50%`, tsd_pred_FE$median, tolerance = tol) expect_equivalent(mix_pred_FE$`2.5%`, tsd_pred_FE$lower, tolerance = tol) expect_equivalent(mix_pred_FE$`97.5%`, tsd_pred_FE$upper, tolerance = tol) }) mix_pred_RE <- as.data.frame(mix_pred_RE) test_that("RE predicted probabilities at 3 years", { expect_equivalent(mix_pred_RE$mean, tsd_pred_RE$mean, tolerance = tol) expect_equivalent(mix_pred_RE$`50%`, tsd_pred_RE$median, tolerance = tol) skip_on_ci() expect_equivalent(mix_pred_RE$sd, tsd_pred_RE$sd, tolerance = tol) expect_equivalent(mix_pred_RE$`2.5%`, tsd_pred_RE$lower, tolerance = tol) expect_equivalent(mix_pred_RE$`97.5%`, tsd_pred_RE$upper, tolerance = tol) }) # FE DIC test_that("FE DIC", { expect_equivalent(mix_dic_FE$resdev, 9.3, tolerance = tol_dic) expect_equivalent(mix_dic_FE$pd, 7.0, tolerance = tol_dic) expect_equivalent(mix_dic_FE$dic, 16.3, tolerance = tol_dic) }) # RE DIC test_that("RE DIC", { expect_equivalent(mix_dic_RE$resdev, 9.6, tolerance = tol_dic) expect_equivalent(mix_dic_RE$pd, 8.5, tolerance = tol_dic) expect_equivalent(mix_dic_RE$dic, 18.1, tolerance = tol_dic) }) # --- Check UME models --- arm_fit_FE_UME <- nma(arm_net, trt_effects = "fixed", consistency = "ume", prior_intercept = normal(scale = 100), prior_trt = normal(scale = 10), iter = 10000) # NOTE: RE UME models are not really tractable - lots of divergences - not # enough data to estimate tau mix_fit_FE_UME <- nma(mix_net, trt_effects = "fixed", consistency = "ume", prior_intercept = normal(scale = 100), prior_trt = normal(scale = 100), iter = 10000) test_that("Arm/Mixed UME results are the same", { expect_equivalent_nma_summary(summary(arm_fit_FE_UME, pars = "d"), summary(mix_fit_FE_UME, pars = "d"), tolerance = tol) }) # NOTE: Contrast (and reordered contrast) results are different, since the # choice of independent contrasts is based on the "observed" baseline arms in # the data rather than the treatment ordering contr_fit_FE_UME <- nma(contr_net, trt_effects = "fixed", consistency = "ume", prior_trt = normal(scale = 100), iter = 10000) reorder_fit_FE_UME <- nma(reorder_net, seed = 1878819627, trt_effects = "fixed", consistency = "ume", prior_trt = normal(scale = 100), iter = 10000) test_that("Contr/Reordered UME results are the same", { expect_equivalent_nma_summary(summary(contr_fit_FE_UME, pars = "d"), summary(reorder_fit_FE_UME, pars = "d"), tolerance = tol) }) # --- Test regression models with contrast data --- park_dummy <- dplyr::mutate(parkinsons, x1 = rnorm(dplyr::n(), 0, 1)) arm_net_reg <- set_agd_arm(park_dummy, study = studyn, trt = trtn, y = y, se = se, sample_size = n) arm_fit_FE_reg <- nma(arm_net_reg, trt_effects = "fixed", regression = ~x1:.trt, prior_intercept = normal(scale = 100), prior_trt = normal(scale = 10), prior_reg = normal(scale = 1), iter = 10000) contr_net_reg <- set_agd_contrast(park_dummy, study = studyn, trt = trtn, y = diff, se = se_diff, sample_size = n) contr_fit_FE_reg <- nma(contr_net_reg, trt_effects = "fixed", regression = ~x1:.trt, prior_trt = normal(scale = 100), prior_reg = normal(scale = 1), iter = 10000) mix_net_reg <- combine_network(set_agd_arm(park_dummy[studies %in% 1:3, ], study = studyn, trt = trtn, y = y, se = se, sample_size = n), set_agd_contrast(park_dummy[studies %in% 4:7, ], study = studyn, trt = trtn, y = diff, se = se_diff, sample_size = n)) mix_fit_FE_reg <- nma(mix_net_reg, trt_effects = "fixed", regression = ~x1:.trt, prior_intercept = normal(scale = 100), prior_trt = normal(scale = 100), prior_reg = normal(scale = 1), iter = 10000) reorder_net_reg <- set_agd_contrast(sample_frac(park_dummy, size = 1, replace = FALSE), study = studyn, trt = trtn, y = diff, se = se_diff, sample_size = n) reorder_fit_FE_reg <- nma(reorder_net_reg, seed = 1878819627, trt_effects = "fixed", regression = ~x1:.trt, prior_trt = normal(scale = 100), prior_reg = normal(scale = 1), iter = 10000) test_that("Regression models work with contrast data", { expect_equivalent_nma_summary(summary(arm_fit_FE_reg, pars = c("d", "beta")), summary(contr_fit_FE_reg, pars = c("d", "beta")), tolerance = tol) expect_equivalent_nma_summary(summary(arm_fit_FE_reg, pars = c("d", "beta")), summary(mix_fit_FE_reg, pars = c("d", "beta")), tolerance = tol) expect_equivalent_nma_summary(summary(contr_fit_FE_reg, pars = c("d", "beta")), summary(reorder_fit_FE_reg, pars = c("d", "beta")), tolerance = tol) }) test_that("Error message when using predict() on only contrast data", { expect_error(predict(contr_fit_FE_reg), "No arm-based data \\(IPD or AgD\\) in network") }) # ML-NMR regression models too park_dummy_mlnmr <- dplyr::mutate(parkinsons, x1_mean = rnorm(dplyr::n(), 0, 1), x1_sd = runif(dplyr::n(), 0.1, 1)) arm_net_mlnmr <- set_agd_arm(park_dummy_mlnmr, study = studyn, trt = trtn, y = y, se = se, sample_size = n) %>% add_integration(x1 = distr(qnorm, x1_mean, x1_sd), n_int = 10) arm_fit_FE_mlnmr <- nma(arm_net_mlnmr, trt_effects = "fixed", regression = ~x1:.trt, prior_intercept = normal(scale = 100), prior_trt = normal(scale = 10), prior_reg = normal(scale = 1), iter = 20000, save_warmup = FALSE) contr_net_mlnmr <- set_agd_contrast(park_dummy_mlnmr, study = studyn, trt = trtn, y = diff, se = se_diff, sample_size = n) %>% add_integration(x1 = distr(qnorm, x1_mean, x1_sd), n_int = 10) contr_fit_FE_mlnmr <- nma(contr_net_mlnmr, trt_effects = "fixed", regression = ~x1:.trt, prior_trt = normal(scale = 100), prior_reg = normal(scale = 1), iter = 20000, save_warmup = FALSE) mix_net_mlnmr <- combine_network(set_agd_arm(park_dummy_mlnmr[studies %in% 1:3, ], study = studyn, trt = trtn, y = y, se = se, sample_size = n), set_agd_contrast(park_dummy_mlnmr[studies %in% 4:7, ], study = studyn, trt = trtn, y = diff, se = se_diff, sample_size = n)) %>% add_integration(x1 = distr(qnorm, x1_mean, x1_sd), n_int = 10) mix_fit_FE_mlnmr <- nma(mix_net_mlnmr, trt_effects = "fixed", regression = ~x1:.trt, prior_intercept = normal(scale = 100), prior_trt = normal(scale = 100), prior_reg = normal(scale = 1), iter = 20000, save_warmup = FALSE) reorder_net_mlnmr <- set_agd_contrast(sample_frac(park_dummy_mlnmr, size = 1, replace = FALSE), study = studyn, trt = trtn, y = diff, se = se_diff, sample_size = n) %>% add_integration(x1 = distr(qnorm, x1_mean, x1_sd), n_int = 10) reorder_fit_FE_mlnmr <- nma(reorder_net_mlnmr, seed = 1878819627, trt_effects = "fixed", regression = ~x1:.trt, prior_trt = normal(scale = 100), prior_reg = normal(scale = 1), iter = 20000, save_warmup = FALSE) test_that("ML-NMR models work with contrast data", { expect_equivalent_nma_summary(summary(arm_fit_FE_mlnmr, pars = c("d", "beta")), summary(contr_fit_FE_mlnmr, pars = c("d", "beta")), tolerance = tol) expect_equivalent_nma_summary(summary(arm_fit_FE_mlnmr, pars = c("d", "beta")), summary(mix_fit_FE_mlnmr, pars = c("d", "beta")), tolerance = tol) expect_equivalent_nma_summary(summary(contr_fit_FE_mlnmr, pars = c("d", "beta")), summary(reorder_fit_FE_mlnmr, pars = c("d", "beta")), tolerance = tol) }) test_that("Error message when using predict() on only contrast data", { expect_error(predict(contr_fit_FE_mlnmr), "No arm-based data \\(IPD or AgD\\) in network") }) test_that("Robust to custom options(contrasts) settings", { arm_fit_FE_SAS <- withr::with_options(list(contrasts = c(ordered = "contr.SAS", unordered = "contr.SAS")), nma(arm_net, trt_effects = "fixed", prior_intercept = normal(scale = 100), prior_trt = normal(scale = 10), iter = 10000)) expect_equivalent_nma_summary(summary(arm_fit_FE_SAS), summary(arm_fit_FE), tolerance = tol) expect_equivalent_nma_summary(relative_effects(arm_fit_FE_SAS), relative_effects(arm_fit_FE), tolerance = tol) contr_fit_FE_SAS <- withr::with_options(list(contrasts = c(ordered = "contr.SAS", unordered = "contr.SAS")), nma(contr_net, trt_effects = "fixed", prior_intercept = normal(scale = 100), prior_trt = normal(scale = 10), iter = 10000)) expect_equivalent_nma_summary(summary(contr_fit_FE_SAS), summary(contr_fit_FE), tolerance = tol) expect_equivalent_nma_summary(relative_effects(contr_fit_FE_SAS), relative_effects(contr_fit_FE), tolerance = tol) mix_fit_FE_SAS <- withr::with_options(list(contrasts = c(ordered = "contr.SAS", unordered = "contr.SAS")), nma(mix_net, trt_effects = "fixed", prior_intercept = normal(scale = 100), prior_trt = normal(scale = 10), iter = 10000)) expect_equivalent_nma_summary(summary(mix_fit_FE_SAS), summary(mix_fit_FE), tolerance = tol) expect_equivalent_nma_summary(relative_effects(mix_fit_FE_SAS), relative_effects(mix_fit_FE), tolerance = tol) })
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.