Once the posterior sampling is done, `metapack`

provides graphical tools to assess the model. Currently, a function for the *surface under the cumulative ranking (SUCRA) plot* is included. Model selection tools like the deviance information criterion (DIC) and log pseudo-marginal likelihood (LPML) are also available.

For a network meta-analysis, it is usually the case that ranking the treatments is of great interest. Based on the Bayesian posterior distributions, for each treatment \(t\) for \(t=1,\ldots, T\), \(P(t,r)\) defined as the probability that treatment \(t\) is ranked \(1\leq r\leq T\). Denoting by \(\mathbf{P}\) the discrete rank probability matrix of size \(T\times T\), it is easily observed that \(\mathbf{P}\) is row-stochastic. That is, $$ \sum_{r=1}^T P(t,r) = 1. $$ Calculating the cumulative probabilities by $$ F(t,x) = \sum_{r=1}^x P(t,r), $$ \(F(t,x)\) is the probability that the \(t\)-th treatment is ranked \(x\) or better. Since \(F(t,T)=1\) for every \(t\), the surface under the cumulative ranking distribution for \(t\)-th treatment is given by $$ \text{SUCRA}(t) = \dfrac{1}{T-1}\sum_{x=1}^{T-1}F(t,x). $$

`sucra`

To use `sucra`

or `plot.sucra`

, you must have an instance of `bayes.nmr`

. Let `fit`

be an object from `bayes.nmr`

.

```
ss <- sucra(fit)
plot(ss)
```

The function `sucra`

will return a list of class `sucra`

that contains the SUCRA, \(P(t,r)\) for all \((t,r)\), and the treatment labels given by the fitted object `fit`

. To directly handle the numerical values of the ranking probabilities, it should be stored in a separate object `ss`

, and then passed to the `plot`

function to generate the plots.

If you’re not interested in the intermediate values used to generate the plots, simply run the following:

```
plot(sucra(ss))
```

With multiple candidate models, it is essential to determine which one is the best. The deviance information criterion is a useful tool for such model selections. Let \(L_{oy}(\theta)\) be the observed-data likelihood from the response variable \(\boldsymbol{y}\) for which the subscript is omitted to accommodate arbitrary number of subscripts. Then, define the deviance function by \(D(\theta) = -2\log L_{oy}(\theta)\). Then, denoting the observed data by \(D_o\), $$ \mathrm{DIC} = D(\bar{\theta}) + 2 P_D,\quad P_D=\bar{D}-D(\bar{\theta}), $$ where \(\bar{D}=\mathrm{E}[D(\theta)\mid D_o]\) and \(\bar{\theta} = \mathrm{E}[\theta\mid D_o ]\).

`dic`

```
dic <- model.comp(fit, type = 'dic')
```

Akin to the DIC, the log pseudo-marginal likelihood (LPML) serves as another useful model selection tool. First, define the conditional predictive ordinate by

$$\mathrm{CPO}_k = \int f(\boldsymbol{y}_k \mid \theta)\pi(\theta\mid D_{oy}^{-k}, D_{os})\,\mathrm{d}\theta,$$

where \(D_{oy}^{-k}\) is \(D_{oy}\) with the \(k\)-th trial removed, and \(\pi(\theta\mid D_{oy}^{-k}, D_{os})\) is the posterior distribution. Then, the LPML is defined by $$ \mathrm{LPML} = \dfrac{1}{K}\sum_{k=1}^K \log (\mathrm{CPO}_k). $$

`lpml`

```
dic <- model.comp(fit, type = 'lpml')
```