An introduction to the Stochastic Search Inconsistency Factor Selection method using the ssifs package

Georgios Seitidis

Quick Revision

Stochastic Search Inconsistency Factor Selection (SSIFS) is the extension of Stochastic Search Variable Selection (SSVS) (George and McCulloch 1993) for identifying inconsistencies in Network Meta-Analysis (NMA). SSIFS is a two-step method, where in the first step inconsistency factors are specified, and in the second step, variable selection on inconsistency factors is performed using the SSVS method.

Inverse-variance NMA model

The inverse-variance random-effects NMA model adjusted to include $\ell =1, 2, \ldots, p$ inconsistency factors is described by the following equation $$\boldsymbol{y = X\mu + \beta +bZ + \epsilon}, \quad \boldsymbol{\epsilon} \sim N(\boldsymbol{0}, \boldsymbol{\Sigma} ) \quad and \quad \boldsymbol{\beta} \sim N(\boldsymbol{0}, \boldsymbol{\Delta})$$

where $\boldsymbol{y}$ denotes the treatments’ effect, $\boldsymbol{X}$ the design matrix, $\boldsymbol{\mu}$ the underlying basic contrasts, $\boldsymbol{\beta}$ the normally distributed random-effects, $\boldsymbol{b}$ the effect of inconsistency factors, $\boldsymbol{Z}$ the inconsistency factor’s index matrix, and $\boldsymbol{\epsilon}$ the normally distributed sampling errors. Correlation matrix $\boldsymbol{\Delta}$ is a block diagonal matrix, assuming common heterogeneity across treatment comparisons, while covariance matrix $\boldsymbol{\Sigma}$ is assumed known and obtained from the data based on Franchini et al. (2012). Matrix $\boldsymbol{Z}$ contains as elements values 1, -1 and 0, indicating in which comparisons inconsistency factor is added. Among the choices that may be considered for the specification of the $\boldsymbol{Z}$ matrix are the Lu and Ades model (Lu and Ades 2006), the design-by-treatment model (Higgins et al. 2012), and the random-effects implementation of the design-by-treatment model (Jackson et al. 2014).

Variable Selection

In SSIFS the effect of an inconsistency factor $\ell$ is described from a mixture of two normal distributions, which can be written as $$b_{\ell} | \gamma_{\ell} \sim (1-\gamma_{\ell}) N(0, \psi_{\ell}^{2}) + \gamma_\ell N(0, c^{2}\psi_\ell^{2})$$ where $b_{\ell}$ is the effect of the inconsistency factor, $\gamma_{\ell}$ is an auxiliary variable indicating if the inconsistency factor is included in the NMA model, and $c$, $\psi_{\ell}$ are tuning parameters controlling the mixing ability of the method.

