In dmphillippo/multinma: Bayesian Network Meta-Analysis of Individual and Aggregate Data
library(multinma)
options(mc.cores = parallel::detectCores())
library(ggplot2)
library(multinma)
nc <- switch(tolower(Sys.getenv("_R_CHECK_LIMIT_CORES_")),
"true" =, "warn" = 2,
parallel::detectCores())
options(mc.cores = nc)
library(ggplot2)
This vignette describes the analysis of 22 trials comparing beta blockers to control for preventing mortality after myocardial infarction [@Carlin1992; @TSD2].
The data are available in this package as blocker:
head(blocker)
Setting up the network
We begin by setting up the network - here just a pairwise meta-analysis.
We have arm-level count data giving the number of deaths (r) out of the total (n) in each arm, so we use the function set_agd_arm().
We set "Control" as the reference treatment.
blocker_net <- set_agd_arm(blocker,
study = studyn,
trt = trtc,
r = r,
n = n,
trt_ref = "Control")
blocker_net
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.
By default, this will use a Binomial likelihood and a logit link function, auto-detected from the data.
blocker_fit_FE <- nma(blocker_net,
trt_effects = "fixed",
prior_intercept = normal(scale = 100),
prior_trt = normal(scale = 100))
Basic parameter summaries are given by the print() method:
blocker_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(blocker_fit_FE, pars = c("d", "mu"))
The prior and posterior distributions can be compared visually using the plot_prior_posterior() function:
plot_prior_posterior(blocker_fit_FE, prior = "trt")
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
blocker_fit_RE <- nma(blocker_net,
trt_effects = "random",
prior_intercept = normal(scale = 100),
prior_trt = normal(scale = 100),
prior_het = half_normal(scale = 5))
Basic parameter summaries are given by the print() method:
blocker_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(blocker_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(blocker_fit_RE, prior = c("trt", "het"))
Model comparison
Model fit can be checked using the dic() function:
(dic_FE <- dic(blocker_fit_FE))
(dic_RE <- dic(blocker_fit_RE))
The residual deviance is lower under the RE model, which is to be expected as this model is more flexible.
However, this comes with an increased effective number of parameters (note the increase in $p_D$).
As a result, the DIC of both models is very similar and the FE model may be preferred for parsimony.
We can also examine the residual deviance contributions with the corresponding plot() method.
plot(dic_FE)
plot(dic_RE)
There are a number of points which are not very well fit by the FE model, having posterior mean residual deviance contributions greater than 1.
Study 14 is a particularly poor fit under the FE model, but its residual deviance is reduced (although still high) under the RE model.
The evidence should be given further careful examination, and consideration given to other issues such as the potential for effect-modifying covariates [@TSD2].
Leverage plots can also be produced with the plot() method, with type = "leverage".
plot(dic_FE, type = "leverage") +
# Add labels for points outside DIC=3
geom_text(aes(label = parameter), data = ~subset(., dic > 3), vjust = -0.5)
plot(dic_RE, type = "leverage") +
# Add labels for points outside DIC=3
geom_text(aes(label = parameter), data = ~subset(., dic > 3), vjust = -0.5)
These plot the leverage for each data point (i.e. its contribution to model complexity $p_D$) against the square root of the residual deviance.
The sign of the square root residual deviance is given by sign of the difference between the observed and model-predicted values, indicating whether each data point is under- or over-estimated by the model.
Contours are displayed which indicate lines of constant contribution to the DIC.
Points contributing more than 3 to the DIC are generally considered to be contributing to poor fit; here we have labelled these points using the ggplot2 function geom_text().
As with the residual deviance plots above, we again see that both arms of study 14 are fit poorly under the FE model, and their fit is improved (but still poor) under the RE model.
Further results
@TSD2 produce absolute predictions of the probability of mortality on beta blockers and control, assuming a Normal distribution on the baseline logit-probability of mortality with mean $-2.2$ and precision $3.3$.
We can replicate these results using the predict() method.
The baseline argument takes a distr() distribution object, with which we specify the corresponding Normal distribution.
We set type = "response" to produce predicted probabilities (type = "link" would produce predicted log odds).
pred_FE <- predict(blocker_fit_FE,
baseline = distr(qnorm, mean = -2.2, sd = 3.3^-0.5),
type = "response")
pred_FE
plot(pred_FE)
pred_RE <- predict(blocker_fit_RE,
baseline = distr(qnorm, mean = -2.2, sd = 3.3^-0.5),
type = "response")
pred_RE
plot(pred_RE)
If instead of information on the baseline logit-probability of mortality we have event counts, we can use these to construct a Beta distribution for the baseline probability of mortality.
For example, if 4 out of 36 individuals died on control treatment in the target population of interest, the appropriate Beta distribution for the probability would be $\textrm{Beta}(4, 36-4)$.
We can specify this Beta distribution for the baseline response using the baseline_type = "reponse" argument (the default is "link", used above for the baseline logit-probability).
pred_FE_beta <- predict(blocker_fit_FE,
baseline = distr(qbeta, 4, 36-4),
baseline_type = "response",
type = "response")
pred_FE_beta
plot(pred_FE_beta)
pred_RE_beta <- predict(blocker_fit_RE,
baseline = distr(qbeta, 4, 36-4),
baseline_type = "response",
type = "response")
pred_RE_beta
plot(pred_RE_beta)
Notice that these results are nearly equivalent to those calculated above using the Normal distribution for the baseline logit-probability, since these event counts correspond to approximately the same distribution on the logit-probability.
References
#--- Test against TSD 2 results ---
library(testthat)
library(dplyr)
tol <- 0.05
tol_dic <- 0.1
# Relative effects
blocker_FE_releff <- as.data.frame(relative_effects(blocker_fit_FE))
test_that("FE relative effects", {
expect_equivalent(blocker_FE_releff$mean, -0.26, tolerance = tol)
expect_equivalent(blocker_FE_releff$sd, 0.050, tolerance = tol)
expect_equivalent(blocker_FE_releff$`2.5%`, -0.36, tolerance = tol)
expect_equivalent(blocker_FE_releff$`50%`, -0.26, tolerance = tol)
expect_equivalent(blocker_FE_releff$`97.5%`, -0.16, tolerance = tol)
})
blocker_RE_releff <- as.data.frame(relative_effects(blocker_fit_RE))
test_that("RE relative effects", {
expect_equivalent(blocker_RE_releff$mean, -0.25, tolerance = tol)
expect_equivalent(blocker_RE_releff$sd, 0.066, tolerance = tol)
expect_equivalent(blocker_RE_releff$`2.5%`, -0.38, tolerance = tol)
expect_equivalent(blocker_RE_releff$`50%`, -0.25, tolerance = tol)
expect_equivalent(blocker_RE_releff$`97.5%`, -0.12, tolerance = tol)
})
# RE heterogeneity SD
blocker_RE_sd <- as.data.frame(summary(blocker_fit_RE, pars = "tau"))
test_that("RE heterogeneity SD", {
expect_equivalent(blocker_RE_sd$mean, 0.14, tolerance = tol)
expect_equivalent(blocker_RE_sd$sd, 0.082, tolerance = tol)
expect_equivalent(blocker_RE_sd$`2.5%`, 0.01, tolerance = tol)
expect_equivalent(blocker_RE_sd$`50%`, 0.13, tolerance = tol)
expect_equivalent(blocker_RE_sd$`97.5%`, 0.32, tolerance = tol)
})
# DIC
test_that("FE DIC", {
expect_equivalent(dic_FE$resdev, 46.8, tolerance = tol_dic)
expect_equivalent(dic_FE$pd, 23.0, tolerance = tol_dic)
expect_equivalent(dic_FE$dic, 69.8, tolerance = tol_dic)
})
test_that("RE DIC", {
expect_equivalent(dic_RE$resdev, 41.9, tolerance = tol_dic)
expect_equivalent(dic_RE$pd, 28.1, tolerance = tol_dic)
expect_equivalent(dic_RE$dic, 70.0, tolerance = tol_dic)
})
# Predictions
blocker_pred_FE <- as.data.frame(pred_FE)
test_that("FE predicted probabilities", {
expect_equivalent(blocker_pred_FE$mean, c(0.11, 0.09), tolerance = tol_dic)
expect_equivalent(blocker_pred_FE$sd, c(0.055, 0.045), tolerance = tol_dic)
expect_equivalent(blocker_pred_FE$`2.5%`, c(0.04, 0.03), tolerance = tol_dic)
expect_equivalent(blocker_pred_FE$`50%`, c(0.10, 0.08), tolerance = tol_dic)
expect_equivalent(blocker_pred_FE$`97.5%`, c(0.25, 0.20), tolerance = tol_dic)
})
blocker_pred_RE <- as.data.frame(pred_RE)
test_that("RE predicted probabilities", {
expect_equivalent(blocker_pred_RE$mean, c(0.11, 0.09), tolerance = tol_dic)
expect_equivalent(blocker_pred_RE$sd, c(0.055, 0.046), tolerance = tol_dic)
expect_equivalent(blocker_pred_RE$`2.5%`, c(0.04, 0.03), tolerance = tol_dic)
expect_equivalent(blocker_pred_RE$`50%`, c(0.10, 0.08), tolerance = tol_dic)
expect_equivalent(blocker_pred_RE$`97.5%`, c(0.25, 0.20), tolerance = tol_dic)
})
# Check predictions with Beta distribution on baseline probability
blocker_predbeta_FE <- as.data.frame(pred_FE_beta)
test_that("FE predicted probabilities (Beta distribution)", {
expect_equal(blocker_pred_FE$mean, blocker_predbeta_FE$mean, tolerance = tol)
expect_equal(blocker_pred_FE$sd, blocker_predbeta_FE$sd, tolerance = tol)
expect_equal(blocker_pred_FE$`2.5%`, blocker_predbeta_FE$`2.5%`, tolerance = tol)
expect_equal(blocker_pred_FE$`50%`, blocker_predbeta_FE$`50%`, tolerance = tol)
expect_equal(blocker_pred_FE$`97.5%`, blocker_predbeta_FE$`97.5%`, tolerance = tol)
})
blocker_predbeta_RE <- as.data.frame(pred_RE_beta)
test_that("RE predicted probabilities (Beta distribution)", {
expect_equal(blocker_pred_RE$mean, blocker_predbeta_RE$mean, tolerance = tol)
expect_equal(blocker_pred_RE$sd, blocker_predbeta_RE$sd, tolerance = tol)
expect_equal(blocker_pred_RE$`2.5%`, blocker_predbeta_RE$`2.5%`, tolerance = tol)
expect_equal(blocker_pred_RE$`50%`, blocker_predbeta_RE$`50%`, tolerance = tol)
expect_equal(blocker_pred_RE$`97.5%`, blocker_predbeta_RE$`97.5%`, tolerance = tol)
})
# Test that ordered multinomial model is equivalent
blocker_ord_net <- set_agd_arm(blocker,
study = studyn,
trt = trtc,
r = multi(nonevents = n, events = r, inclusive = TRUE),
trt_ref = "Control")
blocker_ord_fit_FE <- nma(blocker_ord_net,
trt_effects = "fixed",
link = "logit",
prior_intercept = normal(scale = 100),
prior_trt = normal(scale = 100),
prior_aux = flat())
blocker_ord_fit_RE <- nma(blocker_ord_net,
trt_effects = "random",
link = "logit",
prior_intercept = normal(scale = 100),
prior_trt = normal(scale = 100),
prior_het = half_normal(scale = 5),
prior_aux = flat())
blocker_ord_FE_releff <- as.data.frame(relative_effects(blocker_ord_fit_FE))
test_that("Equivalent ordered multinomial FE relative effects", {
expect_equivalent(blocker_ord_FE_releff$mean, -0.26, tolerance = tol)
expect_equivalent(blocker_ord_FE_releff$sd, 0.050, tolerance = tol)
expect_equivalent(blocker_ord_FE_releff$`2.5%`, -0.36, tolerance = tol)
expect_equivalent(blocker_ord_FE_releff$`50%`, -0.26, tolerance = tol)
expect_equivalent(blocker_ord_FE_releff$`97.5%`, -0.16, tolerance = tol)
})
blocker_ord_RE_releff <- as.data.frame(relative_effects(blocker_ord_fit_RE))
test_that("Equivalent ordered multinomial RE relative effects", {
expect_equivalent(blocker_ord_RE_releff$mean, -0.25, tolerance = tol)
expect_equivalent(blocker_ord_RE_releff$sd, 0.066, tolerance = tol)
expect_equivalent(blocker_ord_RE_releff$`2.5%`, -0.38, tolerance = tol)
expect_equivalent(blocker_ord_RE_releff$`50%`, -0.25, tolerance = tol)
expect_equivalent(blocker_ord_RE_releff$`97.5%`, -0.12, tolerance = tol)
})
blocker_ord_RE_sd <- as.data.frame(summary(blocker_ord_fit_RE, pars = "tau"))
test_that("Equivalent ordered multinomial RE heterogeneity SD", {
expect_equivalent(blocker_ord_RE_sd$mean, 0.14, tolerance = tol)
expect_equivalent(blocker_ord_RE_sd$sd, 0.082, tolerance = tol)
expect_equivalent(blocker_ord_RE_sd$`2.5%`, 0.01, tolerance = tol)
expect_equivalent(blocker_ord_RE_sd$`50%`, 0.13, tolerance = tol)
expect_equivalent(blocker_ord_RE_sd$`97.5%`, 0.32, tolerance = tol)
})
test_that("Equivalent ordered multinomial FE DIC", {
expect_equivalent(dic_FE$resdev, 46.8, tolerance = tol_dic)
expect_equivalent(dic_FE$pd, 23.0, tolerance = tol_dic)
expect_equivalent(dic_FE$dic, 69.8, tolerance = tol_dic)
})
test_that("Equivalent ordered multinomial RE DIC", {
expect_equivalent(dic_RE$resdev, 41.9, tolerance = tol_dic)
expect_equivalent(dic_RE$pd, 28.1, tolerance = tol_dic)
expect_equivalent(dic_RE$dic, 70.0, tolerance = tol_dic)
})
blocker_ord_pred_FE <- as.data.frame(predict(blocker_ord_fit_FE,
baseline = distr(qnorm, mean = -2.2, sd = 3.3^-0.5),
type = "response"))
test_that("Equivalent ordered multinomial FE predicted probabilities", {
expect_equivalent(blocker_ord_pred_FE$mean, c(0.11, 0.09), tolerance = tol_dic)
expect_equivalent(blocker_ord_pred_FE$sd, c(0.055, 0.045), tolerance = tol_dic)
expect_equivalent(blocker_ord_pred_FE$`2.5%`, c(0.04, 0.03), tolerance = tol_dic)
expect_equivalent(blocker_ord_pred_FE$`50%`, c(0.10, 0.08), tolerance = tol_dic)
expect_equivalent(blocker_ord_pred_FE$`97.5%`, c(0.25, 0.20), tolerance = tol_dic)
})
blocker_ord_pred_RE <- as.data.frame(predict(blocker_ord_fit_RE,
baseline = distr(qnorm, mean = -2.2, sd = 3.3^-0.5),
type = "response"))
test_that("Equivalent ordered multinomial RE predicted probabilities", {
expect_equivalent(blocker_ord_pred_RE$mean, c(0.11, 0.09), tolerance = tol_dic)
expect_equivalent(blocker_ord_pred_RE$sd, c(0.055, 0.046), tolerance = tol_dic)
expect_equivalent(blocker_ord_pred_RE$`2.5%`, c(0.04, 0.03), tolerance = tol_dic)
expect_equivalent(blocker_ord_pred_RE$`50%`, c(0.10, 0.08), tolerance = tol_dic)
expect_equivalent(blocker_ord_pred_RE$`97.5%`, c(0.25, 0.20), tolerance = tol_dic)
})
# Check predictions with Beta distribution on baseline probability
blocker_ord_predbeta_FE <- predict(blocker_ord_fit_FE,
baseline = distr(qbeta, 4, 36 - 4), #3.66565, 36.74819 - 3.66565),
baseline_type = "response",
type = "response") %>%
as.data.frame()
test_that("FE ordered multinomial predicted probabilities (Beta distribution)", {
expect_equal(blocker_ord_pred_FE$mean, blocker_ord_predbeta_FE$mean, tolerance = tol)
expect_equal(blocker_ord_pred_FE$sd, blocker_ord_predbeta_FE$sd, tolerance = tol)
expect_equal(blocker_ord_pred_FE$`2.5%`, blocker_ord_predbeta_FE$`2.5%`, tolerance = tol)
expect_equal(blocker_ord_pred_FE$`50%`, blocker_ord_predbeta_FE$`50%`, tolerance = tol)
expect_equal(blocker_ord_pred_FE$`97.5%`, blocker_ord_predbeta_FE$`97.5%`, tolerance = tol)
})
blocker_ord_predbeta_RE <- predict(blocker_ord_fit_RE,
baseline = distr(qbeta, 4, 36 - 4), #3.66565, 36.74819 - 3.66565),
baseline_type = "response",
type = "response") %>%
as.data.frame()
test_that("RE ordered multinomial predicted probabilities (Beta distribution)", {
expect_equal(blocker_ord_pred_RE$mean, blocker_ord_predbeta_RE$mean, tolerance = tol)
expect_equal(blocker_ord_pred_RE$sd, blocker_ord_predbeta_RE$sd, tolerance = tol)
expect_equal(blocker_ord_pred_RE$`2.5%`, blocker_ord_predbeta_RE$`2.5%`, tolerance = tol)
expect_equal(blocker_ord_pred_RE$`50%`, blocker_ord_predbeta_RE$`50%`, tolerance = tol)
expect_equal(blocker_ord_pred_RE$`97.5%`, blocker_ord_predbeta_RE$`97.5%`, 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.
library(multinma) options(mc.cores = parallel::detectCores()) library(ggplot2)
library(multinma) nc <- switch(tolower(Sys.getenv("_R_CHECK_LIMIT_CORES_")), "true" =, "warn" = 2, parallel::detectCores()) options(mc.cores = nc) library(ggplot2)
This vignette describes the analysis of 22 trials comparing beta blockers to control for preventing mortality after myocardial infarction [@Carlin1992; @TSD2].
The data are available in this package as blocker:
head(blocker)
Setting up the network
We begin by setting up the network - here just a pairwise meta-analysis.
We have arm-level count data giving the number of deaths (r) out of the total (n) in each arm, so we use the function set_agd_arm().
We set "Control" as the reference treatment.
blocker_net <- set_agd_arm(blocker, study = studyn, trt = trtc, r = r, n = n, trt_ref = "Control") blocker_net
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.
By default, this will use a Binomial likelihood and a logit link function, auto-detected from the data.
blocker_fit_FE <- nma(blocker_net, trt_effects = "fixed", prior_intercept = normal(scale = 100), prior_trt = normal(scale = 100))
Basic parameter summaries are given by the print() method:
blocker_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(blocker_fit_FE, pars = c("d", "mu"))
The prior and posterior distributions can be compared visually using the plot_prior_posterior() function:
plot_prior_posterior(blocker_fit_FE, prior = "trt")
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
blocker_fit_RE <- nma(blocker_net, trt_effects = "random", prior_intercept = normal(scale = 100), prior_trt = normal(scale = 100), prior_het = half_normal(scale = 5))
Basic parameter summaries are given by the print() method:
blocker_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(blocker_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(blocker_fit_RE, prior = c("trt", "het"))
Model comparison
Model fit can be checked using the dic() function:
(dic_FE <- dic(blocker_fit_FE))
(dic_RE <- dic(blocker_fit_RE))
The residual deviance is lower under the RE model, which is to be expected as this model is more flexible. However, this comes with an increased effective number of parameters (note the increase in $p_D$). As a result, the DIC of both models is very similar and the FE model may be preferred for parsimony.
We can also examine the residual deviance contributions with the corresponding plot() method.
plot(dic_FE)
plot(dic_RE)
There are a number of points which are not very well fit by the FE model, having posterior mean residual deviance contributions greater than 1. Study 14 is a particularly poor fit under the FE model, but its residual deviance is reduced (although still high) under the RE model. The evidence should be given further careful examination, and consideration given to other issues such as the potential for effect-modifying covariates [@TSD2].
Leverage plots can also be produced with the plot() method, with type = "leverage".
plot(dic_FE, type = "leverage") + # Add labels for points outside DIC=3 geom_text(aes(label = parameter), data = ~subset(., dic > 3), vjust = -0.5)
plot(dic_RE, type = "leverage") + # Add labels for points outside DIC=3 geom_text(aes(label = parameter), data = ~subset(., dic > 3), vjust = -0.5)
These plot the leverage for each data point (i.e. its contribution to model complexity $p_D$) against the square root of the residual deviance.
The sign of the square root residual deviance is given by sign of the difference between the observed and model-predicted values, indicating whether each data point is under- or over-estimated by the model.
Contours are displayed which indicate lines of constant contribution to the DIC.
Points contributing more than 3 to the DIC are generally considered to be contributing to poor fit; here we have labelled these points using the ggplot2 function geom_text().
As with the residual deviance plots above, we again see that both arms of study 14 are fit poorly under the FE model, and their fit is improved (but still poor) under the RE model.
Further results
@TSD2 produce absolute predictions of the probability of mortality on beta blockers and control, assuming a Normal distribution on the baseline logit-probability of mortality with mean $-2.2$ and precision $3.3$.
We can replicate these results using the predict() method.
The baseline argument takes a distr() distribution object, with which we specify the corresponding Normal distribution.
We set type = "response" to produce predicted probabilities (type = "link" would produce predicted log odds).
pred_FE <- predict(blocker_fit_FE, baseline = distr(qnorm, mean = -2.2, sd = 3.3^-0.5), type = "response") pred_FE plot(pred_FE)
pred_RE <- predict(blocker_fit_RE, baseline = distr(qnorm, mean = -2.2, sd = 3.3^-0.5), type = "response") pred_RE plot(pred_RE)
If instead of information on the baseline logit-probability of mortality we have event counts, we can use these to construct a Beta distribution for the baseline probability of mortality.
For example, if 4 out of 36 individuals died on control treatment in the target population of interest, the appropriate Beta distribution for the probability would be $\textrm{Beta}(4, 36-4)$.
We can specify this Beta distribution for the baseline response using the baseline_type = "reponse" argument (the default is "link", used above for the baseline logit-probability).
pred_FE_beta <- predict(blocker_fit_FE, baseline = distr(qbeta, 4, 36-4), baseline_type = "response", type = "response") pred_FE_beta plot(pred_FE_beta)
pred_RE_beta <- predict(blocker_fit_RE, baseline = distr(qbeta, 4, 36-4), baseline_type = "response", type = "response") pred_RE_beta plot(pred_RE_beta)
Notice that these results are nearly equivalent to those calculated above using the Normal distribution for the baseline logit-probability, since these event counts correspond to approximately the same distribution on the logit-probability.
References
#--- Test against TSD 2 results --- library(testthat) library(dplyr) tol <- 0.05 tol_dic <- 0.1 # Relative effects blocker_FE_releff <- as.data.frame(relative_effects(blocker_fit_FE)) test_that("FE relative effects", { expect_equivalent(blocker_FE_releff$mean, -0.26, tolerance = tol) expect_equivalent(blocker_FE_releff$sd, 0.050, tolerance = tol) expect_equivalent(blocker_FE_releff$`2.5%`, -0.36, tolerance = tol) expect_equivalent(blocker_FE_releff$`50%`, -0.26, tolerance = tol) expect_equivalent(blocker_FE_releff$`97.5%`, -0.16, tolerance = tol) }) blocker_RE_releff <- as.data.frame(relative_effects(blocker_fit_RE)) test_that("RE relative effects", { expect_equivalent(blocker_RE_releff$mean, -0.25, tolerance = tol) expect_equivalent(blocker_RE_releff$sd, 0.066, tolerance = tol) expect_equivalent(blocker_RE_releff$`2.5%`, -0.38, tolerance = tol) expect_equivalent(blocker_RE_releff$`50%`, -0.25, tolerance = tol) expect_equivalent(blocker_RE_releff$`97.5%`, -0.12, tolerance = tol) }) # RE heterogeneity SD blocker_RE_sd <- as.data.frame(summary(blocker_fit_RE, pars = "tau")) test_that("RE heterogeneity SD", { expect_equivalent(blocker_RE_sd$mean, 0.14, tolerance = tol) expect_equivalent(blocker_RE_sd$sd, 0.082, tolerance = tol) expect_equivalent(blocker_RE_sd$`2.5%`, 0.01, tolerance = tol) expect_equivalent(blocker_RE_sd$`50%`, 0.13, tolerance = tol) expect_equivalent(blocker_RE_sd$`97.5%`, 0.32, tolerance = tol) }) # DIC test_that("FE DIC", { expect_equivalent(dic_FE$resdev, 46.8, tolerance = tol_dic) expect_equivalent(dic_FE$pd, 23.0, tolerance = tol_dic) expect_equivalent(dic_FE$dic, 69.8, tolerance = tol_dic) }) test_that("RE DIC", { expect_equivalent(dic_RE$resdev, 41.9, tolerance = tol_dic) expect_equivalent(dic_RE$pd, 28.1, tolerance = tol_dic) expect_equivalent(dic_RE$dic, 70.0, tolerance = tol_dic) }) # Predictions blocker_pred_FE <- as.data.frame(pred_FE) test_that("FE predicted probabilities", { expect_equivalent(blocker_pred_FE$mean, c(0.11, 0.09), tolerance = tol_dic) expect_equivalent(blocker_pred_FE$sd, c(0.055, 0.045), tolerance = tol_dic) expect_equivalent(blocker_pred_FE$`2.5%`, c(0.04, 0.03), tolerance = tol_dic) expect_equivalent(blocker_pred_FE$`50%`, c(0.10, 0.08), tolerance = tol_dic) expect_equivalent(blocker_pred_FE$`97.5%`, c(0.25, 0.20), tolerance = tol_dic) }) blocker_pred_RE <- as.data.frame(pred_RE) test_that("RE predicted probabilities", { expect_equivalent(blocker_pred_RE$mean, c(0.11, 0.09), tolerance = tol_dic) expect_equivalent(blocker_pred_RE$sd, c(0.055, 0.046), tolerance = tol_dic) expect_equivalent(blocker_pred_RE$`2.5%`, c(0.04, 0.03), tolerance = tol_dic) expect_equivalent(blocker_pred_RE$`50%`, c(0.10, 0.08), tolerance = tol_dic) expect_equivalent(blocker_pred_RE$`97.5%`, c(0.25, 0.20), tolerance = tol_dic) }) # Check predictions with Beta distribution on baseline probability blocker_predbeta_FE <- as.data.frame(pred_FE_beta) test_that("FE predicted probabilities (Beta distribution)", { expect_equal(blocker_pred_FE$mean, blocker_predbeta_FE$mean, tolerance = tol) expect_equal(blocker_pred_FE$sd, blocker_predbeta_FE$sd, tolerance = tol) expect_equal(blocker_pred_FE$`2.5%`, blocker_predbeta_FE$`2.5%`, tolerance = tol) expect_equal(blocker_pred_FE$`50%`, blocker_predbeta_FE$`50%`, tolerance = tol) expect_equal(blocker_pred_FE$`97.5%`, blocker_predbeta_FE$`97.5%`, tolerance = tol) }) blocker_predbeta_RE <- as.data.frame(pred_RE_beta) test_that("RE predicted probabilities (Beta distribution)", { expect_equal(blocker_pred_RE$mean, blocker_predbeta_RE$mean, tolerance = tol) expect_equal(blocker_pred_RE$sd, blocker_predbeta_RE$sd, tolerance = tol) expect_equal(blocker_pred_RE$`2.5%`, blocker_predbeta_RE$`2.5%`, tolerance = tol) expect_equal(blocker_pred_RE$`50%`, blocker_predbeta_RE$`50%`, tolerance = tol) expect_equal(blocker_pred_RE$`97.5%`, blocker_predbeta_RE$`97.5%`, tolerance = tol) }) # Test that ordered multinomial model is equivalent blocker_ord_net <- set_agd_arm(blocker, study = studyn, trt = trtc, r = multi(nonevents = n, events = r, inclusive = TRUE), trt_ref = "Control") blocker_ord_fit_FE <- nma(blocker_ord_net, trt_effects = "fixed", link = "logit", prior_intercept = normal(scale = 100), prior_trt = normal(scale = 100), prior_aux = flat()) blocker_ord_fit_RE <- nma(blocker_ord_net, trt_effects = "random", link = "logit", prior_intercept = normal(scale = 100), prior_trt = normal(scale = 100), prior_het = half_normal(scale = 5), prior_aux = flat()) blocker_ord_FE_releff <- as.data.frame(relative_effects(blocker_ord_fit_FE)) test_that("Equivalent ordered multinomial FE relative effects", { expect_equivalent(blocker_ord_FE_releff$mean, -0.26, tolerance = tol) expect_equivalent(blocker_ord_FE_releff$sd, 0.050, tolerance = tol) expect_equivalent(blocker_ord_FE_releff$`2.5%`, -0.36, tolerance = tol) expect_equivalent(blocker_ord_FE_releff$`50%`, -0.26, tolerance = tol) expect_equivalent(blocker_ord_FE_releff$`97.5%`, -0.16, tolerance = tol) }) blocker_ord_RE_releff <- as.data.frame(relative_effects(blocker_ord_fit_RE)) test_that("Equivalent ordered multinomial RE relative effects", { expect_equivalent(blocker_ord_RE_releff$mean, -0.25, tolerance = tol) expect_equivalent(blocker_ord_RE_releff$sd, 0.066, tolerance = tol) expect_equivalent(blocker_ord_RE_releff$`2.5%`, -0.38, tolerance = tol) expect_equivalent(blocker_ord_RE_releff$`50%`, -0.25, tolerance = tol) expect_equivalent(blocker_ord_RE_releff$`97.5%`, -0.12, tolerance = tol) }) blocker_ord_RE_sd <- as.data.frame(summary(blocker_ord_fit_RE, pars = "tau")) test_that("Equivalent ordered multinomial RE heterogeneity SD", { expect_equivalent(blocker_ord_RE_sd$mean, 0.14, tolerance = tol) expect_equivalent(blocker_ord_RE_sd$sd, 0.082, tolerance = tol) expect_equivalent(blocker_ord_RE_sd$`2.5%`, 0.01, tolerance = tol) expect_equivalent(blocker_ord_RE_sd$`50%`, 0.13, tolerance = tol) expect_equivalent(blocker_ord_RE_sd$`97.5%`, 0.32, tolerance = tol) }) test_that("Equivalent ordered multinomial FE DIC", { expect_equivalent(dic_FE$resdev, 46.8, tolerance = tol_dic) expect_equivalent(dic_FE$pd, 23.0, tolerance = tol_dic) expect_equivalent(dic_FE$dic, 69.8, tolerance = tol_dic) }) test_that("Equivalent ordered multinomial RE DIC", { expect_equivalent(dic_RE$resdev, 41.9, tolerance = tol_dic) expect_equivalent(dic_RE$pd, 28.1, tolerance = tol_dic) expect_equivalent(dic_RE$dic, 70.0, tolerance = tol_dic) }) blocker_ord_pred_FE <- as.data.frame(predict(blocker_ord_fit_FE, baseline = distr(qnorm, mean = -2.2, sd = 3.3^-0.5), type = "response")) test_that("Equivalent ordered multinomial FE predicted probabilities", { expect_equivalent(blocker_ord_pred_FE$mean, c(0.11, 0.09), tolerance = tol_dic) expect_equivalent(blocker_ord_pred_FE$sd, c(0.055, 0.045), tolerance = tol_dic) expect_equivalent(blocker_ord_pred_FE$`2.5%`, c(0.04, 0.03), tolerance = tol_dic) expect_equivalent(blocker_ord_pred_FE$`50%`, c(0.10, 0.08), tolerance = tol_dic) expect_equivalent(blocker_ord_pred_FE$`97.5%`, c(0.25, 0.20), tolerance = tol_dic) }) blocker_ord_pred_RE <- as.data.frame(predict(blocker_ord_fit_RE, baseline = distr(qnorm, mean = -2.2, sd = 3.3^-0.5), type = "response")) test_that("Equivalent ordered multinomial RE predicted probabilities", { expect_equivalent(blocker_ord_pred_RE$mean, c(0.11, 0.09), tolerance = tol_dic) expect_equivalent(blocker_ord_pred_RE$sd, c(0.055, 0.046), tolerance = tol_dic) expect_equivalent(blocker_ord_pred_RE$`2.5%`, c(0.04, 0.03), tolerance = tol_dic) expect_equivalent(blocker_ord_pred_RE$`50%`, c(0.10, 0.08), tolerance = tol_dic) expect_equivalent(blocker_ord_pred_RE$`97.5%`, c(0.25, 0.20), tolerance = tol_dic) }) # Check predictions with Beta distribution on baseline probability blocker_ord_predbeta_FE <- predict(blocker_ord_fit_FE, baseline = distr(qbeta, 4, 36 - 4), #3.66565, 36.74819 - 3.66565), baseline_type = "response", type = "response") %>% as.data.frame() test_that("FE ordered multinomial predicted probabilities (Beta distribution)", { expect_equal(blocker_ord_pred_FE$mean, blocker_ord_predbeta_FE$mean, tolerance = tol) expect_equal(blocker_ord_pred_FE$sd, blocker_ord_predbeta_FE$sd, tolerance = tol) expect_equal(blocker_ord_pred_FE$`2.5%`, blocker_ord_predbeta_FE$`2.5%`, tolerance = tol) expect_equal(blocker_ord_pred_FE$`50%`, blocker_ord_predbeta_FE$`50%`, tolerance = tol) expect_equal(blocker_ord_pred_FE$`97.5%`, blocker_ord_predbeta_FE$`97.5%`, tolerance = tol) }) blocker_ord_predbeta_RE <- predict(blocker_ord_fit_RE, baseline = distr(qbeta, 4, 36 - 4), #3.66565, 36.74819 - 3.66565), baseline_type = "response", type = "response") %>% as.data.frame() test_that("RE ordered multinomial predicted probabilities (Beta distribution)", { expect_equal(blocker_ord_pred_RE$mean, blocker_ord_predbeta_RE$mean, tolerance = tol) expect_equal(blocker_ord_pred_RE$sd, blocker_ord_predbeta_RE$sd, tolerance = tol) expect_equal(blocker_ord_pred_RE$`2.5%`, blocker_ord_predbeta_RE$`2.5%`, tolerance = tol) expect_equal(blocker_ord_pred_RE$`50%`, blocker_ord_predbeta_RE$`50%`, tolerance = tol) expect_equal(blocker_ord_pred_RE$`97.5%`, blocker_ord_predbeta_RE$`97.5%`, 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.