`bayes.nmr`

is a function that fits the model introduced in *Bayesian Network Meta-Regression Models Using Heavy-Tailed Multivariate Random Effects with Covariate-Dependent Variances (submitted)*. This page serves as a self-contained tutorial for `bayes.nmr`

.

The function `bayes.nmr`

assumes a (aggregate-data and univariate) network meta-analysis model where there are more than two treatments. Trials are free to have different numbers of trials. Let $\mathcal{T}={1,\ldots,T}$ be the set of all treatments across the $K$ trials. Consider the model
\[
\overline{y}_{\cdot tk} = \boldsymbol{x}_{tk}^\prime\boldsymbol{\beta} + \tau_{tk}\gamma_{tk} + \overline{\epsilon}_{\cdot tk},
\]
where $\overline{y}_{tk}$ is the univeriate aggreate response for the $t$th treatment of $k$th trial, and $\overline{\epsilon}_{tk}\sim \mathcal{N}(0,\sigma_{tk}^2/n_{tk})$. The data include the sample standard deviations of all treatment arms, modeled by
\[
\dfrac{(n_{tk}-1)s_{tk}^2}{\sigma_{tk}^2} \sim \chi_{n_{tk}-1}^2.
\]
By the Basu’s theorem, $\overline{y}_{tk}$ and $s_{tk}^2$ are independent. $\tau_{tk}^2$ is the key component in the regression model, encoding the variance of the random effect for the $t$th treatment in the $k$th trial. Instead of treating it as unknown, each $\tau_{tk}$ is modeled as a function of a known covariate vector, i.e., $\log \tau_{tk} = \boldsymbol{z}_{tk}^\prime\boldsymbol{\phi}$, where $\boldsymbol{z}_{tk}$ is the associated aggregate covariate vector, and $\boldsymbol{\phi}$ is the corresponding coefficient vector.

Since the variance of $\gamma_{tk}$ has been absorbed into $\tau_{tk}$, $\gamma_{tk}$ is only left with correlations. This decoupling of variance and correlation facilitates fast convergence of the Markov chain Monte Carlo algorithm. The scaled random effects, denoted by $\boldsymbol{\gamma}_k = (\gamma_{1k},\ldots,\gamma_{Tk})^\prime$ follow a multivariate $t$ distribution, i.e., $\boldsymbol{\gamma}_k \sim t_T(\boldsymbol{\gamma},\boldsymbol{\rho},\nu)$, where $t_T(\boldsymbol{\mu},\Sigma,\nu)$ is a multivariate $t$ distribution with $\nu$ degrees of freedom, a location vector $\boldsymbol{\mu}$, and a scale matrix $\Sigma$. Define a noncentered reparametrized random effects $\boldsymbol{\gamma}_{k,o} = E_k^\prime(\boldsymbol{\gamma}_k - \boldsymbol{\gamma})$. Denoting $\boldsymbol{X}_k^* = (\boldsymbol{X}_k, E_k^\prime)$ and $\boldsymbol{\theta}=(\boldsymbol{\beta}^\prime,\boldsymbol{\gamma}^\prime)^\prime$, \[ \overline{\boldsymbol{y}}_k = \boldsymbol{X}_k^*\boldsymbol{\theta} + \boldsymbol{Z}_k(\boldsymbol{\phi})\boldsymbol{\gamma}_{k,o}+ \overline{\boldsymbol{\epsilon}}_k, \] where $\overline{\boldsymbol{\epsilon}}_k \sim \mathcal{N}(\boldsymbol{0},\Sigma_k)$, $\Sigma_k=\mathrm{diag}(\sigma_{kt_{k1}}^2/n_{kt_{k1}}, \ldots, \sigma_{kt_{kT_k}}^2/n_{kt_{kT_k}})$, and $\boldsymbol{Z}_k(\boldsymbol{\phi}) = \mathrm{diag}(\exp(\boldsymbol{z}_{kt_{k1}}^\prime\boldsymbol{\phi}), \ldots, \exp(\boldsymbol{z}_{kt_{kT_k}}^\prime\boldsymbol{\phi}))$. The multivariate $t$ random effects are reparametrized as a scale mixture of normals: \[ (\boldsymbol{\gamma}_{k,o}\mid \lambda_k) \overset{\text{ind}}{\sim}\mathcal{N}(\boldsymbol{0},\lambda_k^{-1}(E_k^\prime\boldsymbol{\rho}E_k)),\quad \lambda_k \overset{\text{iid}}{\sim}\mathcal{G}a\left(\frac{\nu}{2},\frac{\nu}{2} \right), \] where $\mathcal{G}a(a,b)$ indicates the gamma distribution with mean $a/b$.

