---
title: "crossnma: Cross-Design & Cross-Format Network Meta-Analysis and Regression"
author: "Tasnim Hamza, Guido Schwarzer and Georgia Salanti"
output:
  rmarkdown::pdf_document:
   number_sections: true
  rmarkdown::html_vignette:
   toc: true
   number_sections: true
bibliography: references.bib 
vignette: >
  %\VignetteIndexEntry{crossnma: Cross-Design & Cross-Format Network Meta-Analysis and Regression}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  message = TRUE,
  warning = FALSE
  )
options(knitr.kable.NA = ".")
```

```{r setup, message = FALSE, echo = FALSE}
library("crossnma")
set.seed(1910)
settings.meta(digits = 3)
cilayout("(", " to ")
```


# Introduction

In network meta-analysis we synthesize all relevant available evidence
about health outcomes from competing treatments. That evidence might
come from different study designs and in different formats: from
non-randomized studies (NRS) or randomized controlled trials (RCT) as
individual participant data (IPD) or as aggregate data (AD). We set up
the package **crossnma** to synthesize all available evidence for a
binary outcome with the odds ratio as effect measure.

This document demonstrates how to use **crossnma** to synthesize
cross-design evidence and cross-format data via Bayesian network
meta-analysis and meta-regression (NMA and NMR). All models are
implemented in JAGS [@plummer_jags].

We describe the workflow within the package using a worked example
from a network meta-analysis of studies for treatments in relapsing
remitting multiple sclerosis (RRMS). The primary outcome is the
occurrence of relapses in two years (binary outcome). In the analysis,
the relative effect will be the odds ratio (OR). The aim is to compare
the efficacy of four treatments using the data from 6 different
studies in different formats and different designs.

# The synthesis models

We first introduce the model that synthesizes studies with
individual-level (IPD) or/and aggregate data (AD) ignoring their
design (unadjusted synthesis). Then, we present three possible models
that account for the different study designs. In the table below we
set the notation that will be used in the description of the four
synthesis models.

| Notation | Description | Argument in `crossnma.model()` | 
|:---------- |:---------- | :---------- |
|$i=1, ..., np_j$ | participant id|  |
|$j=1, ..., ns$ | study id| `study` |
|$k=1, ..., K$ | treatment index| `trt` |
|$ns_{IPD}, ns_{AD}, ns_{RCT}, ns_{NRS}$| the number of studies. The index refers to the design or format of the study|  |
|$y_{ijk}$ | binary outcome (0/1)| `outcome` |
|$p_{ijk}$ | probability of the event to occur| |
|$r_{jk}$ | the number of events per arm| `outcome` |
|$n_{jk}$ | the sample size per arm| `n` |
|$b$ |the study-specific reference|*|
|$u_{jb}$ | The treatment effect of the study-specific reference $b$  when $x_{ijk}=\bar{x}_{j}=0$ | |
|$\delta_{jbk}$|log(OR) of treatment $k$ relative to $b$||
|$x_{ijk}$|the covariate|`cov1`, `cov2`, `cov3`|
|$\bar{x}_{j}$|the mean covariate for study $j$||
|$d_{Ak}$| the basic parameters. Here, $d_{AA}=0$ when A is set as the reference in the network|use `reference` to assign the reference treatment|
|$z_j$| study characteristics to estimate the bias probability $\pi_j$| `bias.covariate` |
|$w$| common inflation factor of variance for the NRS estimates | the element `var.infl` in `run.nrs`|
|$\zeta$| common mean shift of the NRS estimates | the element `mean.shift` in `run.nrs`|
*The study-specific reference $b$ is assigned automatically to be the
network reference for studies that have the network reference
treatment. If not, it is assigned to the first alphabetically ordered
treatment on the study.


## Unadjusted network meta-regression (NMR)

We synthesize the evidence from RCT and NRS without acknowledging the
differences between them. We combine the IPD data from RCT and NRS in
one model and we do the same in another model with the AD
information. Then, we combine the estimates from both parts as
described in Section 2.5.

**NMR model for IPD studies**

$$
y_{ijk} \sim Bernoulli(p_{ijk})
$$
\begin{equation}
  logit(p_{ijk}) =
    \begin{cases}
      u_{jb} +\beta_{0j} x_{ijk} & \text{if $k=b$}\\
      u_{jb} +\delta_{jbk} + \beta_{0j}x_{ijk}+\beta^w_{1,jbk}x_{ijk} +
 (\beta^B_{1,jbk}-\beta^w_{1,jbk}) \bar{x}_{j} & \text{if $k\ne b$}
    \end{cases}       
\end{equation}

**NMR model for AD studies**

$$
r_{jk} \sim Binomial(p_{.jk},n_{jk})
$$
\begin{equation}
  logit(p_{.jk}) =
    \begin{cases}
      u_{jb}  & \text{if $k=b$}\\
      u_{jb} +\delta_{jbk} +\beta^B_{1,jbk} \bar{x}_{j} & \text{if $k\ne b$}
    \end{cases}       
\end{equation}


## Using non-randomized studies (NRS) to construct priors for the treatment effects

First, the (network) meta-regression with only NRS data estimates the
relative treatment effects with posterior distribution of mean
$\tilde{d}^{NRS}_{Ak}$ and variance $V^{NRS}_{Ak}$ (use `run.nrs` in
`crossnma.model()` to control this process). The posteriors of NRS
results are then used as priors for the corresponding basic parameters
in the RCT model, $d_{Ak} \sim
\mathcal{N}(\tilde{d}^{NRS}_{Ak},V^{NRS}_{Ak})$. We can adjust for
potential biases associated with NRS by either shifting the mean of
the prior distribution with a bias term $\zeta$ or by dividing the
prior variance with a common inflation factor $w, 0<w<1$ controls NRS
contribution. The assigned priors become $d_{Ak} \sim
\mathcal{N}(\tilde{d}^{NRS}_{Ak}+\zeta,V^{NRS}_{Ak}/w)$.


## Bias-adjusted model 1

We incorporate judgments about study risk of bias (RoB) in
bias-adjusted model 1 and model 2. Each judgment about the risk of
bias in a study is summarized by the index $R_j$ which takes binary
values 0 (no bias) or 1 (bias). In bias-adjusted model 1, we extend
the method introduced by @dias_2010 by adding a treatment-specific
bias term $\gamma_{2,jbk} R_j$ to the relative treatment effect on
both the AD and IPD parts of the model. A multiplicative model can
also be employed, where treatment effects are multiplied by
$\gamma_{1,jbk}^{R_j}$. We can add either multiplicative bias effects,
additive bias effects, or both (in this case, $\delta_{jbk}$ should be
dropped from the additive part). The models in previous section are
extended to adjust for bias as follows.

**NMR model for IPD studies**

\begin{equation}
  logit(p_{ijk}) =
    \begin{cases}
      u_{jb} +\beta_{0j} x_{ijk} & \text{if $k=b$}\\
      u_{jb} +\overbrace{\delta_{jbk} \gamma_{1,jbk}^{R_j}}^{\text{multiplicative}}+\overbrace{\delta_{jbk}+\gamma_{2,jbk} R_j}^{\text{additive}}+ \beta_{0j}x_{ijk}+\beta^w_{1,jbk} x_{ijk}+
 (\beta^B_{1,jbk}-\beta^w_{1,jbk}) \bar{x}_{j} & \text{if $k\ne b$}
    \end{cases}       
\end{equation}

**NMR model for AD studies**

\begin{equation}
  logit(p_{.jk}) =
    \begin{cases}
      u_{jb}  & \text{if $k=b$}\\
      u_{jb} +\overbrace{\delta_{jbk} \gamma_{1,jbk}^{R_j}}^{\text{multiplicative}}+\overbrace{\delta_{jbk}+\gamma_{2,jbk} R_j}^{\text{additive}}+
 \beta^B_{1,jbk} \bar{x}_{j} & \text{if $k\ne b$}
    \end{cases}       
\end{equation}

The bias indicator $R_j$ follows the following distribution

$$
R_j \sim Bernoulli(\pi_j)
$$
The bias probabilities $\pi_j$ are study-specific and can be estimated in two different ways. They are either given informative beta priors (${Beta(a_1,a_2)}$) that are set according to the risk of bias for each study.
$$
\pi_j \sim Beta(a_1, a_2)
$$

The hyperparameters $a_1$ and $a_2$ should be chosen in a way that
reflects the risk of bias for each study. The degree of skewness in
beta distribution can be controlled by the ratio $a_1/a_2$ . When
$a_1/a_2$ equals 1 (or $a_1=a_2$), there is no skewness in the beta
distribution (the distribution is reduced to a uniform distribution),
which is appropriate for studies with unclear risk of bias. When the
ratio $a_1/a_2$ is closer to 1, the more the mean of probability of
bias (expected value of $\pi_j=a_1/(a_1+a_2))$ gets closer to 1 and
the study acquires 'major' bias adjustment. The default beta priors
are as follows: high bias RCT `prior.pi.high.rct='dbeta(10, 1)'`, low
bias RCT `prior.pi.low.rct = 'dbeta(1, 10)'`, high bias NRS
`prior.pi.high.nrs = 'dbeta(30, 1)'` and low bias NRS
`prior.pi.low.nrs = 'dbeta(1, 30)'`. Alternatively, we can use the
study characteristics $z_j$ to estimate $\pi_j$ through a logistic
transformation (internally coded).

We combine the multiplicative and the additive treatment-specific bias
effects across studies by assuming they are exchangeable
$\gamma_{1,jbk}\sim \mathcal{N}(g_{1,bk},\tau_{1,\gamma}^2
)$,$\gamma_{2,jbk}\sim \mathcal{N}(g_{2,bk},\tau_{2,\gamma}^2 )$) or
common $\gamma_{1,jbk}=g_{1,bk}$ and
$\gamma_{2,jbk}=g_{2,bk}$. @dias_2010 proposed to model the mean bias
effect $(g_{1,bk}, g_{2,bk})$ based on the treatments being compared.

\begin{equation}
  g_{m,bk} =
    \begin{cases}
      g_m  & \text{if $b$ is inactive treatment}\\
      0 \text{ or } (-1)^{dir_{bk}} g_m^{act} & \text{if $b$ and $k$ are active treatments}
    \end{cases}
\end{equation}
where $m={1,2}$. This approach assumes a common mean bias for studies
that compare active treatments with an inactive treatment (placebo,
standard or no treatment). For active vs active comparisons, we could
assume either a zero mean bias effect or a common bias effect
$g_m^{act}$. The direction of bias $dir_{bk}$ in studies that compare
active treatments with each other should be defined in the data. That
is set to be either 0, meaning that bias favors $b$ over $k$, or 1 ,
meaning that $k$ is favored to $b$. In `crossnma.model()`, the bias
direction is specified by providing the unfavoured treatment for each
study, `unfav`. To select which mean bias effect should be applied,
the user can provide the `bias.group` column as data. Its values can
be 0 (no bias adjustment), 1 (to assign for the comparison mean bias
effect $g_m$) or 2 (to set bias $g_m^{act}$).

Another parameterisation of the logistic model with additive bias
effect is

**NMR model for IPD studies**

\begin{equation}
  logit(p_{ijk}) =
    \begin{cases}
      u_{jb} +\beta_{0j} x_{ijk} & \text{if $k=b$}\\
      u_{jb} +(1-R_j)\delta_{jbk}+\delta_{jbk}^{bias}R_j+ \beta_{0j}x_{ijk}+\beta^w_{1,jbk} x_{ijk}+
 (\beta^B_{1,jbk}-\beta^w_{1,jbk}) \bar{x}_{j} & \text{if $k\ne b$}
    \end{cases}
\end{equation}

**NMR model for AD studies**

\begin{equation}
  logit(p_{.jk}) =
    \begin{cases}
      u_{jb}  & \text{if $k=b$}\\
      u_{jb} +(1-R_j)\delta_{jbk}+\delta_{jbk}^{bias}R_j+
 \beta^B_{1,jbk} \bar{x}_{j} & \text{if $k\ne b$}
    \end{cases}
\end{equation}

Then the bias-adjusted relative treatment effect
($\delta_{jbk}^{bias}=\delta_{jbk}+\gamma_{jbk}$) can be assumed
exchangeable across studies $\delta_{jbk}^{bias} \sim
\mathcal{N}(g_{bk}+d_{Ak}-d_{Ab}, \tau^2 /q_j)$ or fixed as
$\delta_{jbk}^{bias}=g_{bk} + d_{Ak}-d_{Ab}$. In this
parameterisation, instead of assigning prior to the between-study
heterogeneity in bias effect $\tau_{\gamma}$, we model the RoB weight
$q_j=\tau^2/(\tau^2+\tau_{\gamma}^2$) for each study. This quantity
$0<q_j<1$ quantifies the proportion of the between-study heterogeneity
that is not explained by accounting for risk of bias. The values of
$v$ determine the extent studies at high risk of bias will be
down-weighted on average. Setting $v=1$ gives $E(q_j )=v/(v+1)=0.5$,
which means that high risk of bias studies will be penalized by 50% on
average. In `crossnma.model()`, the user can assign the average
down-weight $E(q_j )$ to the argument `down.wgt`.


## Bias-adjusted model 2

Another way to incorporate the RoB of the study is by replacing
$\delta_{jbk}$ by a "bias-adjusted" relative treatment effect
$\theta_{jbk}$. Then $\theta_{jbk}$ is modeled with a bimodal normal
distribution as described in Section 2.5. For more details see
@verde_2020.

**NMR model for IPD studies**

\begin{equation}
  logit(p_{ijk}) =
    \begin{cases}
      u_{jb} +\beta_{0j} x_{ijk} & \text{if $k=b$}\\
      u_{jb} +\theta_{jbk} + \beta_{0j} x_{ijk}+\beta^w_{1,jbk} x_{ijk}+
 (\beta^B_{1,jbk}-\beta^w_{1,jbk}) \bar{x}_{j} & \text{if $k\ne b$}
    \end{cases}       
\end{equation}

**NMR model for AD studies**

\begin{equation}
  logit(p_{jk}) =
    \begin{cases}
      u_{jb} & \text{if $k=b$}\\
      u_{jb} +\theta_{jbk} +\beta^B_{1,jbk} \bar{x}_{j} & \text{if $k\ne b$}
    \end{cases}       
\end{equation}

where the bias-adjusted relative treatment effect ($\theta_{jk}$) are
modeled via random-effects model with a mixture of two normal
distributions.

$$
\theta_{jbk} \sim (1-\pi_j) \mathcal{N}(d_{Ak}-d_{Ab}, \tau^2) +  \pi_j \mathcal{N}(d_{Ak}-d_{Ab}+\gamma_{jbk}, \tau^2+\tau_\gamma^2)
$$

Alternatively, we can summarize these relative effects assuming a
common-effect model

$$
\theta_{jbk}= d_{Ak}-d_{Ab}+\pi_j \gamma_{jbk}
$$


## Assumptions about the model parameters

The table below summarizes the different assumptions implemented in
the package about combining the parameters in the models described
above.

| Parameter 	| Assumptions| Argument in `crossnma.model()`|  	
|:---------- |:---------- | :------------ |
|Relative treatment effect ($\delta_{jbk}$)| Random-effects: $\delta_{jbk}\sim \mathcal{N}(d_{Ak}-d_{Ab}, \tau^2)$| `trt.effect='random'` |
| |Common-effect: $\delta_{jbk}=d_{Ak}-d_{Ab}$| `trt.effect='common'`|
|Covariate effect ($\beta_{0j}$) | Independent effects: $\beta_{0j} \sim \mathcal{N}(0, 10^2)$| `reg0.effect='independent'` |
|  |Random-effects: $\beta_{0j} \sim \mathcal{N}(B_0, \tau^2_{0})$| `reg0.effect='random'`|
|Within-study covariate-treatment | Independent effects: $\beta_{1,jbk}^W \sim \mathcal{N}(0, 10^2)$| `regw.effect='independent'` |
|interaction ($\beta_{1,jbk}^W$) | Random-effects: $\beta_{1,jbk}^W \sim \mathcal{N}(B_{1,Ak}^W-B_{1,Ab}^W, \tau^2_{W})$| `regw.effect='random'` |
| 	|Common-effect: $\beta_{1,jbk}^W = B_{1, Ak}^W-B_{1, Ab}^W$| `regw.effect='common'`|
|Between-study covariate-treatment| Independent effects: $\beta_{1,jbk}^B \sim \mathcal{N}(0, 10^2)$| `regb.effect='independent'` |
|interaction ($\beta_{1,jbk}^B$) | Random-effects: $\beta_{1,jbk}^B \sim \mathcal{N}(B_{1, Ak}^B-B_{1, Ab}^B, \tau_B^2)$| `regb.effect='random'` |
|  	|Common-effect: $\beta_{1,jbk}^B = B_{1, Ak}^B-B_{1, Ab}^B$| `regb.effect='common'`|
|Bias effect ($\gamma_{m,jbk}$), $m={1,2}$| Random-effects: $\gamma_{m,jbk} \sim \mathcal{N}(g_{m, bk}, \tau_{m,\gamma}^2)$| `bias.effect='random'` |
| 	|Common-effect: $\gamma_{m,jbk}=g_{m,bk}$| `bias.effect='common'`|
|Mean bias effect $g_{m,bk}$|The treatment $k$ is active. \ $g_{m,bk}=g_m$ ($b$ inactive), \ $g_{m,bk}=0$ ($b$ active \& no bias) \ $g_{m,bk}=g_m^{act}$($b$ active \& bias)| `unfav=0`, `bias.group=1` \ `unfav=1`, `bias.group=0` \ `unfav=1`, `bias.group=2`|
|Bias probability ($\pi_j$)| $\pi_j \sim Beta(a_1,a_2)$| `pi.high.nrs`, `pi.low.nrs`, `pi.high.rct`, `pi.low.rct`|
|| $\pi_j = e+fz_j$| `bias.covariate`|


# Synthesis of studies comparing drugs for relapsing-remitting multiple sclerosis


## Description of the data

The data we use are fictitious but resemble real RCTs with IPD and
aggregate data included in @Tramacere15. The studies provide either
individual participant data `ipddata` (3 RCTs and 1 cohort study) or
aggregate data `stddata` (2 RCTs). In total, four drugs are compared
which are anonymized.

The `ipddata` contains 1944 participants / rows. We display the first
few rows of the data set:

```{r}
dim(ipddata)
head(ipddata)
```

For each participant, we have information for the `outcome` relapse (0
= no, 1 = yes), the treatment label `treat`, the `age` (in years) and
`sex` (0 = female, 1 = male) of the participant. The following columns
are set on study-level (it is repeated for each participant in each
study): the `id`, the `design` of the study (needs to be either
`"rct"` or `"nrs"`), the risk of bias `rob` on each study (can be set
as low, high or unclear), the `year` of publication, the `bias.group`
for the study comparison and the study unfavoured treatment
`unfavored`.

The aggregate meta-analysis data must be in long arm-based format with
the exact same variable names and an additional variable with the
sample sizes:

```{r}
stddata
```


## Analysis

There are two steps to run the NMA/NMR model. The first step is to
create a JAGS model using `crossnma.model()` which produces the JAGS
code and the data. In the second step, the output of that function
will be used in `crossnma()` to run the analysis through JAGS.


### Unadjusted network meta-analysis

We start by providing the essential variables which - as stated
earlier - must have equal names in both data sets. Next, we give the
names of the datasets on participant-level (argument `prt.data`) and
aggregate data (argument `std.data`). By default, binary data is
analyzed using the odds ratio as a summary measure (`sm = "OR"`). The
`reference` treatment can be assigned which by default is the first
treatment (here: drug A). By default (`trt.effect = "random"`), we are
assigning a normal distribution to each relative treatment effect to
allow the synthesis across studies, see the table in Section
2.1. The different designs; RCT and NRS are combined with the
information taken at face-value as `method.bias = "naive"`.

Optionally, we can specify a prior to the common heterogeneity of the
treatment effect across studies. We indicate that distribution in the
argument `prior.tau.trt = "dunif(0, 3)"`.

Finally, we calculate the Surface Under the Cumulative Ranking (SUCRA)
in order to rank the treatments. It is essential to specify argument
`small.values` to get the correct ranking. For the RRMS studies, a
small number of relapses is desirable.


```{r}
# JAGS model: code + data
mod1 <- crossnma.model(treat, id, relapse, n, design,
  prt.data = ipddata, std.data = stddata,
  #---------- bias adjustment ----------
  method.bias = "naive",
  #---------- assign a prior ----------
  prior.tau.trt = "dunif(0, 3)",
  #---------- SUCRA ----------
  sucra = TRUE, small.values = "desirable"
  )
