bayes.nmr

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.

Models

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)$.

Implementation

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

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.

Example

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)

References

• Li, H., Chen, M. H., Ibrahim, J. G., Kim, S., Shah, A. K., Lin, J., & Tershakovec, A. M. (2019). Bayesian inference for network meta-regression using multivariate random effects with applications to cholesterol lowering drugs. Biostatistics, 20(3), 499-516.