2.2.1 Hierarchical model

To allow and capture heterogeneity between studies in prevalence, sensitivities, specificities, and correlations among them, we include multivariate random effects in the meta-regression model and writeReference Ma, Lian, Chu, Ibrahim and Chen ⁸

(2.3) $$ \begin{align} \pi_{i} = \Phi(\eta + \epsilon_{i}), \text{Se}_{ik} = \Phi(\alpha_{k} + \mu_{ik}), \text{Sp}_{ik} = \Phi(\beta_{k} + v_{ik}),i = 1,\ldots,N,k = 1,\ldots,K, \end{align} $$

where $\Phi (\cdot )$ is the standard normal cumulative distribution function. The parameter $\eta $ is a fixed effect for prevalence and the parameters $\alpha _{k}$ and $\beta _{k}$ are fixed effects for the sensitivity and specificity of $T_k$ , respectively. The $\epsilon _i$ is the study-specific random effect for prevalence in study i, and $\mu _{ik}$ and $\nu _{ik}$ are the study-specific random effects for the sensitivity and specificity of test $T_k$ in study i, respectively, capturing between-study heterogeneity and inducing correlation among prevalence, sensitivity, and specificity across tests.

We assume within-study random effects follow the multivariate normal distribution to account for dependency between test parameters and correlation between prevalence and test parameters if applicable, that is,

$$ \begin{align*} \boldsymbol{\theta}_i = (\epsilon_i,\mu_{i1},\nu_{i1},\ldots,\mu_{iK},\nu_{iK})^\top\sim MVN(0,\boldsymbol{\Sigma}),i = 1,\ldots,N \end{align*} $$

as in Ma et al.Reference Ma, Lian, Chu, Ibrahim and Chen ⁸ The covariance matrix $\boldsymbol {\Sigma } = \mathbf {S}\boldsymbol {\Omega }\mathbf {S}$ , where $\mathbf {S}$ is a diagonal matrix; its diagonal entries are standard deviations describing between-study heterogeneity, while $\boldsymbol {\Omega }$ has diagonal entries equal to 1 and off-diagonal entries measuring correlations between prevalence and sensitivities and specificities.

2.2.2 Prior specifications

The precision matrix $\boldsymbol {\Sigma }^{-1}$ follows a Wishart( $\mathbf {R}, v$ ) distribution, where the degrees of freedom parameter $v=2K+1$ , the dimension of $\boldsymbol {\Sigma }$ . The matrix $\mathbf {R}$ acts as a scale (inverse prior mean) for $\boldsymbol {\Sigma }^{-1}$ , so it determines the prior magnitude of between-study variability and the strength of prior dependence across random effects. Since the Wishart distribution has support only on symmetric positive-definite matrices, placing a Wishart prior on the precision matrix $\boldsymbol {\Sigma }^{-1}$ guarantees that the implied covariance matrix $\boldsymbol {\Sigma }$ is also positive definite. Consequently, all required covariance constraints, such as $a_{ij}^2 \le a_{ii} a_{jj}$ for all $i \neq j$ , are automatically satisfied. The choice of $\mathbf {R}$ can vary according to different situations; different choices of $\mathbf {R}$ give priors containing different levels of information. Section 4.2 describes different choices of priors on $\mathbf {R}$ available in NMADTA. Vague normal priors with mean 0 and variance 10 are placed on $\eta ,\alpha _k$ , and $\beta _k$ ( $k=1,\ldots ,K$ ).

2.2.3 Posterior distribution

Given the likelihood and prior specifications above, the posterior can be written asReference Ma, Lian, Chu, Ibrahim and Chen ⁸