```

The network should be checked for its connectivity before running the
analysis. This is a vital step as the model will run even if the
network is not connected.

```{r, fig.width=4.5, fig.height=5,fig.show='hold',fig.align='center'}
netgraph(mod1, cex.points = n.trts, adj = 0.5, plastic = FALSE,
 number = TRUE, pos.number.of.studies = c(0.5, 0.4, 0.5, 0.5, 0.6, 0.5))
```

We are using argument `number.of.studies = TRUE` in order to print the number of studies in each direct comparison. The position of the number of studies is set by argument `pos.number.of.studies`.

Next, we fit the NMA model using `crossnma()`. We change the default settings for 
the number of iterations, burn-in and thinning.

```{r}
# Run JAGS
jagsfit1 <- crossnma(mod1, n.iter = 5000, n.burnin = 2000, thin = 1)
jagsfit1
```

By default (argument `backtransf = TRUE`), estimated odds ratios,
i.e., exp(d.B), exp(d.C) and exp(d.D), are printed. The value of tau
refers to the estimates of the heterogeneity standard deviation in the
relative treatment effects d.B, d.C and d.D across studies. The SUCRA
values are probabilities with treatment D being notably superior to
the other treatments.

We summarize the estimated parameters (argument `backtransf = FALSE`)
in the following table.

```{r, echo = FALSE}
knitr::kable(summary(jagsfit1, backtransf = FALSE), digits = 3)
```

We need also to assess the convergence of the MCMC chains either by
checking the Gelman and Rubin statistic, Rhat (it should be
approximately 1) in the table above or visually inspect the trace
plot.

```{r}
par(mar = rep(2, 4), mfrow = c(2, 3))
plot(jagsfit1)
```


### Unadjusted network meta-regression 

In this part, we set argument `cov1 = age` to run a NMR model with one
covariate. Again, datasets `ipddata` and `stddata` must use the same
variable name.

```{r}
# JAGS model: code + data
mod2 <- crossnma.model(treat, id, relapse, n, design,
  prt.data = ipddata, std.data = stddata,
  #---------- bias adjustment ----------
  method.bias = "naive",
  #----------  meta-regression ----------
  cov1 = age,
  split.regcoef = FALSE
  )