In matrix notation SSIFS is written as $$\boldsymbol{b | \gamma} \sim N( \boldsymbol{0}, \boldsymbol{D_\gamma R D_\gamma} )$$ where $$\boldsymbol{D}_{\boldsymbol{\gamma}}= \begin{pmatrix} a_1 \psi_1 & 0 & \dots & 0 \\ 0 & a_2 \psi_2 & \dots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \dots & a_\text{p} \psi_\text{p} \end{pmatrix} ,\qquad a_\ell= \begin{cases} 1, \quad \gamma_\ell=0\\ c, \quad \gamma_\ell=1 \end{cases}.$$ Matrix $\boldsymbol{R}$ denotes the prior correlation between the inconsistency factors. We can assume that inconsistency factors are independent by setting $\boldsymbol{R} = \boldsymbol{I}$, or we can assume a dependency between inconsistency factors by using a Zellner g-prior as described bellow $$\boldsymbol{R}=g \boldsymbol{ (Z^{'}Z) } \sigma^{2}, \quad \pi (\sigma^{2} ) \propto \frac{1}{\sigma^2}.$$ For the specification of the parameter $g$, the unit information criterion (Kass and Wasserman 1995) is used, which translates in SSIFS to the total number of observed comparisons in the network.

Prior inclusion probabilities of the inconsistency factors are specified as $$\gamma_\ell \sim Bernoulli(1 - \pi_{con}^{ \frac{1}{p}}), \quad \ell = 1, 2, \ldots, p$$ where $\pi_{con}$ is the probability to have a consistent network and reflects our prior believes on how likely is to have a consistent network. In a review of 201 networks, 44 networks were found to be globally inconsistent (Veroniki et al. 2021). Thus, $\pi_{con} \sim Beta(157, 44)$ is proposed.

Tuning

Tuning is crucial in SSIFS in order to ensure a good mixing of the method. Ideally, the effect of an inconsistency factor when it is included in the NMA model should lie in an area close to zero, and far away from this area when it is not included in the NMA model. Regarding parameter $c$, values between 10 and 100 usually perform well in most cases (George and McCulloch 1993; Perrakis and Ntzoufras 2015). Possible values of parameter $\psi_\ell$ could be obtained from a pilot MCMC run of the NMA model as the standard deviation of the inconsistency factors.

Minimum value of inconsistency

By properly tuning parameters $c$ and $\psi_\ell$, a difference between direct and indirect evidence that is of practical significance (say $\omega$) can be defined. Thus, an inconsistency factor with a coefficient larger than $\omega$ in absolute values ($|b_{\ell}| > \omega$), should be included in the NMA model ($\gamma_\ell = 1$). In the case where $\boldsymbol{R = I}$, the inconsistency factor will have higher probability to be included in the NMA model when $$|b_\ell |> \psi_\ell \sqrt{\xi(c)}, \quad \xi(c) = \frac{2c^2 \log{c}}{c^2 -1}.$$ For example, if a difference above 0.2 is considered important, one possible parameterization is to set $c = 10$ and $\psi_\ell = \frac{0.2}{\sqrt{\xi(10)}} \approx 0.1.$

Inconsistency Detection

Inconsistency in SSIFS is evaluated by examining the posterior inclusion probabilities of the inconsistency factors, the posterior model probabilities, the posterior model odds and the Bayes factor of the consistent NMA model over the inconsistent NMA model.

Posterior Inclusion probabilities

Posterior inclusion probabilities estimated as the average of times the inconsistency factor was included in the NMA model in the MCMC draws. Estimates above 0.5 indicates local inconsistency, which cause global inconsistency to the network.

Posterior Model Odds

Posterior model odds are obtained as the ratio of the posterior model probabilities which are estimated as $$f\hat{(m|\boldsymbol{y})}=\frac{1}{M-B} \sum_{t=B+1}^{M}{I(m^{(t)}=m)}, \quad m(\boldsymbol{\gamma})=\sum_{\ell = 1}^{p}{\gamma_\ell 2^{\ell-1}}$$ where $M$ is the number of MCMC iterations, $B$ the burn-in period and $m^{(t)}$ a model indicator which transforms the $\boldsymbol{\gamma}$ to a unique decimal number. By examined the posterior odds of the consistent NMA model ($m(\boldsymbol{\gamma}) = 0$) over the inconsistent NMA models, we can evaluate the consistency assumption. Also, the comparison between the consistent NMA model ($m(\boldsymbol{\gamma}) = 0$) over all the other observed inconsistent NMA models ($m(\boldsymbol{\gamma}) \ne 0$), indicates if the NMA model is globally consistent.

Implementation through ssifs

Installation

You can install the development version of ssifs like so:

install.packages("devtools")
devtools::install_github("georgiosseitidis/ssifs")

Data

ssifs requires the contrast-based data used for the NMA model. Also, in the multi-arm studies, all possible comparisons must be provided. In the case where the network is disconnected, ssifs keeps only those studies that belong to the largest sub-network in order to maintain one connected network.

Example

Load the brief alcohol intervention dataset from the ssifs package. The dataset is from a published NMA (Seitidis et al. 2022; Hennessy et al. 2019) and contains 37 studies evaluating the comparative effectiveness of brief alcohol interventions on preventing hazardous drinking in college students.

library(ssifs)
data("Alcohol", package = "ssifs")

Prepare the data for the ssifs.

TE <- Alcohol$TE
seTE <- Alcohol$seTE
study <- Alcohol$studyid
treat1 <- Alcohol$treat1
treat2 <- Alcohol$treat2

Run the ssifs function, using AO-CT as a reference intervention.

set.seed(12)
m <- ssifs(TE = TE, seTE = seTE, studlab = study, treat1 = treat1, treat2 = treat2, ref = "AO-CT")

The function by default for the specification of the matrix $\boldsymbol{Z}$ uses the design-by-treatment model. You can use the Lu & Ades model by setting the argument method = "LuAdes", or the random-effects implementation of the design-by-treatment model by setting method = "Jackson". Also, the function by default specifies the correlation matrix $\boldsymbol{R}$ by using a Zellner g-prior. You can assume that inconsistency factors are independent by setting the argument zellner = FALSE.

Regarding the prior inclusion probabilities, the function assumes that $\pi_{con} \sim Beta(157, 44)$ (argument rpcon = TRUE). By setting the argument rpcon = FALSE you can set the probability of $\pi_{con}$ fixed. If rpcon = FALSE, the function assumes that $\pi_{con} = 0.5$. You can change this probability from the argument pcons.

Detection of inconsistency

Posterior inclusion probabilities

Posterior inclusion probabilities can be obtained like so:

m$Posterior_inclusion_probabilities
#>            Comparison            Design    PIP       b    b.lb   b.ub
#> 1     Alc101 ; BASICS      Alc101BASICS 0.0234  0.0006 -0.1243 0.1480
#> 2      AO-CT ; Alc101       AO-CTAlc101 0.0208  0.0027 -0.1354 0.2078
#> 3     e-CHUG ; BASICS      e-CHUGBASICS 0.0236 -0.0011 -0.1974 0.1716
#> 4  Active-CT ; THRIVE   Active-CTTHRIVE 0.0237 -0.0019 -0.1798 0.0978
#> 5      AO-CT ; AlcEdu AO-CTAlcEdue-CHUG 0.0226  0.0021 -0.1023 0.1670
#> 6      AO-CT ; e-CHUG AO-CTAlcEdue-CHUG 0.0241  0.0015 -0.1064 0.1527
#> 7      AO-CT ; e-CHUG AO-CTe-CHUGBASICS 0.0216 -0.0016 -0.1522 0.1170
#> 8      AO-CT ; BASICS AO-CTe-CHUGBASICS 0.0235  0.0034 -0.0808 0.1778
#> 9      AO-CT ; AlcEdu AlcEduAO-CTAlc101 0.0242 -0.0024 -0.1250 0.0724
#> 10     AO-CT ; Alc101 AlcEduAO-CTAlc101 0.0204  0.0035 -0.1172 0.2037

The first two columns refer to comparisons where inconsistency factors are added. For example, the first row refers to the inconsistency factor that added to the comparisons between interventions Alc101 and BASICS, obtained from the two-arm studies that compare these interventions. The fifth row refers to the inconsistency factor that added to the comparison between interventions AO-CT and AlcEdu, obtained from the multi-arm studies that compare the interventions AO-CT, AlcEdu and e-CHUG. Column PIP refers to posterior inclusion probability, while columns b, b.lb and b.ub to the inconsistency factors effect estimates with the corresponding 95% credible interval.

If method = "LuAdes", the column Design is NA. This is because the Lu & Ades model accounts only for loop inconsistencies.

In this example, the posterior inclusion probabilities suggest that the network is globally and locally consistent, since there are not any significant local inconsistencies that causes global inconsistency to the network (posterior inclusion probabilities $\approx 0 < 0.5$). Also, note that the corresponding effect estimates are not significant and close to zero.

Posterior model odds

The posterior model odds can be obtained like so:

head(m$Posterior_Odds)
#>                                IFs  Freq f(m|y) PO_IFCONS
#> 1                           No IFs 63573 0.7947    1.0000
#> 2 AO-CT ; AlcEdu_AlcEduAO-CTAlc101  1556 0.0194   40.8567
#> 3               Active-CT ; THRIVE  1547 0.0193   41.0944
#> 4                  e-CHUG ; BASICS  1536 0.0192   41.3887
#> 5 AO-CT ; BASICS_AO-CTe-CHUGBASICS  1532 0.0192   41.4967
#> 6 AO-CT ; e-CHUG_AO-CTAlcEdue-CHUG  1522 0.0190   41.7694

Column IFs refers to the model observed in the MCMC draws (Inconsistency factors are separated by the symbol ,), Freq refers to the number of times the model was observed in the MCMC draws, f(m|y) denotes the posterior model probability and PO_IFCONS to the posterior odds of the consistent NMA model (NO IFs) over the model in the corresponding row. An estimate over 1 favors the consistent NMA model.

In this example, the first row refers to the consistent NMA model (NO IFs), showing that the posterior model probability of the consistent NMA model is 0.79. The posterior odds is 1 as expected, since $$PO = \frac{f(m_{con}|\boldsymbol{y})}{f(m_{con}|\boldsymbol{y})}=\frac{0.7943}{0.7943} = 1.$$ The third row refers to the inconsistent NMA model (say $m_{3}$) where inconsistency factor was added in the comparisons between interventions AO-CT and e-CHUG, obtained from the multi-arm studies that compare the interventions AO-CT, e-CHUG and BASICS. The posterior odds calculated as $$PO_{m_{con}m_{3}} = \frac{f(m_{con}|\boldsymbol{y})}{f(m_{3}|\boldsymbol{y})} =\frac{0.7942875}{0.0193125} = 41.1282.$$

By looking the posterior model probabilities and the posterior model odds, we conclude that the consistent NMA is the most dominant model since $\hat{f(m|y)}=0.79$. Also, the posterior odds clearly favors the consistent NMA model. Thus, we conclude that the network is both globally and locally consistent, because significant local inconsistencies that causes globally inconsistency to network were not observed.

Global test of inconsistency

The global test of inconsistency is conducted by calculating the Bayes factor of the consistent NMA model over the rest inconsistent NMA models. Thus, the posterior model probabilities of the inconsistent NMA models are summed. An estimate above 1 indicates that the network is globally consistent.

The global test of inconsistency can be obtained like so:

m$Bayes_Factor
#> [1] 1.0783

In the example, the Bayes factor estimated above 1, suggesting that the network is globally consistent.

Model’s mixing ability

You can test the mixing ability of the ssifs model like so:

spike.slab(m)

References

Franchini, Angelo J, Sofia Dias, Anthony E Ades, Jeroen P Jansen, and Nicky J Welton. 2012. “Accounting for Correlation in Network Meta-Analysis with Multi-Arm Trials.” Research Synthesis Methods 3 (2): 142–60.

George, Edward I, and Robert E McCulloch. 1993. “Variable Selection via Gibbs Sampling.” Journal of the American Statistical Association 88 (423): 881–89.

Hennessy, Emily Alden, Emily E Tanner-Smith, Dimitris Mavridis, and Sean P Grant. 2019. “Comparative Effectiveness of Brief Alcohol Interventions for College Students: Results from a Network Meta-Analysis.” Prevention Science 20 (5): 715–40.

Higgins, JPT, D Jackson, JK Barrett, G Lu, AE Ades, and IR White. 2012. “Consistency and Inconsistency in Network Meta-Analysis: Concepts and Models for Multi-Arm Studies.” Research Synthesis Methods 3 (2): 98–110.

Jackson, Dan, Jessica K Barrett, Stephen Rice, Ian R White, and Julian PT Higgins. 2014. “A Design-by-Treatment Interaction Model for Network Meta-Analysis with Random Inconsistency Effects.” Statistics in Medicine 33 (21): 3639–54.

Kass, Robert E, and Larry Wasserman. 1995. “A Reference Bayesian Test for Nested Hypotheses and Its Relationship to the Schwarz Criterion.” Journal of the American Statistical Association 90 (431): 928–34.

Lu, Guobing, and AE Ades. 2006. “Assessing Evidence Inconsistency in Mixed Treatment Comparisons.” Journal of the American Statistical Association 101 (474): 447–59.

Perrakis, Konstantinos, and Ioannis Ntzoufras. 2015. “Stochastic Search Variable Selection (SSVS).” In Wiley StatsRef: Statistics Reference Online, 1–6. John Wiley & Sons, Ltd.

Seitidis, G, S Nikolakopoulos, EA Hennessy, EE Tanner-Smith, and D Mavridis. 2022. “Network Meta-Analysis Techniques for Synthesizing Prevention Science Evidence.” Prevention Science 23 (3): 415–24.

Veroniki, Areti Angeliki, Sofia Tsokani, Ian R White, Guido Schwarzer, Gerta Rücker, Dimitris Mavridis, Julian Higgins, and Georgia Salanti. 2021. “Prevalence of Evidence of Inconsistency and Its Association with Network Structural Characteristics in 201 Published Networks of Interventions.” BMC Medical Research Methodology 21 (1): 1–10.