$$ \begin{align*} \begin{aligned} &P\!\left(\eta,\boldsymbol{\alpha},\boldsymbol{\beta},\{\epsilon_i,\boldsymbol{\mu}_i,\boldsymbol{\nu}_i\}_{i=1}^N \,\middle|\, \boldsymbol{y}\right) \propto\\ & \prod_{i=1}^N \left[\prod_{j=1}^{J_i}(l_{ij1})^{\delta_{i0}} (l_{ij0})^{1-\delta_{i0}}|\boldsymbol{\Sigma}|^{-1/2} e^{-\frac{1}{2}\theta_i'\boldsymbol{\Sigma}^{-1}\theta_i} \right]\times \left\{\prod_{k=1}^K f(\alpha_k) f(\beta_k) \right\} f(\eta) f(\boldsymbol{\Sigma}), \end{aligned} \end{align*} $$

where $J_i$ is the number of subjects in study i, and $l_{ij1}$ and $l_{ij0}$ are defined above in equations (2.1) and (2.2).

Posterior medians for prevalence, sensitivity, and specificity of $T_k$ can be estimated using the MCMC sampling and these definitions: $\pi = \Phi (\eta )$ , $\text {Se}_k = \Phi (\alpha _k)$ , and $\text {Sp}_k = \Phi (\beta _k)$ . Posterior medians of other measures, such as positive and negative predictive values (PPV and NPV) and positive and negative likelihood (LR $+$ and LR $-$ ) for each candidate test, can also be calculated using

$$ \begin{align*} \begin{aligned} \text{PPV}_k = \frac{\text{Se}_k\pi}{\text{Se}_k\pi+(1-\text{Sp}_k)(1-\pi)}; &\quad \text{NPV}_k = \frac{\text{Sp}_k(1-\pi)}{\text{Sp}_k(1-\pi)+(1-\text{Se}_k)\pi};\\ \text{LR}_{+k} = \frac{\text{Se}_k}{1-\text{Sp}_k}; &\quad \text{LR}_{-k} = \frac{1-\text{Se}_k}{\text{Sp}_k}. \end{aligned} \end{align*} $$

The SROC curve is obtained by using equation (2.3) to obtain a regression line of $\alpha _k$ on $\beta _k$ Reference Arends, Hamza, Van Houwelingen, Heijenbrok-Kal, Hunink and Stijnen ³² to estimate the median probit transformed sensitivity given several values of ${1}\!{-}$ specificity:

$$ \begin{align*} \text{Se}_k = \Phi \left\{ -\alpha_k + \frac{\boldsymbol{\Sigma}_{3k-1,3k}}{\boldsymbol{\Sigma}_{3k,3k}} (\Phi^{-1}\left(1 - \text{Sp}_k) - \beta_k\right) \right\}, \end{align*} $$

where $\boldsymbol {\Sigma }$ is specified above as the covariance matrix, $\boldsymbol {\Sigma }_{3k-1,3k}$ is the correlation between $\text {Se}_k$ and $\text {Sp}_k$ , and $\boldsymbol {\Sigma }_{3k,3k}$ captures the variance of $\text {Sp}_k$ .

2.3.1 The HSROC model

The HSROC model has three layers that capture variability within a study, variation between studies, and correlations among tests in the same study.

We assume that the outcomes of all tests are dichotomous. Within a study, let the latent variable $Z_{ijk}$ follow a normal distributionReference Lian, Hodges and Chu ⁹ :

$$ \begin{align*} \begin{aligned} Z_{ijk}\sim \left\{ \begin{array}{lcl} N\left(-\alpha_{ik}/2,\exp(-\beta_k/2)\right) & \mbox{for} & y_{ij0}=0, \\ N\left(\alpha_{ik}/2,\exp(\beta_k/2)\right) & \mbox{for} & y_{ij0}=1. \end{array}\right. \end{aligned} \end{align*} $$

If $Z_{ijk}$ is greater than the study-specific cutoff value $\theta _{ik}$ , subject j in study i has a positive result for $T_k$ ; otherwise, the result is negative.

The parameter $\alpha _{ik}$ measures the difference between the distributions of $Z_{ijk}$ when disease is present versus absent in study i. The shape parameter $\beta _k$ measures the difference in the variance in the outcomes of diseased and non-diseased populations. The $\beta _k$ ’s are assumed to be the same across all studies because otherwise the model is not identified in typical datasets. Although we use the same notation $\alpha $ and $\beta $ in both model formulations to denote test-specific fixed effects, they enter the models through different parameterizations. In the SROC-based model, $\alpha _k$ and $\beta _k$ act directly on the link functions of sensitivity and specificity, whereas in the HSROC formulation, the fixed effects are reparameterized into accuracy and threshold components and enter through different linear predictors. Thus, while the underlying notion of test-specific fixed effects is shared, their interpretations depend on the chosen parameterization.

With these assumptions, test $T_k$ in study i has sensitivity and specificityReference Lian, Hodges and Chu ⁹ :

$$ \begin{align*} \text{Se}_{ik} = \Phi\left(-\frac{\theta_{ik} - \alpha_{ik}/2}{\exp(\beta_k/2)}\right), \text{Sp}_{ik} = \Phi\left(\frac{\theta_{ik} + \alpha_{ik}/2}{\exp(-\beta_k/2)}\right). \end{align*} $$

Let $\theta _{i0}$ be the cutoff for the gold standard in study i and let the true disease status $D_{ij}$ be positive if the latent variable $Z_{ij0}>\theta _{i0}$ and negative otherwise. This accounts for heterogeneity of prevalence between studies; the study-specific prevalence can be written as $\pi _i = \Phi (-\theta _{i0})$ .

For each study, let $\boldsymbol {\alpha }_i = (\alpha _{i1}, \ldots , \alpha _{iK})^\top $ and $\boldsymbol {\theta }_i = (\theta _{i0}, \theta _{i1}, \ldots , \theta _{iK})^\top $ ( $i =1,2, \ldots , N$ ) be mutually independent.Reference Rutter and Gatsonis ^{11 ,} Reference Dendukuri, Schiller, Joseph and Pai ¹⁸ Let the study-specific cutoff values $\boldsymbol {\theta }_i$ and values $\boldsymbol {\alpha }_i$ follow multivariate normal distributions to account for heterogeneity between studies and correlations between different index testsReference Lian, Hodges and Chu ⁹ :

$$ \begin{align*} \begin{aligned} \begin{pmatrix} \theta_{i0} \\ \theta_{i1} \\ \vdots \\ \theta_{iK} \end{pmatrix} &\overset{\text{iid}}{\sim} \text{MVN} \begin{pmatrix}\boldsymbol{\Theta} = \begin{pmatrix} \Theta_{0} \\ \Theta_{1} \\ \vdots \\ \Theta_{K} \end{pmatrix}, \boldsymbol{\Sigma}_{\Theta} = \begin{pmatrix} \Sigma_{0}^2 & \Sigma_{01} & \dots & \Sigma_{0K}\\ {\Sigma_{01}} & \Sigma_{1}^2 & \dots & \vdots\\ \vdots& \vdots& \ddots &\vdots \\ \Sigma_{0K} & \dots&\dots & \Sigma_{K}^2\end{pmatrix}\end{pmatrix} \ \text{and} \\ \begin{pmatrix} \alpha_{i1} \\ \alpha_{i2} \\ \vdots \\ \alpha_{iK} \end{pmatrix} &\overset{\text{iid}}{\sim} \text{MVN} \begin{pmatrix} \boldsymbol{\Lambda} = \begin{pmatrix}\Lambda_{1} \\ \Lambda_{2} \\\vdots \\\Lambda_{K} \end{pmatrix}, \boldsymbol{\Sigma}_{\Lambda} = \begin{pmatrix} \tau_{1}^2 & \tau_{12} & \dots & \tau_{1K}\\ \tau_{12} & \tau_{2}^2 & \dots & \vdots\\ \vdots& \vdots& \ddots &\vdots \\ \tau_{1K} & \dots&\dots & \tau_{K}^2 \end{pmatrix} \end{pmatrix}, \ i=1, \ldots,N, \end{aligned} \end{align*} $$

where the covariance matrices $\boldsymbol {\Sigma }_{\Theta }$ and $\boldsymbol {\Sigma }_{\Lambda }$ are positive definite with diagonal entries capturing between-study variability in the cutoff and accuracy values, and off-diagonal entries measuring covariance between tests of pairs of cutoff and accuracy values, respectively. Also, $\boldsymbol {\Sigma }_{\Theta }$ and $\boldsymbol {\Sigma }_{\Lambda }$ are independent and jointly determine the correlation between index tests’ sensitivities and specificities.Reference Lian, Hodges and Chu ⁹ This model assumes that $\boldsymbol {\Sigma }_{\Theta }$ and $\boldsymbol {\Sigma }_{\Lambda }$ are unstructured, and the data will permit inferences about these correlations.

2.3.2 Prior specifications

We assume that the mean cutoff values $\Theta _k$ and accuracy values $\Lambda _k$ , $k=1, 2, \ldots , K$ have vague normal priors with large variance. The shape parameters $\beta _k$ have a Uniform $(b_1, b_2)$ prior where $b_1$ and $b_2$ depend on prior knowledge about the index test k and $(b_1, b_2)$ should allow all possible $\beta _k$ .

We decompose the covariance matrix $\boldsymbol {\Sigma }_{\Theta }$ of the cutoff values as diag( $\mathbf {S}$ ) ${\cdot} \mathbf {R}{\cdot} $ diag( $\mathbf {S}$ ), where diag( $\mathbf {S}$ ) is a diagonal matrix with diagonal entries $\mathbf {S}= (s_0,s_1,\ldots ,s_K$ ) and $\boldsymbol {R}$ is a $(K+1) \times (K+1)$ matrix that captures the correlations between cutoff values. This decomposition, using a separation technique, can simplify specification and simplify calculations by not forcing $\mathbf {R}$ defined above to be a correlation matrix,Reference Barnard, McCulloch and Meng ³³ as follows. The cutoff values have standard deviations $\boldsymbol {\Sigma }_k = \sqrt {s_k^2 \cdot \mathbf {R}_{k,k}}$ ( $k = 0,1, \ldots ,K$ ), while the correlation matrix can be written as $\boldsymbol {\text {Corr}}_{\Theta } = \text {diag}(\mathbf {R})^{-1/2}\cdot \mathbf {R}\ \text {diag}(\mathbf {R})^{-1/2}$ , so that these items are identified even though $\mathbf {S}$ and $\mathbf {R}$ are not. Similarly, decompose the covariance matrix of the accuracy values $\boldsymbol {\Sigma }_{\Lambda }$ as $\boldsymbol {\Sigma }_{\Lambda }=\text {diag}(\mathbf {P})\cdot \boldsymbol {\Omega }\cdot \text {diag}(\mathbf {P})$ , where $\mathbf {P}=(p_1,\ldots ,p_K)$ is a diagonal matrix and the standard deviations and correlation matrix for the accuracy values are $\tau _{k} = \sqrt {p_k^2 \cdot \boldsymbol {\Omega }_{k,k}}$ and $\text {Corr}_{\Lambda } = \text {diag}(\boldsymbol {\Omega })^{-1/2}\cdot \boldsymbol {\Omega }\cdot \text {diag}(\boldsymbol {\Omega })^{-1/2}$ , respectively. Let each element of $\log (\mathbf {P})$ and $\log (\mathbf {S})$ follow an $N(\eta , \zeta ^2)$ prior, where $\zeta ^2$ is a variance, and assign values to $\eta $ and $\xi ^2$ based on prior knowledge about their standard deviations. When prior information does not suggest different magnitudes of between-study variability for cutoff and accuracy components, we can use the same weakly informative hyperparameters. Let $\mathbf {R}$ and $\boldsymbol {\Omega }$ follow Inverse-Wishart priors $\text {IW}(\mathbf {I}_{K+1},K+2)$ and $IW(\mathbf {I}_K,K+1),$ respectively, where $\mathbf {I}_{K+1}$ and $\mathbf {I}_K$ are identity matrices of dimensions $K+1$ and $K,$ respectively. These prior specifications imply that the marginal priors of all correlation parameters are approximately uniform. Studies including all tests contribute to the full matrices ${\boldsymbol {\Sigma }_{\Theta }}$ and $\boldsymbol {\Sigma }_{\Lambda }$ while studies with missing tests contribute only to submatrices of them. We note that identifiability can be challenging in Bayesian NMA models for diagnostic tests, particularly when the network is sparse (e.g., many studies include only a subset of tests) or when the reference gold standard is frequently missing. In such settings, posterior inferences for variance–covariance components and correlation structures may be more sensitive to prior choices, and MCMC mixing can be slower. We therefore recommend routine convergence diagnostics and prior sensitivity analyses in case studies.

2.3.3 Posterior distribution

Using the priors and likelihood functions specified above, the joint posterior distribution can be written asReference Lian, Hodges and Chu ⁹

$$ \begin{align*} \begin{aligned} & L(\boldsymbol{\theta}_{i}, \boldsymbol{\alpha}_{i},\boldsymbol{\beta},\boldsymbol{\Theta},\boldsymbol{\Sigma}_{\Theta},\boldsymbol{\Lambda},\boldsymbol{\Sigma}_{\Lambda}|\boldsymbol{y}^o)\\ & \propto \prod_{i} \left[\prod_{j} \Big( (l_{ij1})^{\delta_{i0}} (l_{ij0})^{1-\delta_{i0}}\Big) \times |\boldsymbol{\Sigma}_{\Theta}|^{-1/2} e^{-\frac{1}{2}(\boldsymbol{\theta}_i-\boldsymbol{\Theta})^{'} \boldsymbol{\Sigma}_{\Theta}^{-1}(\boldsymbol{\theta}_i-\boldsymbol{\Theta})} |\boldsymbol{\Sigma}_{\Lambda}|^{-1/2} e^{-\frac{1}{2}(\boldsymbol{\alpha}_i-\boldsymbol{\Lambda})^{'} {\boldsymbol{\Sigma}_{\Lambda}}^{-1}(\boldsymbol{\alpha}_i-\boldsymbol{\Lambda})}\right]\nonumber\\ & \times f_{\boldsymbol{\Theta}}(\boldsymbol{\Theta}) f_{\boldsymbol{\Lambda}}(\boldsymbol{\Lambda}) f_{\mathbf{P}}(\mathbf{P}) f_{\mathbf{S}}(\mathbf{S}) f_{\mathbf{R}}(\mathbf{R}) f_{\boldsymbol{\Omega}}(\boldsymbol{\Omega})f_{\boldsymbol{\beta}}(\boldsymbol{\beta}). \end{aligned} \end{align*} $$

The total number of degrees of freedom of this model varies with both the number of index tests included in a study and whether a gold standard test is included. When results from a gold standard test are available, a study contributes $2^{(K_i+1)}-1$ degrees of freedom; otherwise, it contributes $2^{K_i}-1$ degrees of freedom, where $K_i$ is the number of candidate tests included in study i. Thus, the total degrees of freedom could vary from $3N$ to $N\times (2^{K+1}-1)$ depending on the designs of the N studies. Since the total number of parameters is at least $K^2+5K+2$ , the minimum number of studies N required to estimate this model without setting informative priors is between $ \frac {K^2+5K+2}{2^{K+1}-1} $ and $ \frac {K^2+5K+2}{3} $ depending on the study designs.Reference Lian, Hodges and Chu ⁹

We can sample from the joint posterior distribution using MCMC sampling algorithm implemented in JAGS.Reference Plummer, Stukalov and Denwood ³⁴ The overall prevalence of the disease ( $\Pi $ ) and the sensitivity and specificity of each test can be written asReference Lian, Hodges and Chu ⁹

$$ \begin{align*} \begin{aligned} \Pi = \Phi\left(-\Theta_0\right), \text{Se}_{k} = \Phi\left(-\dfrac{\Theta_{k}-\Lambda_{k}/2}{\exp(\beta_k/2)}\right), \text{Sp}_{k} = \Phi\left(\dfrac{\Theta_{k} + \Lambda_{k}/2}{\exp(-\beta_k/2)}\right). \end{aligned} \end{align*} $$

We will use medians and equal-tailed credible intervals of these posterior samples for inferences.

Other accuracy measures, such as a diagnostic test’s PPV and NPV, can also be calculated from the MCMC samples.Reference Lian, Hodges and Chu ⁹ The SROC curve can be derived using posterior samples of $\Lambda $ , $\Theta ,$ and $\beta $ . As in Rutter and Gatsonis,Reference Rutter and Gatsonis ¹¹ we can use equation (2.4) to express the test-specific sensitivities in terms of given values of specificities by writing $\Theta _k$ as an expression of $\Lambda _k$ , $\beta _k,$ and $Sp_k$ and eliminating $\Theta _k$ in the expression of $\text {Se}_k$ Reference Lian, Hodges and Chu ⁹ :

$$ \begin{align*} \text{Se}_k = \Phi \left\{-\dfrac{\Phi^{-1}(\text{Sp}_k)\exp(-\beta_k/2)-\Lambda_k}{\exp(\beta_k/2)}\right\}. \end{align*} $$

4.2.1 Function nmadt.hierarchical() for the regular Bayesian hierarchical model

The arguments of the function nmadt.hierarchical are as follows:

Long description

The code block consists of five lines.

Line 1: nmadt dot hierarchical open parenthesis n s t u, K, data, testname, directory equals N U L L,

Line 2: plus diag equals 5, off underscore diag equals 0.05, digits equals 4, mu underscore alpha equals 0, mu underscore beta equals 0,

Line 3: plus mu underscore eta equals minus 0, preci underscore alpha equals 0.1, preci underscore beta equals 0.1, preci underscore eta equals 0.1,

Line 4: plus n dot adapt equals 5000, n dot iter equals 50000, n dot chains equals 3,

Line 5: plus conv dot diag equals F A L S E, trace equals N U L L, dic equals F A L S E, m c m c dot samples equals F A L S E close parenthesis.

Users need to input a dataset with structure similar to dat.kang, described above. The arguments nstu and K specify the number of studies and the number of candidate diagnostic tests included in the dataset. The argument data is the dataset to be meta-analyzed. The argument testname is a string vector that specifies the names of the candidate diagnostic tests. The argument directory is a string that specifies the directory to save trace plots or coda samples when trace or mcmc.samples are set to TRUE. The arguments diag and off_diag specify the values of diagonal and off-diagonal entries of the scale matrix $\mathbf {R}$ of the precision matrix $\Sigma $ . The arguments diag and off_diag determine the scale matrix $\mathbf {R}$ in the Wishart prior for the precision matrix, $\boldsymbol {\Sigma }^{-1}\sim \text {Wishart}(\mathbf {R},\nu )$ . In NMADTA, $\mathbf {R}$ is constructed with diagonal entries equal to diag and off-diagonal entries equal to off_diag. Larger values of diag typically imply a strong prior shrinkage, favoring smaller between-study variability in prevalence, sensitivity, and specificity, whereas off_diag controls the strength of prior dependence among these random effects (e.g., correlations between sensitivity and specificity and their association with prevalence). The default setting (diag=5, off_diag=0.05) is weakly informative and works well in most applications; if strong correlations are expected, users may increase off_diag, and if substantially larger heterogeneity is anticipated, users may decrease diag. We recommend assessing sensitivity of posterior summaries to these choices in conjunction with standard MCMC convergence diagnostics. The argument digits is used to specify the number of digits to the right of the decimal point to keep for the results; the default is digits = 4.

We assume that the priors for $\boldsymbol {\alpha }$ , $\boldsymbol {\beta }$ , and $\eta $ follow a vague normal distribution with mean $0$ and variance $10$ (the default), which corresponds to $95\%$ prior credible intervals of approximately $(0,1)$ for test specific prevalence, sensitivity, and specificity.Reference Ma, Lian, Chu, Ibrahim and Chen ⁸

The arguments n.adapt, n.iter, n.burnin, n.chains, and n.thin specify the MCMC algorithm run by rjags. The argument n.adapt is the number of iterations for adaptation (the default is 5,000). The argument n.iter is the number of iterations in each MCMC chain, and n.burnin is the number of burn-in iterations at the beginning of each chain that are not saved. The argument n.chains=3 (the default) is the number of MCMC chains. The argument n.thin is the thinning rate for each chain, which can be used to save memory and computation time when n.iter is large. For example, the algorithm saves only one sample in every $n^{\text {th}}$ iteration, where n is given by n.thin.

The argument conv.diag specifies whether to compute potential scale reduction factors (PSRFs) proposed by Gelman and RubinReference Gelman and Rubin ³⁸ for convergence diagnostics. The argument trace is a string vector describing the subset of quantities for which trace plots can be drawn, which can include prevalence ("prev"), sensitivity ("Se"), specificity ("Sp"), positive and negative predictive values ("ppv" and "npv", respectively), positive likelihood ("LRpos"), and negative likelihood ("LRneg"). Trace plots for the specified quantities are saved in the user-specified directory as .png files. Both PSRF statistics and trace plots can be used to judge whether the MCMC chains converged eventually. Finally, if dic = TRUE, the function will output the deviance information criterion (DIC) statistic proposed by Spiegelhalter et al.Reference Spiegelhalter, Best, Carlin and Van Der Linde ³⁹ The argument mcmc.samples is a logical value indicating whether the coda samples generated in the meta-analysis are saved in the output object. By specifying the argument mcmc.samples = TRUE in nmadt.hierarchical, nmadt.hierarchical.MNAR, nmadt.hsroc, and nmadt.hsroc.MNAR, the MCMC samples are saved in the output objects. Users can then use packages such as mcmcse Reference Flegal, Hughes and Vats ⁴⁰ to determine whether the chains were run long enough to give sufficiently accurate estimates of posterior quantities.

The function nmadt.hierarchical returns a list with the raw output for graphing the results, the effect size estimates, including the posterior mean, standard deviation, median, and a $95$ % equal-tailed credible interval for the median.

Long description

The code consists of two main R commands.

Line 1: R prompt followed by set dot seed open parenthesis 10 close parenthesis.

Line 2 through 7: A multi-line function call assigning a value to the object kang dot out. The function is n m a d t dot hierarchical with the following arguments:

* n s t u equals 12

* K equals 2

* data equals dat dot kang

* directory equals tempdir open parenthesis close parenthesis

* testname equals c open parenthesis double quote D dash dimer double quote comma double quote Ultrasonography double quote close parenthesis

* diag equals 5

* off underscore diag equals 0.05

* digits equals 4

* mu underscore alpha equals 0

* mu underscore beta equals 0

* mu underscore eta equals 0

* preci underscore alpha equals 0.1

* preci underscore beta equals 0.1

* preci underscore eta equals 0.1

* n dot adapt equals 5000

* n dot iter equals 50000

* n dot chains equals 3

* conv dot diag equals TRUE

* trace equals double quote prev double quote

* dic equals TRUE

* m c m c dot samples equals FALSE

The following messages are output as the function runs:

Long description

The log begins with the line Start running M C M C dot dot dot followed by a series of indented processes.

Compiling data graph includes

* Resolving undeclared variables

* Allocating nodes

* Initializing

* Reading data back into data table

Compiling model graph includes

* Resolving undeclared variables

* Allocating nodes

Graph information provides the following data

* Observed stochastic nodes 48

* Unobserved stochastic nodes 18

* Total graph size 1898

Initializing model is followed by three progress bars. The first is filled with plus signs and marked 100 percent. The second and third are filled with asterisks and marked 100 percent.

The final section lists three concluding actions

* Start calculating M C M C convergence diagnostic statistics dot dot dot

* Start calculating deviance information criterion statistics dot dot dot followed by a progress bar of asterisks at 100 percent

* Start saving trace plots dot dot dot

When the proposed model is compiled, it requires an initial sampling phase during which the samplers adapt their behavior to maximize their efficiency. If the warning “adaptation incomplete” occurs, the users might need to increase n.adapt to make the process reach its maximum efficiency.

The results are saved in the list object kang.out. It contains the raw output for prevalence, sensitivity, specificity, positive and negative predictive values, positive likelihood, negative likelihood, and DIC. The trace plots are saved in the user-specified directory. The effect size names can be used to display the estimates. For example, the estimates of sensitivities and specificities (posterior mean and standard deviation, and posterior median and 95% credible intervals) can be displayed as

Long description

The document begins with the header Model type hierarchical. It is divided into four sections.

1. Sensitivity Mean plus or minus S D.

* D-dimer has a mean of 0.8376 and an S D of 0.0739.

* Ultrasonography has a mean of 0.9328 and an S D of 0.0466.

2. Sensitivity Median and 95 percent C r I.

* D-dimer has a median of 0.8451 with a 95 percent C r I of 0.6743 to 0.9540.

* Ultrasonography has a median of 0.9378 with a 95 percent C r I of 0.8277 to 0.9981.

3. Specificity Mean plus or minus S D.

* D-dimer has a mean of 0.8609 and an S D of 0.0675.

* Ultrasonography has a mean of 0.7606 and an S D of 0.1305.

4. Specificity Median and 95 percent C r I.

* D-dimer has a median of 0.8663 with a 95 percent C r I of 0.7101 to 0.9735.

* Ultrasonography has a median of 0.7724 with a 95 percent C r I of 0.4529 to 0.9658.

As shown above, each list for the effect sizes in the object kang.out consists of posterior means with sample standard deviations specified in the parenthesis and posterior medians with 95% equal tail credible interval. These results indicate that, on average, Ultrasonography shows higher sensitivity (median = 0.94, 95% CrI: 0.83–1.00) but lower specificity (median = 0.77, 95% CrI: 0.45–0.97) compared with D-dimer. This suggests that Ultrasonography is more effective in identifying true positive cases but may yield more false positives than D-dimer, reflecting a trade-off between sensitivity and specificity that is typical in diagnostic accuracy analyses.

In addition to these test-specific summaries, NMADTA now provides direct relative comparisons between tests. Specifically, the output includes pairwise differences in sensitivity and specificity (returned as $dSe and $dSp), which are computed by subtracting posterior samples across MCMC iterations and summarized using the same posterior Mean(SD) and Median(95% CrI) format. For the current example with two tests, kang.out$dSe and kang.out$dSp quantify the magnitude and uncertainty of the sensitivity and specificity differences between Ultrasonography and D-dimer. A positive value of $\Delta Se$ (or $\Delta Sp$ ) indicates higher sensitivity (or specificity) for the first-listed test, and a 95% CrI excluding zero provides strong posterior evidence of a meaningful difference under the fitted model.

Long description

The output is divided into two main sections based on R commands.

First section: R greater than kang dot out dollar d S e.

This section contains two sub-blocks for sensitivity data.

1. Under the header dollar Mean underscore S D, the row D-dimer vs Ultrasonography shows a Mean of negative 0.0952 and a S D of 0.0869.

2. Under the header dollar Median underscore C r I, the row D-dimer vs Ultrasonography shows a Median of negative 0.0910 and a 95 percent C r I ranging from negative 0.2798 to 0.0600.

Second section: R greater than kang dot out dollar d S p.

This section contains two sub-blocks for specificity data.

1. Under the header dollar Mean underscore S D, the row D-dimer vs Ultrasonography shows a Mean of 0.1003 and a S D of 0.1360.

2. Under the header dollar Median underscore C r I, the row D-dimer vs Ultrasonography shows a Median of 0.0934 and a 95 percent C r I ranging from negative 0.1379 to 0.4054.

In the example code, dic was specified as TRUE so the DIC statistic would be calculated. Users can access the DIC statistic and its components using

Trace plots for prevalence were generated because trace = "prev" was specified. Figure 1 shows the trace plots of the prevalence of DVT in the case study. The argument trace is a string vector describing the subset of quantities for which trace plots can be drawn, which can include prevalence ("prev"), sensitivity ("Se"), specificity ("Sp"), positive and negative predictive values ("ppv" and "npv", respectively), positive likelihood ("LRpos"), and negative likelihood ("LRneg"). Since we used the default n.chains = 3, three trace plots are drawn. Each trace plot shows evidence that the posterior samples of prevalence are drawn from the stationary distribution. In addition to visual inspection of trace plots, NMADTA provides a quantitative convergence diagnostic based on the PSRF proposed by Gelman and Rubin. When conv.diag = TRUE, PSRF values are computed for the monitored parameters and saved in a text file (e.g., ConvergenceDiagnostic.txt) in the user-specified directory. The PSRF compares between-chain and within-chain variability, and values approaching 1 indicate convergence. In practice, PSRF values below 1.1 are commonly regarded as evidence of acceptable convergence, while values closer to 1.0 (e.g., < 1.05) provide stronger reassurance of adequate mixing.Reference Gelman and Rubin ^{38 ,} Reference Brooks and Gelman ⁴¹ Upper confidence limits substantially above these thresholds may indicate lack of convergence and the need for additional iterations.

Figure 1

Trace plots generated by the R function nmadt.hierarchical for the prevalence of DVT in dat.kang

Figure 1 Long description

A multi-panel figure containing three vertically stacked line graphs labeled Chain 1, Chain 2, and Chain 3. All three graphs share identical axes. The horizontal x-axis represents Iterations, ranging from 0 to 50,000 with major ticks every 10,000. The vertical y-axis represents prevalence, ranging from 0.30 to 0.60.

* Chain 1 at the top shows a red jagged line fluctuating rapidly between 0.35 and 0.55. The trace appears stable with no significant upward or downward trend, indicating convergence.

* Chain 2 in the middle displays a similar red trace plot. The density of the fluctuations is consistent across the 50,000 iterations, centered around a prevalence of approximately 0.42.

* Chain 3 at the bottom mirrors the behavior of the previous two chains, with the red line oscillating within the same range of 0.30 to 0.60.

All three chains exhibit a fuzzy caterpillar appearance, which is characteristic of well-mixed Markov Chain Monte Carlo simulations.

4.2.2 Function nmadt.hierarchical.MNAR() under MNAR assumptions

In this context, MNAR refers specifically to outcome-dependent test performance availability within the evidence network. That is, the probability that a diagnostic test appears in a study (or that its accuracy results are reported) may depend on its underlying performance. This differs from standard MAR assumptions, under which test availability is assumed independent of unobserved accuracy parameters. While publication bias can be viewed as one manifestation of MNAR, our formulation more generally captures selection processes that operate at the test level within multi-test diagnostic studies. In the regular NMA-DT model, we assume that tests are missing from studies at random, but this assumption can be questionable in some situations. Sometimes tests can be MNAR, and we would need a different missingness mechanism for the meta-analysis. The function nmadt.hierarchical.MNAR permits a sensitivity analysis in which tests are MNAR. We define an $N\times K$ matrix $\mathbf {M}$ denoting the missingness status of the K candidate tests in each of the N studies. Entry $m_{ik}$ ( $i=1,\ldots , N$ , $k=1, \ldots , K$ ) of $\mathbf {M}$ indicates whether $T_k$ is missing in study i, that is, $m_{ik}=1$ if missing and $m_{ik}=0$ otherwise. We model $m_{ik}\sim Bern(p_{ik})$ , where $p_{ik}$ is the probability that $T_k$ is missing in study i, and specify

(4.1) $$ \begin{align} \begin{aligned} \operatorname{logit}(p_{ik}) = \gamma_{k}+\gamma_{1k}\times \operatorname{logit}(\text{Se}_{ik})+\gamma_{0k}\times \operatorname{logit}(\text{Sp}_{ik}), \end{aligned} \end{align} $$

where $\gamma _{1k}$ and $\gamma _{0k}$ represent the strength of association between missingness and the study-specific sensitivity and specificity, respectively. Note that $\gamma _{1k}$ and $\gamma _{0k}$ are regression coefficients in the logistic missingness model on the latent scale, rather than correlation coefficients; therefore, they are not constrained to lie within $[-1,1]$ .

The arguments of nmadt.hierarchical.MNAR are similar to those of nmadt.hierarchical. The major difference is that in nmadt.hierarchical.MNAR, users need to specify values of $\gamma _{1k}$ and $\gamma _{0k}$ , described above, as vectors of length K in the arguments gamma1 and gamma0. Users can also specify the prior mean and precision for $\gamma _k$ using mu_gamma and preci_gamma; the default is a vague $N(0,10)$ . The dataset dat.kang is used again as an example:

Long description

The code begins with the command set dot seed 10. The main function call assigns the result to kang M N A R dot out using the function nmadt dot hierarchical dot M N A R. The function parameters are listed across several lines starting with a plus sign.

Line 1 includes nstu equals 12, K equals 2, and data equals dat dot kang.

Line 2 defines testname as a vector containing D-dimer and Ultrasonography, and directory as tempdir.

Line 3 sets diag equals 5, off underscore diag equals 0.05, mu underscore alpha equals 0, mu underscore beta equals 0, and mu underscore eta equals negative 0.

Line 4 sets preci underscore alpha, preci underscore beta, and preci underscore eta all to 0.1.

Line 5 sets gamma 1 and gamma 0 to vectors of negative 0.5 and negative 0.5, and mu underscore gamma equals 0.

Line 6 sets preci underscore gamma equals 1, n dot adapt equals 5000, n dot iter equals 50000, and n dot chains equals 3.

Line 7 concludes with conv dot diag equals FALSE, trace equals NULL, dic equals FALSE, and mcmc dot samples equals FALSE.

For the dat.kang example, we set $\gamma _{1k}$ and $\gamma _{0k}$ to non-positive values to represent the assumption that $T_k$ is more likely to be missing if $T_k$ has low accuracy. In the MNAR missingness model, $\text {logit}(p_{ik})=\gamma _k+\gamma _{1k}\text {logit}(\mathrm {Se}_{ik})+\gamma _{0k}\text {logit}(\mathrm {Sp}_{ik})$ , $\gamma _{1k}$ and $\gamma _{0k}$ are logistic regression coefficients (not correlation parameters) and are therefore not constrained to $[-1,1]$ . A useful interpretation is that $\exp (\gamma _{1k})$ (resp. $\exp (\gamma _{0k})$ ) is the multiplicative change in the odds of missingness associated with a one-unit increase in $\text {logit}(\mathrm {Se}_{ik})$ (resp. $\text {logit}(\mathrm {Sp}_{ik})$ ), holding other terms fixed. We recommend treating $(\gamma _{1k},\gamma _{0k})$ as sensitivity-analysis parameters: $(0,0)$ corresponds to MAR with respect to accuracy, while negative values (e.g., $\{0,-0.5,-1,-2\}$ ) represent the common assumption that tests with poorer accuracy are more likely to be missing. Sensitivity analyses should be reported alongside routine convergence diagnostics. Here, more negative values imply a stronger MNAR dependence, that is, a larger increase in the odds of missingness as study-specific accuracy decreases on the logit scale. In practice, we recommend exploring a small grid, reporting how pooled sensitivity/specificity and relative test rankings change, and avoiding excessively large magnitudes that may push missingness probabilities toward 0/1 and slow MCMC mixing.

4.2.3 Functions nmadt.hsroc() and nmadt.hsroc.MNAR() for the HSROC model

Functions nmadt.hsroc and nmadt.hsroc.MNAR implement the HSROC models introduced previously, respectively, under the MAR assumption and the MNAR assumption, under which a test’s missingness is believed to be related to its accuracy. The arguments of these two functions are similar to those of nmadt.hierarchical. The arguments eta and xi_precision specify the prior mean and precision of $\log (\mathbf {S})$ and $\log (\mathbf {P})$ , the covariance matrices of the cutoff values and accuracy values, respectively.

MNAR mechanisms for the HSROC analysis use a missingness model similar to the one used for the hierarchical analysis. For nmadt.hsroc.MNAR, as for nmadt.hierarchical.MNAR, users need to specify values of $\gamma _{1k}$ and $\gamma _{0k}$ , described above, as vectors of length K in arguments gamma1 and gamma0; they can also specify a prior mean and precision for $\gamma _k$ other than the default vague $N(0,10)$ . Since the two functions are similar, and they give output in the same format as the analogous functions for the hierarchical models, we use the dataset dat.kang to illustrate the usage of nmadt.hsroc. The function is called as follows:

Long description

The code consists of two main commands.

First line: R prompt followed by set dot seed open parenthesis 10 close parenthesis.

Second line through sixth line: A multi-line function call assigned to the variable kang dot h s r o c.

- The function is n m a d t dot h s r o c with the following arguments:

- n s t u equals 12

- K equals 2

- data equals dat dot kang

- testname equals c open parenthesis double quote D dash dimer double quote comma double quote Ultrasonography double quote close parenthesis

- directory equals tempdir open parenthesis close parenthesis

- eta equals 0

- xi underscore preci equals 1.25

- digits equals 4

- n dot adapt equals 5000

- n dot iter equals 3000

- n dot chains equals 3

- conv dot diag equals FALSE

- trace equals NULL

- dic equals FALSE

- m c m c dot samples equals FALSE

Because the HSROC model is more complicated and has more parameters to estimate, the warning “adaptation incomplete” may occur more frequently. To avoid this warning, users can increase n.adapt.

It should be noted that the network diagnostic meta-analysis model is computationally intensive by design. This is primarily due to three factors. First, the model adopts a Bayesian hierarchical network structure and relies on MCMC sampling for inference, which requires joint sampling of high-dimensional random effects. Second, the incorporation of MNAR mechanisms introduces additional latent variables and regression parameters associated with the missingness process, substantially increasing model complexity. Third, jointly modeling multiple diagnostic tests across studies induces strong posterior dependence among parameters, which can slow down MCMC mixing and necessitate a larger number of iterations to achieve stable posterior inference. Consequently, longer runtime is generally expected for this model compared with standard diagnostic test meta-analysis models. For practical use, we recommend that users first run the model with a reduced number of MCMC iterations to verify data formatting and model behavior, and then increase the number of iterations for final analyses after convergence has been assessed. For illustration purposes and to reduce computation time, we use a smaller number of MCMC iterations (e.g., n.iter = 3000) in the example code. This setting is intended only to demonstrate the workflow and basic functionality. For real-world analyses, users should increase n.iter (and assess convergence diagnostics) to obtain stable posterior summaries. Corresponding analyses based on the HSROC model are provided in the Supplementary Material.