```

We could add two more covariates to the NMR model using arguments
`cov2` and `cov3`.

The MCMC is run under the same set up as in the network meta-analysis.

```{r}
# Run JAGS
jagsfit2 <- crossnma(mod2, n.iter = 5000, n.burnin = 2000, thin = 1)
```

and the output table is presented below

```{r, echo = FALSE}
knitr::kable(summary(jagsfit2, backtransf = FALSE), digits = 3)
```

Now, we additionally estimate b_1 which indicates the mean effect of
age and tau.b_1 which refers to the heterogeneity standard deviation
in the effect of age across studies. Here, we obtain a single estimate
because we choose to not split the within- and between-study age
coefficients $(\beta^w_{1,jbk} = \beta^B_{1,jbk}=\beta_{1,jbk})$ to
improve the convergence of MCMC.

The league table summarizes the relative effect with the 95% credible
interval of each treatment on the top compared to the treatment on the
left. All estimates are computed for participant age 38. We can
display the table in wide format

```{r, fig.width=6, fig.height=5,fig.show='hold',fig.align='center'}
league(jagsfit2, cov1.value = 38, digits = 2)
```

or in long format

```{r, fig.width=6, fig.height=5,fig.show='hold',fig.align='center'}
league(jagsfit2, cov1.value = 38, digits = 2, direction = "long")
```


### Using non-randomized studies (NRS) to construct priors for the treatment effects

To run NMA with a prior constructed from NRS, two additional arguments
are needed: we indicate using NRS as a prior by setting `
method.bias = "prior"`. That means that the model runs internally NMA
with only NRS data which are then used to construct informative
priors. This requires defining MCMC settings (the number of
adaptations, iterations, burn-ins, thinning and chains) in the
arguments starts with `run.nrs`.

In this method, the prior for the basic parameters is set to a normal
distribution. For basic parameters not examined in the NRS, the code
sets a minimally informative prior `d~dnorm(0, {15*ML}^2)`, where ML
is the largest maximum likelihood estimates of all relative treatment
effects in all studies. To account for possible bias, the means of the
distribution can be shifted by `run.nrs.mean.shift` and/or the
variance can be inflated by `run.nrs.var.infl` to control the
influence of NRS on the final estimation.

```{r}
# JAGS model: code + data
mod3 <- crossnma.model(treat, id, relapse, n, design,
  prt.data = ipddata, std.data = stddata,
  reference = "D",
  #----------  meta-regression ----------
  cov1 = age,
  split.regcoef = FALSE,
  #---------- bias adjustment ----------
  method.bias = "prior",
  run.nrs.trt.effect= "common",
  run.nrs.var.infl = 0.6, run.nrs.mean.shift = 0,
  run.nrs.n.iter = 10000, run.nrs.n.burnin = 4000,
  run.nrs.thin = 1, run.nrs.n.chains = 2
  )
