Model Assessment and Treatment Ranking

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.

Surface under the cumulative ranking (SUCRA)

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

Example code - `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))

Deviance Information Criterion

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

Example - `dic`

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

Logarithm of the Pseudo-marginal Likelihood

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

Example - `lpml`

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

Docs