The prior distributions are as follows: $\boldsymbol{\theta}\sim \mathcal{N}(\boldsymbol{0}, c_{01}\boldsymbol{I})$, $\boldsymbol{\phi}\sim \mathcal{N}(\boldsymbol{0}, c_{02}\boldsymbol{I})$, and $\pi(\sigma_{kt_{kl}}^2) \propto \sigma_{kt_{kl}}^{-2}$. For the degrees of freedom of the random-effects distribution, $\nu$, it is by default assumed to be fixed. In case the degrees of freedom are treated as random, the following hierarchical prior is used: $(\nu\mid \nu_a,\nu_b) \sim \mathcal{G}a(\nu_a,\nu_b/\nu_a)$, $\nu_a \sim \mathcal{G}a(a_4,1/b_4)$, and $\nu_b \sim \mathcal{IG}(a_5,b_5)$ where $\mathcal{IG}(a,b)$ indicates the inverse gamma distribution with mean $\mu=b/(a-1)$ and variance $\sigma^2 = \mu^2/(a-2)$.

The following table visualizes the required data format for `bayes.nmr`

:

Outcome | SD | XCovariate | ZCovariate | Trial | Treat | Npt |
---|---|---|---|---|---|---|

$\overline{y}_{\cdot 13}$ | $s_{13}$ | $\boldsymbol{x}_{13}^\prime$ | $\boldsymbol{z}_{13}^\prime$ | 1 | 3 | 1000 |

$\overline{y}_{\cdot 10}$ | $s_{10}$ | $\boldsymbol{x}_{10}^\prime$ | $\boldsymbol{z}_{10}^\prime$ | 1 | 0 | 1000 |

$\overline{y}_{\cdot 20}$ | $s_{20}$ | $\boldsymbol{x}_{20}^\prime$ | $\boldsymbol{z}_{20}^\prime$ | 2 | 0 | 1000 |

$\vdots$ | $\vdots$ | $\vdots$ | $\vdots$ | $\vdots$ | $\vdots$ | $\vdots$ |

The variables included in the table are all required except `ZCovariate`

, the covariate matrix associated with the variance of the random effects. If missing, `ZCovariate`

will default to a vector of ones. That is $\boldsymbol{Z}_k = \boldsymbol{1}$ for all $k$. This reduces to the model where $\mathrm{Var}(\bar{y}_{tk}) = \sigma_{tk}^2/n_{tk} + \exp(2\phi)$ for every $(t,k)$.

The remaining optional arguments are as follows:

`mcmc`

- a list for MCMC specification.`ndiscard`

is the number of burn-in iterations.`nskip`

configures the thinning of the MCMC. For instance, if`nskip=5`

,`bayes.nmr`

will save the posterior every 5 iterations.`nkeep`

is the size of the posterior sample. The total number of iterations will be`ndiscard + nskip * nkeep`

.`prior`

- a list of hyperparameters. The hyperparameters include`df`

,`c01`

,`c02`

,`a4`

,`b4`

,`a5`

, and`b5`

.`df`

indicates the degrees of freedom whose default value is 20. The hyperparameters`a*`

and`b*`

will take effect only if`sample_df=TRUE`

. See`control`

.`control`

- a list of tuning parameters for the Metropolis-Hastings algorithm.`lambda`

,`phi`

, and`Rho`

are sampled through the localized Metropolis algorithm.`*_stepsize`

with the asterisk replaced with one of the names above specifies the stepsize for determining the sample evaluation points in the localized Metropolis algorithm.`sample_Rho`

can be set to`FALSE`

to suppress the sampling of`Rho`

. When`sample_Rho`

is`FALSE`

, $\boldsymbol{\rho}$ will instantiated using the value given by the`init`

argument, which defaults to $\boldsymbol{\rho} = 0.5\boldsymbol{I}+0.5\boldsymbol{1}\boldsymbol{1}^\prime$ where $\boldsymbol{1}$ is the vector of ones. When`sample_df`

is`TRUE`

, $\nu$ (`df`

) will be sampled.`scale_x`