```


```{r}
# Run JAGS
jagsfit3 <- crossnma(mod3, n.iter = 5000, n.burnin = 2000, thin = 1)
```

The heat plot summarizes the relative effect with the 95% credible
interval of each treatment on the top compared to the treatment on the
left. All estimates are computed for participant age 38.

```{r, fig.width=6, fig.height=5,fig.show='hold',fig.align='center'}
heatplot(jagsfit3, cov1.value = 38,
  size = 6, size.trt = 20, size.axis = 12)
```


### Bias-adjusted model 1

In this part, the overall relative treatment effects are estimated
from both NRS and RCT with adjustment to study-specific bias.

To fit the model, we set `method.bias = "adjust1"` and we need to
provide the bias variable `bias = rob` in the datasets. The direction
of bias is determined by the column `unfav = unfavored` which
indicates the unfavoured treatment. The mean bias effect can be
indicated by `bias.group`, $0$ (`bias.group = 0`), $g$ (`bias.group =
1`) or $g^{act}$ (`bias.group = 2`). By default, the effect of bias is
assumed to be additive `bias.type = "add"` and equal across studies
`bias.effect = "common"`. We also use the `year` of study publication
to estimate the study-probability of bias, `bias.covariate = year`.


```{r}
# JAGS model: code + data
mod4 <- crossnma.model(treat, id, relapse, n, design,
  prt.data = ipddata, std.data = stddata,
  #---------- bias adjustment ----------
  method.bias = "adjust1",
  bias.type = "add",
  bias.effect = "common",
  bias = rob,
  unfav = unfavored,
  bias.group = bias.group,
  bias.covariate = year
)
```

```{r}
# Run JAGS
jagsfit4 <- crossnma(mod4, n.iter = 5000, n.burnin = 2000, thin = 1)
```

The results are presented below

```{r, echo = FALSE}
knitr::kable(summary(jagsfit4, backtransf = FALSE), digits = 3)
```

The parameter `g` refers to the mean bias effect, common for all studies.


### Bias-adjusted model 2

The arguments for `method.bias = "adjust2"` are similar to the ones
used before in `method.bias = "adjust1"`.

```{r}
# JAGS model: code + data
mod5 <- crossnma.model(treat, id, relapse, n, design,
  prt.data = ipddata, std.data = stddata,
  #---------- bias adjustment ----------
  method.bias = "adjust2",
  bias.type = "add",
  bias = rob,
  unfav = unfavored,
  bias.group = bias.group
)
```


```{r}
# Run JAGS
jagsfit5 <- crossnma(mod5, n.iter = 5000, n.burnin = 2000, thin = 1)
```


```{r, echo = FALSE}
knitr::kable(summary(jagsfit5, backtransf = FALSE), digits = 3)
```

```{r, echo = FALSE, message = FALSE}
tools::compactPDF(path = ".", gs_quality = "ebook")
```

# References