- a Boolean whether`XCovariate`

should be scaled/standardized. The effect of setting this to`TRUE`

is not limited to merely standardizing`XCovariate`

. The following generic functions will scale the posterior sample of`theta`

back to its original unit:`plot`

,`fitted`

,`summary`

, and`print`

. That is, $\theta_j \gets \theta_j/\mathrm{sd}(X^*[:,j])$ where $X^*[:,j]$ indicates the $j$th column of $\boldsymbol{X}^*$.`init`

- a list of initial values for the parameters to be sampled:`theta`

,`phi`

,`sig2`

, and`Rho`

.`Treat_order`

- a vector of unique treatments to be used for renumbering the`Treat`

vector. The first element will be assigned treatment zero, potentially indicating placebo. If not provided, the numbering will default to an alphabetical/numerical order.`Trial_order`

- a vector of unique trials. The first element will be assigned trial zero. If not provided, the numbering will default to an alphabetical/numerical order.`verbose`

- a Boolean indicating whether to print the progress bar during the MCMC sampling.

**metapack** contains a data set, `TNM`

which stands for the Triglycerides Network Meta data, consisting of 29 studies. This data set has 73 observations and 15 variables.

`Trial`

- trial identifiers from 1 to 29.`Treat`

- the treatment indicators for placebo (`PBO`

), simvastatin (`S`

), atorvastatin (`A`

), lovastatin (`L`

), rosuvastatin (`R`

), pravastatin (`P`

), ezetimibe (`E`

), simvastatin+ezetimibe (`SE`

), atorvastatin+ezetimibe (`AE`

), lovastatin+ezetimibe (`LE`

), and pravastatin+ezetimibe (`PE`

).`bldlc`

- baseline low-density lipoprotein cholesterol (LDL-C).`bhdlc`

- baseline high-density lipoprotein cholesterol (LDL-C).`btg`

- baseline triglycerides (TG)`age`

- age in years.`white`

- the proportion of white participants.`male`

- the proportion of male participants.`durat`

- the duration in weeks.`potencymed`

- the proportion of medium statin potency`potencyhigh`

- the proportion of high statin potency`ptg`

- mean percent difference in the triglycerides`sdtg`

- the standard deviation of the percent difference in the triglycerides`Npt`

- the number of participants for the arm of a trial

The construction of `ZCovariate`

requires both the domain knowledge and exploratory data analysis. We present the grouping scheme as a promising method in variance estimation. Run the following code to construct the `ZCovariate`

for the benchmark model:

```
data(TNM)
groupInfo <- list(c("PBO"), c("R"))
nz <- length(groupInfo)
ns <- nrow(TNM)
ZCovariate <- matrix(0, ns, nz)
for (j in 1:length(groupInfo)) {
for (i in 1:ns) {
if (TNM$Treat[i] %in% groupInfo[[j]]) {
ZCovariate[i, j] <- 1
}
}
}
```

The `ZCovariate`

constructed by the grouping is a binary matrix. If there are additional covariates believed to be associated with the treatment variation, they should be adjoined horizontally. The baseline triglycerides (TG) is included to accommodate the baseline effects on the variances. Since the low-density lipoprotein cholesterol (LDL-C) is known to be highly correlated with triglycerides, the baseline LDL-C is also included. Lastly, a column of ones is added to account for the common variability across treatments.

```
addz <- scale(cbind(TNM$bldlc, TNM$btg), center = TRUE, scale = TRUE)
ZCovariate <- cbind(1, ZCovariate, addz)
```

Use `set.seed`

to reproduce the output, and run the following code to fit the model.

```
XCovariate <- model.matrix(~ 0 + bldlc + bhdlc + btg + age + white +
male + bmi + potencymed + potencyhigh + durat, data = TNM)
XCovariate <- scale(XCovariate, center = TRUE, scale = FALSE)
fit <- bayes.nmr(TNM$ptg, TNM$sdtg, XCovariate, ZCovariate, TNM$Trial, TNM$Treat, TNM$Npt,
prior = list(c01 = 1.0e05, c02 = 4, df = 3),
mcmc = list(ndiscard = 2500, nskip = 1, nkeep = 10000),
Treat_order = c("PBO", "S", "A", "L", "R", "P",
"E", "SE", "AE", "LE", "PE"),
scale_x = TRUE, verbose = TRUE)
```