# Example: Plaque psoriasis ML-NMR

    library(multinma)
    #> For execution on a local, multicore CPU with excess RAM we recommend calling
    #> options(mc.cores = parallel::detectCores())
    #>
    #> Attaching package: 'multinma'
    #> The following objects are masked from 'package:stats':
    #>
    #>     dgamma, pgamma, qgamma
    library(dplyr)      # dplyr and tidyr for data manipulation
    #>
    #> Attaching package: 'dplyr'
    #> The following objects are masked from 'package:stats':
    #>
    #>     filter, lag
    #> The following objects are masked from 'package:base':
    #>
    #>     intersect, setdiff, setequal, union
    library(tidyr)
    library(ggplot2)    # ggplot2 for plotting covariate distributions

    options(mc.cores = parallel::detectCores())

Simulated individual patient data (IPD) from trials comparing treatments for plaque psoriasis are found in the data set `plaque_psoriasis_ipd`. Aggregate data (AgD) are available on a further set of trials, found in the data set `plaque_psoriasis_agd`. In this vignette, we recreate the multilevel network meta-regression (ML-NMR) analyses performed by Phillippo et al. (2020; see also Phillippo 2019). We will analyse IPD from three studies, UNCOVER-1, UNCOVER-2, and UNCOVER-3 (Griffiths et al. 2015; Gordon et al. 2016), and AgD from one study, FIXTURE (Langley et al. 2014).

    pso_ipd <- filter(plaque_psoriasis_ipd,
                      studyc %in% c("UNCOVER-1", "UNCOVER-2", "UNCOVER-3"))

    pso_agd <- filter(plaque_psoriasis_agd,
                      studyc == "FIXTURE")

    head(pso_ipd)
    #>      studyc      trtc_long    trtc trtn pasi75 pasi90 pasi100 age  bmi pasi_w0  male bsa weight
    #> 1 UNCOVER-1 Ixekizumab Q2W IXE_Q2W    2      0      0       0  34 32.2    18.2  TRUE  18   98.1
    #> 2 UNCOVER-1 Ixekizumab Q2W IXE_Q2W    2      1      0       0  64 41.9    23.4  TRUE  33  129.6
    #> 3 UNCOVER-1 Ixekizumab Q2W IXE_Q2W    2      1      1       0  42 26.2    12.8  TRUE  33   78.0
    #> 4 UNCOVER-1 Ixekizumab Q2W IXE_Q2W    2      0      0       0  45 52.9    36.0 FALSE  50  139.9
    #> 5 UNCOVER-1 Ixekizumab Q2W IXE_Q2W    2      1      0       0  67 22.9    20.9 FALSE  35   54.2
    #> 6 UNCOVER-1 Ixekizumab Q2W IXE_Q2W    2      1      1       1  57 22.4    18.2  TRUE  29   67.5
    #>   durnpso prevsys   psa
    #> 1     6.7    TRUE  TRUE
    #> 2    14.5   FALSE  TRUE
    #> 3    26.5    TRUE FALSE
    #> 4    25.0    TRUE  TRUE
    #> 5    11.9    TRUE FALSE
    #> 6    15.2    TRUE FALSE
    head(pso_agd)
    #>    studyc          trtc_long    trtc trtn pasi75_r pasi75_n pasi90_r pasi90_n pasi100_r pasi100_n
    #> 1 FIXTURE         Etanercept     ETN    4      142      323       67      323        14       323
    #> 2 FIXTURE            Placebo     PBO    1       16      324        5      324         0       324
    #> 3 FIXTURE Secukinumab 150 mg SEC_150    5      219      327      137      327        47       327
    #> 4 FIXTURE Secukinumab 300 mg SEC_300    6      249      323      175      323        78       323
    #>   sample_size_w0 age_mean age_sd bmi_mean bmi_sd pasi_w0_mean pasi_w0_sd male bsa_mean bsa_sd
    #> 1            326     43.8   13.0     28.7    5.9         23.2        9.8 71.2     33.6   18.0
    #> 2            326     44.1   12.6     27.9    6.1         24.1       10.5 72.7     35.2   19.1
    #> 3            327     45.4   12.9     28.4    5.9         23.7       10.5 72.2     34.5   19.4
    #> 4            327     44.5   13.2     28.4    6.4         23.9        9.9 68.5     34.3   19.2
    #>   weight_mean weight_sd durnpso_mean durnpso_sd prevsys  psa
    #> 1        84.6      20.5         16.4       12.0    65.6 13.5
    #> 2        82.0      20.4         16.6       11.6    62.6 15.0
    #> 3        83.6      20.8         17.3       12.2    64.8 15.0
    #> 4        83.0      21.6         15.8       12.3    63.0 15.3

We consider running a ML-NMR adjusting for five potential effect-modifying covariates: duration of psoriasis `durnpso`, weight `weight`, previous systemic treatment `prevsys`, body surface area `bsa`, and psoriatic arthritis `psa`.

# Setup

## Preparing the data

We need to prepare the data so that it is in an acceptable format to run a ML-NMR model. Firstly, we need to handle the binary covariates `prevsys` and `psa`. In the IPD, these are coded as `TRUE` or `FALSE`, but in the AgD these are coded as percentages (out of 100). We need these to transform both of these sets of variables so that they are numeric and lie in the interval \([0,1]\), so that the variables are compatible across the data sources. Whilst we are here, we also transform body surface area `bsa` (a percentage) to lie in \([0,1]\), since that will make specifying an appropriate marginal distribution easier later, and rescale `weight` and `duration` to aid interpretation of the regression coefficients (in terms of 10 kilos and 10 years respectively). We also add in a `trtclass` variable, indicating which treatments belong to which classes. Finally, we check for missing values in the IPD.

    pso_ipd <- pso_ipd %>%
      mutate(# Variable transformations
             bsa = bsa / 100,
             prevsys = as.numeric(prevsys),
             psa = as.numeric(psa),
             weight = weight / 10,
             durnpso = durnpso / 10,
             # Treatment classes
             trtclass = case_when(trtn == 1 ~ "Placebo",
                                  trtn %in% c(2, 3, 5, 6) ~ "IL blocker",
                                  trtn == 4 ~ "TNFa blocker"),
             # Check complete cases for covariates of interest
             complete = complete.cases(durnpso, prevsys, bsa, weight, psa)
      )

    pso_agd <- pso_agd %>%
      mutate(
        # Variable transformations
        bsa_mean = bsa_mean / 100,
        bsa_sd = bsa_sd / 100,
        prevsys = prevsys / 100,
        psa = psa / 100,
        weight_mean = weight_mean / 10,
        weight_sd = weight_sd / 10,
        durnpso_mean = durnpso_mean / 10,
        durnpso_sd = durnpso_sd / 10,
        # Treatment classes
        trtclass = case_when(trtn == 1 ~ "Placebo",
                                  trtn %in% c(2, 3, 5, 6) ~ "IL blocker",
                                  trtn == 4 ~ "TNFa blocker")
      )

A small number of individuals have missing covariates:

    sum(!pso_ipd$complete)
    #> [1] 4
    mean(!pso_ipd$complete)
    #> [1] 0.001036807

Since the proportion of missing data is so small, we will simply exclude these individuals from the analysis.

    pso_ipd <- filter(pso_ipd, complete)

## Creating the network

Set up the network, setting the IPD with `set_ipd()`, AgD (arm-based) with `set_agd_arm()`, and combining together using `combine_network()`. We specify the binary `pasi75` outcome as `r` in the IPD, and the count outcome `pasi75_r` and denominator `pasi75_n` as `r` and `n` in the AgD. We specify the treatment classes with `trt_class = trtclass`.

    pso_net <- combine_network(
      set_ipd(pso_ipd,
              study = studyc,
              trt = trtc,
              r = pasi75,
              trt_class = trtclass),
      set_agd_arm(pso_agd,
                  study = studyc,
                  trt = trtc,
                  r = pasi75_r,
                  n = pasi75_n,
                  trt_class = trtclass)
    )

    pso_net
    #> A network with 3 IPD studies, and 1 AgD study (arm-based).
    #>
    #> ------------------------------------------------------------------- IPD studies ----
    #>  Study     Treatments
    #>  UNCOVER-1 3: IXE_Q2W | IXE_Q4W | PBO
    #>  UNCOVER-2 4: ETN | IXE_Q2W | IXE_Q4W | PBO
    #>  UNCOVER-3 4: ETN | IXE_Q2W | IXE_Q4W | PBO
    #>
    #>  Outcome type: binary
    #> ------------------------------------------------------- AgD studies (arm-based) ----
    #>  Study   Treatments
    #>  FIXTURE 4: ETN | PBO | SEC_150 | SEC_300
    #>
    #>  Outcome type: count
    #> ------------------------------------------------------------------------------------
    #> Total number of treatments: 6, in 3 classes
    #> Total number of studies: 4
    #> Reference treatment is: PBO
    #> Network is connected

We can produce a network plot with the `plot()` method:

    plot(pso_net, weight_nodes = TRUE, weight_edges = TRUE, show_trt_class = TRUE) +
      ggplot2::theme(legend.position = "bottom", legend.box = "vertical")

## Numerical integration for ML-NMR

ML-NMR models define the meta-regression model at the individual level, in exactly the same manner as a full-IPD meta-regression. ML-NMR then incorporates the AgD into the model by integrating this individual-level model over the covariate distribution in each AgD study (Phillippo et al. 2020; Phillippo 2019). Using integration, instead of simply “plugging-in” mean covariate values for the AgD studies, avoids aggregation bias when the link function is not the identity function.

This package utilises numerical integration to incorporate the aggregate data - specifically, quasi-Monte Carlo (QMC) integration with a Gaussian copula (Phillippo et al. 2020; Phillippo 2019). QMC integration is a very general and flexible integration approach, which typically requires far fewer integration points than standard (pseudo-random) Monte-Carlo integration to achieve the same numerical accuracy.1 A Gaussian copula allows us to account for correlations between covariates, which may have any specified marginal distributions.

We now set up the numerical integration for the network. The five covariates that we will consider adjusting for are body surface area `bsa`, duration of psoriasis `durnpso`, previous systemic treatment `prevsys`, psoriatic arthritis `psa`, and weight `weight`. We need to choose suitable marginal distributions for these covariates to draw the integration points from. `prevsys` and `psa` are binary covariates, so these are given a Bernoulli distribution. `bsa` is a percentage, so we choose a logit-Normal distribution (note, this requires the `logitnorm` package to be installed). We choose Gamma distributions for `durnpso` and `weight` to account for skewness. These choices seem to match well the marginal distributions observed in the IPD:

    # Get mean and sd of covariates in each study
    ipd_summary <- pso_ipd %>%
      group_by(studyc) %>%
      summarise_at(vars(weight, durnpso, bsa), list(mean = mean, sd = sd, min = min, max = max)) %>%
      pivot_longer(weight_mean:bsa_max, names_sep = "_", names_to = c("covariate", ".value")) %>%
      # Assign distributions
      mutate(dist = recode(covariate,
                           bsa = "dlogitnorm",
                           durnpso = "dgamma",
                           weight = "dgamma")) %>%
      # Compute density curves
      group_by(studyc, covariate) %>%
      mutate(value = if_else(dist == "dlogitnorm",
                             list(seq(0, 1, length.out = 101)),
                             list(seq(min*0.8, max*1.2, length.out = 101)))) %>%
      unnest(cols = value) %>%
      mutate(dens = eval(call(first(dist), x = value, mean = first(mean), sd = first(sd))))

    # Plot histograms and assumed densities
    pso_ipd %>%
      pivot_longer(c(weight, durnpso, bsa), names_to = "covariate", values_to = "value") %>%
    ggplot(aes(x = value)) +
      geom_histogram(aes(y = stat(density)),
                     binwidth = function(x) diff(range(x)) / nclass.Sturges(x),
                     boundary = 0,
                     fill = "grey50") +
      geom_line(aes(y = dens), data = ipd_summary,
                colour = "darkred", size = 0.5) +
      facet_wrap(~studyc + covariate, scales = "free", ncol = 3) +
      theme_multinma()

We add integration points to the AgD studies in the network using the `add_integration()` function. Marginal distributions for each covariate are specified using the `distr()` function, which takes a cumulative distribution function corresponding to the chosen marginal distribution, and arguments to that distribution as column names in the aggregate data. Since we do not know the correlations between covariates in the AgD studies, we impute these with the weighted mean of the correlations in the IPD studies (the default option).

    pso_net <- add_integration(pso_net,
      durnpso = distr(qgamma, mean = durnpso_mean, sd = durnpso_sd),
      prevsys = distr(qbern, prob = prevsys),
      bsa = distr(qlogitnorm, mean = bsa_mean, sd = bsa_sd),
      weight = distr(qgamma, mean = weight_mean, sd = weight_sd),
      psa = distr(qbern, prob = psa),
      n_int = 1000
    )
    #> Using weighted average correlation matrix computed from IPD studies.

> Note: This package provides several convenience functions for specifying these distributions, including `qgamma()` which allows for a parameterisation of the Gamma distribution in terms of mean and standard deviation, `qbern()` which provides the Bernoulli distribution, and `qlogitnorm()` which provides the logit-Normal distribution allowing for a parameterisation in terms of mean and standard deviation (requires the `logitnorm` package to be installed).

# ML-NMR models

We fit both fixed effect (FE) and random effects (RE) ML-NMR models.

## Fixed effect ML-NMR

First, we fit a FE ML-NMR model using the function `nma()`. Following (Phillippo et al. 2020) we specify weakly-informative \(N(0, 10^2)\) priors on each parameter. The range of parameter values implied by these prior distributions can be checked using the `summary()` method:

    summary(normal(scale = 10))
    #> A Normal prior distribution: location = 0, scale = 10.
    #> 50% of the prior density lies between -6.74 and 6.74.
    #> 95% of the prior density lies between -19.6 and 19.6.

The regression model is specified with `regression = ~(durnpso + prevsys + bsa + weight + psa)*.trt`, which will include the main (prognostic) effects of each covariate as well as interactions with treatment. We use a probit link function (`link = "probit"`), and specify that the two-parameter Binomial approximation for the aggregate-level likelihood should be used (`likelihood = "bernoulli2"`, where “bernoulli” refers to the individual-level likelihood, and “2” denotes the two-parameter adjustment to the aggregate-level likelihood) (Phillippo et al. 2020). We utilise the shared effect modifier assumption to help identify the model, setting treatment-covariate interactions to be equal within each class (`class_interactions = "common"`). We narrow the possible range for random initial values with `init_r = 0.1` (the default is `init_r = 2`), since probit models in particular are often hard to initialise. Using the QR decomposition (`QR = TRUE`) greatly improves sampling efficiency here, as is often the case for regression models.

    pso_fit_FE <- nma(pso_net,
                      trt_effects = "fixed",
                      link = "probit",
                      likelihood = "bernoulli2",
                      regression = ~(durnpso + prevsys + bsa + weight + psa)*.trt,
                      class_interactions = "common",
                      prior_intercept = normal(scale = 10),
                      prior_trt = normal(scale = 10),
                      prior_reg = normal(scale = 10),
                      init_r = 0.1,
                      QR = TRUE)
    #> Note: Setting "PBO" as the network reference treatment.

Basic parameter summaries are given by the `print()` method:

    print(pso_fit_FE)
    #> A fixed effects ML-NMR with a bernoulli2 likelihood (probit link).
    #> Regression model: ~(durnpso + prevsys + bsa + weight + psa) * .trt.
    #> Centred covariates at the following overall mean values:
    #>   durnpso   prevsys       bsa    weight       psa
    #> 1.8259772 0.6495432 0.2917678 8.9328063 0.2172823
    #> Inference for Stan model: binomial_2par.
    #> 4 chains, each with iter=2000; warmup=1000; thin=1;
    #> post-warmup draws per chain=1000, total post-warmup draws=4000.
    #>
    #>                                         mean se_mean   sd     2.5%      25%      50%      75%
    #> beta[durnpso]                           0.05    0.00 0.06    -0.08     0.00     0.04     0.09
    #> beta[prevsys]                          -0.13    0.00 0.16    -0.45    -0.24    -0.13    -0.03
    #> beta[bsa]                              -0.07    0.01 0.44    -0.97    -0.35    -0.05     0.23
    #> beta[weight]                            0.04    0.00 0.03    -0.02     0.02     0.04     0.06
    #> beta[psa]                              -0.08    0.00 0.17    -0.41    -0.19    -0.07     0.03
    #> beta[durnpso:.trtclassTNFa blocker]    -0.03    0.00 0.07    -0.18    -0.08    -0.03     0.02
    #> beta[durnpso:.trtclassIL blocker]      -0.01    0.00 0.07    -0.16    -0.06    -0.01     0.03
    #> beta[prevsys:.trtclassTNFa blocker]     0.19    0.00 0.19    -0.18     0.06     0.19     0.31
    #> beta[prevsys:.trtclassIL blocker]       0.06    0.00 0.18    -0.28    -0.06     0.06     0.18
    #> beta[bsa:.trtclassTNFa blocker]         0.05    0.01 0.52    -0.94    -0.30     0.05     0.39
    #> beta[bsa:.trtclassIL blocker]           0.29    0.01 0.48    -0.62    -0.05     0.29     0.61
    #> beta[weight:.trtclassTNFa blocker]     -0.17    0.00 0.04    -0.24    -0.19    -0.17    -0.14
    #> beta[weight:.trtclassIL blocker]       -0.10    0.00 0.03    -0.16    -0.12    -0.10    -0.08
    #> beta[psa:.trtclassTNFa blocker]        -0.06    0.00 0.21    -0.46    -0.20    -0.06     0.08
    #> beta[psa:.trtclassIL blocker]           0.00    0.00 0.18    -0.35    -0.12     0.00     0.13
    #> d[ETN]                                  1.55    0.00 0.08     1.40     1.50     1.55     1.61
    #> d[IXE_Q2W]                              2.95    0.00 0.09     2.78     2.90     2.95     3.01
    #> d[IXE_Q4W]                              2.54    0.00 0.08     2.38     2.49     2.54     2.60
    #> d[SEC_150]                              2.15    0.00 0.11     1.93     2.07     2.14     2.22
    #> d[SEC_300]                              2.45    0.00 0.12     2.23     2.37     2.45     2.53
    #> lp__                                -1576.42    0.09 3.47 -1583.90 -1578.63 -1576.08 -1573.92
    #>                                        97.5% n_eff Rhat
    #> beta[durnpso]                           0.17  6237    1
    #> beta[prevsys]                           0.17  5193    1
    #> beta[bsa]                               0.76  4328    1
    #> beta[weight]                            0.10  5050    1
    #> beta[psa]                               0.25  5605    1
    #> beta[durnpso:.trtclassTNFa blocker]     0.11  6338    1
    #> beta[durnpso:.trtclassIL blocker]       0.12  7413    1
    #> beta[prevsys:.trtclassTNFa blocker]     0.54  5881    1
    #> beta[prevsys:.trtclassIL blocker]       0.41  6725    1
    #> beta[bsa:.trtclassTNFa blocker]         1.09  5212    1
    #> beta[bsa:.trtclassIL blocker]           1.27  5211    1
    #> beta[weight:.trtclassTNFa blocker]     -0.10  5500    1
    #> beta[weight:.trtclassIL blocker]       -0.04  5701    1
    #> beta[psa:.trtclassTNFa blocker]         0.36  6190    1
    #> beta[psa:.trtclassIL blocker]           0.37  6477    1
    #> d[ETN]                                  1.71  4358    1
    #> d[IXE_Q2W]                              3.13  4845    1
    #> d[IXE_Q4W]                              2.71  4731    1
    #> d[SEC_150]                              2.37  5522    1
    #> d[SEC_300]                              2.69  6411    1
    #> lp__                                -1570.66  1598    1
    #>
    #> Samples were drawn using NUTS(diag_e) at Tue Mar 09 10:56:26 2021.
    #> For each parameter, n_eff is a crude measure of effective sample size,
    #> and Rhat is the potential scale reduction factor on split chains (at
    #> convergence, Rhat=1).

By default, summaries of the study-specific intercepts \(\mu_j\) are hidden, but could be examined by changing the `pars` argument:

    # Not run
    print(pso_fit_FE, pars = c("d", "beta", "mu"))

The prior and posterior distributions can be compared visually using the `plot_prior_posterior()` function:

    plot_prior_posterior(pso_fit_FE, prior = c("intercept", "trt", "reg"))

Plots of estimated numerical integration error are produced using the `plot_integration_error()` function:

    plot_integration_error(pso_fit_FE)

## Random effects ML-NMR

We now fit a RE model. Again, we specify weakly-informative \(N(0, 10^2)\) priors on each parameter, and now specify a \(\textrm{half-N}(0, 2.5^2)\) prior for the heterogeneity standard deviation \(\tau\). The range of parameter values implied by these prior distributions can be checked using the `summary()` method:

    summary(normal(scale = 10))
    #> A Normal prior distribution: location = 0, scale = 10.
    #> 50% of the prior density lies between -6.74 and 6.74.
    #> 95% of the prior density lies between -19.6 and 19.6.
    summary(half_normal(scale = 2.5))
    #> A half-Normal prior distribution: location = 0, scale = 2.5.
    #> 50% of the prior density lies between 0 and 1.69.
    #> 95% of the prior density lies between 0 and 4.9.

Fitting the model uses the same call to `nma()` as before, except now with `trt_effects = "random"`.

    pso_fit_RE <- nma(pso_net,
                      trt_effects = "random",
                      link = "probit",
                      likelihood = "bernoulli2",
                      regression = ~(durnpso + prevsys + bsa + weight + psa)*.trt,
                      class_interactions = "common",
                      prior_intercept = normal(scale = 10),
                      prior_trt = normal(scale = 10),
                      prior_reg = normal(scale = 10),
                      prior_het = half_normal(scale = 2.5),
                      init_r = 0.1,
                      QR = TRUE)
    #> Note: Setting "PBO" as the network reference treatment.
    #> Warning: There were 8 divergent transitions after warmup. See
    #> http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
    #> to find out why this is a problem and how to eliminate them.
    #> Warning: Examine the pairs() plot to diagnose sampling problems

Basic parameter summaries are given by the `print()` method:

    print(pso_fit_RE)
    #> A random effects ML-NMR with a bernoulli2 likelihood (probit link).
    #> Regression model: ~(durnpso + prevsys + bsa + weight + psa) * .trt.
    #> Centred covariates at the following overall mean values:
    #>   durnpso   prevsys       bsa    weight       psa
    #> 1.8259772 0.6495432 0.2917678 8.9328063 0.2172823
    #> Inference for Stan model: binomial_2par.
    #> 4 chains, each with iter=2000; warmup=1000; thin=1;
    #> post-warmup draws per chain=1000, total post-warmup draws=4000.
    #>
    #>                                         mean se_mean   sd     2.5%      25%      50%      75%
    #> beta[durnpso]                           0.05    0.00 0.06    -0.07     0.01     0.05     0.09
    #> beta[prevsys]                          -0.12    0.00 0.16    -0.44    -0.23    -0.12    -0.01
    #> beta[bsa]                              -0.09    0.01 0.45    -0.98    -0.39    -0.08     0.20
    #> beta[weight]                            0.04    0.00 0.03    -0.01     0.02     0.04     0.06
    #> beta[psa]                              -0.07    0.00 0.17    -0.40    -0.18    -0.06     0.05
    #> beta[durnpso:.trtclassTNFa blocker]    -0.04    0.00 0.07    -0.18    -0.09    -0.04     0.01
    #> beta[durnpso:.trtclassIL blocker]      -0.02    0.00 0.07    -0.15    -0.06    -0.02     0.03
    #> beta[prevsys:.trtclassTNFa blocker]     0.18    0.00 0.19    -0.19     0.05     0.18     0.31
    #> beta[prevsys:.trtclassIL blocker]       0.04    0.00 0.18    -0.30    -0.08     0.05     0.16
    #> beta[bsa:.trtclassTNFa blocker]         0.07    0.01 0.54    -0.98    -0.29     0.07     0.42
    #> beta[bsa:.trtclassIL blocker]           0.33    0.01 0.49    -0.64     0.00     0.34     0.67
    #> beta[weight:.trtclassTNFa blocker]     -0.17    0.00 0.04    -0.24    -0.20    -0.17    -0.15
    #> beta[weight:.trtclassIL blocker]       -0.10    0.00 0.03    -0.17    -0.12    -0.10    -0.08
    #> beta[psa:.trtclassTNFa blocker]        -0.07    0.00 0.20    -0.46    -0.21    -0.07     0.06
    #> beta[psa:.trtclassIL blocker]          -0.01    0.00 0.19    -0.36    -0.14    -0.01     0.11
    #> d[ETN]                                  1.56    0.00 0.14     1.28     1.47     1.56     1.64
    #> d[IXE_Q2W]                              2.98    0.00 0.15     2.68     2.89     2.97     3.06
    #> d[IXE_Q4W]                              2.57    0.00 0.15     2.28     2.48     2.56     2.65
    #> d[SEC_150]                              2.12    0.01 0.24     1.63     2.00     2.12     2.26
    #> d[SEC_300]                              2.43    0.01 0.24     1.91     2.30     2.43     2.56
    #> lp__                                -1580.32    0.15 4.83 -1591.00 -1583.40 -1580.11 -1576.83
    #> tau                                     0.19    0.01 0.12     0.02     0.10     0.17     0.25
    #>                                        97.5% n_eff Rhat
    #> beta[durnpso]                           0.17  4902 1.00
    #> beta[prevsys]                           0.20  4469 1.00
    #> beta[bsa]                               0.77  4351 1.00
    #> beta[weight]                            0.10  4913 1.00
    #> beta[psa]                               0.26  5172 1.00
    #> beta[durnpso:.trtclassTNFa blocker]     0.11  5382 1.00
    #> beta[durnpso:.trtclassIL blocker]       0.12  5849 1.00
    #> beta[prevsys:.trtclassTNFa blocker]     0.55  4955 1.00
    #> beta[prevsys:.trtclassIL blocker]       0.39  4949 1.00
    #> beta[bsa:.trtclassTNFa blocker]         1.15  4568 1.00
    #> beta[bsa:.trtclassIL blocker]           1.27  5034 1.00
    #> beta[weight:.trtclassTNFa blocker]     -0.10  4953 1.00
    #> beta[weight:.trtclassIL blocker]       -0.04  5626 1.00
    #> beta[psa:.trtclassTNFa blocker]         0.35  5444 1.00
    #> beta[psa:.trtclassIL blocker]           0.36  5650 1.00
    #> d[ETN]                                  1.86  2160 1.00
    #> d[IXE_Q2W]                              3.28  1825 1.00
    #> d[IXE_Q4W]                              2.88  2190 1.00
    #> d[SEC_150]                              2.60  1515 1.00
    #> d[SEC_300]                              2.88  1643 1.00
    #> lp__                                -1571.72  1009 1.01
    #> tau                                     0.49   433 1.01
    #>
    #> Samples were drawn using NUTS(diag_e) at Tue Mar 09 11:04:31 2021.
    #> For each parameter, n_eff is a crude measure of effective sample size,
    #> and Rhat is the potential scale reduction factor on split chains (at
    #> convergence, Rhat=1).

By default, summaries of the study-specific intercepts \(\mu_j\) and study-specific relative effects \(\delta_{jk}\) are hidden, but could be examined by changing the `pars` argument:

    # Not run
    print(pso_fit_RE, pars = c("d", "beta", "tau", "mu", "delta"))

There are a number of divergent transitions, which we can investigate using the `pairs()` method:

    pairs(pso_fit_RE, pars = c("delta[UNCOVER-2: ETN]", "d[ETN]", "tau", "lp__"))

The divergent transition errors (red crosses) seem to be concentrated in the upper tail of the heterogeneity standard deviation parameter. This suggests that the information to identify the heterogeneity parameter is weak - we have only four studies in the network - and that a more informative prior distribution might aid estimation.

The prior and posterior distributions can be compared visually using the `plot_prior_posterior()` function:

    plot_prior_posterior(pso_fit_RE, prior = c("intercept", "trt", "reg", "het"))

Plots of estimated numerical integration error are produced using the `plot_integration_error()` function:

    plot_integration_error(pso_fit_RE)

# Model comparison

The model fit under the FE and RE models can be checked using the `dic()` function.

    (pso_dic_FE <- dic(pso_fit_FE))
    #> Residual deviance: 3129.6 (on 3858 data points)
    #>                pD: 24.3
    #>               DIC: 3153.9
    (pso_dic_RE <- dic(pso_fit_RE))
    #> Residual deviance: 3123.7 (on 3858 data points)
    #>                pD: 28.3
    #>               DIC: 3152

The DIC is similar between the FE and RE models, suggesting that there is little evidence for any residual heterogeneity.

# Further results

Parameter estimates can be plotted using the `plot()` method, for example to examine the estimated regression coefficients:

    plot(pso_fit_FE,
         pars = "beta",
         stat = "halfeye",
         ref_line = 0)

Plots of posterior summaries are based on the `ggdist` package, which allows a great degree of flexibility, and can be further customised using `ggplot2` commands. In the above command we specify a `"halfeye"` plot, which shows the posterior density along with posterior medians (points) and 95% Credible Intervals (thin line) with 66% inner bands (thicker line) by default. For more details on the plotting options see `?plot.nma_summary`.

We can produce population-adjusted relative effects for each study population in the network using the `relative_effects()` function.

    (pso_releff_FE <- relative_effects(pso_fit_FE))
    #> ---------------------------------------------------------------- Study: FIXTURE ----
    #>
    #> Covariate values:
    #>  durnpso prevsys  bsa weight  psa
    #>     1.65    0.64 0.34   8.32 0.15
    #>
    #>                     mean   sd 2.5%  25%  50%  75% 97.5% Bulk_ESS Tail_ESS Rhat
    #> d[FIXTURE: ETN]     1.66 0.09 1.49 1.60 1.66 1.73  1.84     3665     3399    1
    #> d[FIXTURE: IXE_Q2W] 3.03 0.10 2.83 2.96 3.03 3.10  3.24     4426     2875    1
    #> d[FIXTURE: IXE_Q4W] 2.62 0.10 2.43 2.55 2.62 2.68  2.81     4454     2538    1
    #> d[FIXTURE: SEC_150] 2.22 0.12 2.00 2.14 2.22 2.30  2.45     4732     3065    1
    #> d[FIXTURE: SEC_300] 2.53 0.12 2.30 2.45 2.52 2.61  2.77     5609     3106    1
    #>
    #> -------------------------------------------------------------- Study: UNCOVER-1 ----
    #>
    #> Covariate values:
    #>  durnpso prevsys  bsa weight  psa
    #>        2    0.73 0.28   9.24 0.28
    #>
    #>                       mean   sd 2.5%  25%  50%  75% 97.5% Bulk_ESS Tail_ESS Rhat
    #> d[UNCOVER-1: ETN]     1.51 0.08 1.34 1.45 1.51 1.56  1.67     5016     3477    1
    #> d[UNCOVER-1: IXE_Q2W] 2.92 0.09 2.75 2.86 2.92 2.98  3.10     5307     2714    1
    #> d[UNCOVER-1: IXE_Q4W] 2.51 0.08 2.35 2.46 2.51 2.56  2.67     5168     2966    1
    #> d[UNCOVER-1: SEC_150] 2.11 0.12 1.89 2.03 2.11 2.19  2.34     5955     3539    1
    #> d[UNCOVER-1: SEC_300] 2.42 0.12 2.20 2.34 2.42 2.50  2.67     6938     3358    1
    #>
    #> -------------------------------------------------------------- Study: UNCOVER-2 ----
    #>
    #> Covariate values:
    #>  durnpso prevsys  bsa weight  psa
    #>     1.87    0.64 0.27   9.17 0.24
    #>
    #>                       mean   sd 2.5%  25%  50%  75% 97.5% Bulk_ESS Tail_ESS Rhat
    #> d[UNCOVER-2: ETN]     1.51 0.08 1.35 1.45 1.51 1.57  1.67     4959     3392    1
    #> d[UNCOVER-2: IXE_Q2W] 2.92 0.09 2.75 2.86 2.92 2.98  3.09     5257     2738    1
    #> d[UNCOVER-2: IXE_Q4W] 2.51 0.08 2.35 2.46 2.51 2.56  2.67     5116     3013    1
    #> d[UNCOVER-2: SEC_150] 2.11 0.11 1.89 2.04 2.11 2.19  2.34     5896     3544    1
    #> d[UNCOVER-2: SEC_300] 2.42 0.12 2.20 2.34 2.42 2.50  2.67     6856     3355    1
    #>
    #> -------------------------------------------------------------- Study: UNCOVER-3 ----
    #>
    #> Covariate values:
    #>  durnpso prevsys  bsa weight psa
    #>     1.78    0.59 0.28   9.01 0.2
    #>
    #>                       mean   sd 2.5%  25%  50%  75% 97.5% Bulk_ESS Tail_ESS Rhat
    #> d[UNCOVER-3: ETN]     1.53 0.08 1.37 1.47 1.53 1.59  1.69     4613     3242    1
    #> d[UNCOVER-3: IXE_Q2W] 2.94 0.09 2.77 2.88 2.94 3.00  3.11     5064     3071    1
    #> d[UNCOVER-3: IXE_Q4W] 2.53 0.08 2.37 2.47 2.53 2.58  2.69     4929     2854    1
    #> d[UNCOVER-3: SEC_150] 2.13 0.11 1.91 2.05 2.13 2.21  2.35     5704     3425    1
    #> d[UNCOVER-3: SEC_300] 2.44 0.12 2.22 2.36 2.44 2.52  2.68     6706     3355    1
    plot(pso_releff_FE, ref_line = 0)

Predicted probabilities of achieving PASI 75 in each study population on each treatment are produced using the `predict()` method. The argument `type = "reponse"` specifies that we want predicted probabilities, rather than probit probabilities.

    (pso_pred_FE <- predict(pso_fit_FE, type = "response"))
    #> ---------------------------------------------------------------- Study: FIXTURE ----
    #>
    #>                        mean   sd 2.5%  25%  50%  75% 97.5% Bulk_ESS Tail_ESS Rhat
    #> pred[FIXTURE: PBO]     0.04 0.01 0.03 0.04 0.04 0.05  0.06     4301     3194    1
    #> pred[FIXTURE: ETN]     0.46 0.02 0.41 0.44 0.46 0.47  0.50     6261     3090    1
    #> pred[FIXTURE: IXE_Q2W] 0.89 0.02 0.85 0.88 0.89 0.90  0.92     6659     2965    1
    #> pred[FIXTURE: IXE_Q4W] 0.80 0.03 0.74 0.78 0.80 0.81  0.84     7048     3083    1
    #> pred[FIXTURE: SEC_150] 0.67 0.03 0.62 0.65 0.67 0.69  0.72     9087     3123    1
    #> pred[FIXTURE: SEC_300] 0.77 0.02 0.72 0.76 0.77 0.79  0.82     9366     3197    1
    #>
    #> -------------------------------------------------------------- Study: UNCOVER-1 ----
    #>
    #>                          mean   sd 2.5%  25%  50%  75% 97.5% Bulk_ESS Tail_ESS Rhat
    #> pred[UNCOVER-1: PBO]     0.06 0.01 0.04 0.05 0.06 0.06  0.08     5731     2394    1
    #> pred[UNCOVER-1: ETN]     0.46 0.03 0.41 0.44 0.46 0.48  0.52     8706     3029    1
    #> pred[UNCOVER-1: IXE_Q2W] 0.90 0.01 0.88 0.89 0.90 0.91  0.92     9824     3034    1
    #> pred[UNCOVER-1: IXE_Q4W] 0.81 0.02 0.78 0.80 0.81 0.82  0.84     9083     3313    1
    #> pred[UNCOVER-1: SEC_150] 0.69 0.04 0.60 0.66 0.69 0.72  0.76     6937     3364    1
    #> pred[UNCOVER-1: SEC_300] 0.78 0.04 0.71 0.76 0.79 0.81  0.85     7919     3079    1
    #>
    #> -------------------------------------------------------------- Study: UNCOVER-2 ----
    #>
    #>                          mean   sd 2.5%  25%  50%  75% 97.5% Bulk_ESS Tail_ESS Rhat
    #> pred[UNCOVER-2: PBO]     0.05 0.01 0.03 0.04 0.05 0.05  0.06     5692     3167    1
    #> pred[UNCOVER-2: ETN]     0.42 0.02 0.38 0.41 0.42 0.43  0.46     8123     3337    1
    #> pred[UNCOVER-2: IXE_Q2W] 0.88 0.01 0.86 0.87 0.88 0.89  0.91     8887     3089    1
    #> pred[UNCOVER-2: IXE_Q4W] 0.78 0.02 0.75 0.77 0.78 0.79  0.81     9978     3434    1
    #> pred[UNCOVER-2: SEC_150] 0.65 0.04 0.57 0.62 0.65 0.68  0.72     7118     3464    1
    #> pred[UNCOVER-2: SEC_300] 0.75 0.04 0.68 0.73 0.75 0.78  0.82     7371     3472    1
    #>
    #> -------------------------------------------------------------- Study: UNCOVER-3 ----
    #>
    #>                          mean   sd 2.5%  25%  50%  75% 97.5% Bulk_ESS Tail_ESS Rhat
    #> pred[UNCOVER-3: PBO]     0.08 0.01 0.06 0.07 0.08 0.08  0.10     5211     3361    1
    #> pred[UNCOVER-3: ETN]     0.53 0.02 0.49 0.52 0.53 0.54  0.57     9464     3018    1
    #> pred[UNCOVER-3: IXE_Q2W] 0.93 0.01 0.91 0.92 0.93 0.93  0.94     9399     2992    1
    #> pred[UNCOVER-3: IXE_Q4W] 0.85 0.01 0.83 0.85 0.85 0.86  0.88     8633     2981    1
    #> pred[UNCOVER-3: SEC_150] 0.75 0.03 0.68 0.72 0.75 0.77  0.81     6904     3423    1
    #> pred[UNCOVER-3: SEC_300] 0.83 0.03 0.77 0.81 0.83 0.85  0.88     7580     3359    1
    plot(pso_pred_FE, ref_line = c(0, 1))

We can produce population-adjusted ranks, rank probabilities, and cumulative rank probabilities in each study population using the `posterior_ranks()` and `posterior_rank_probs()` functions (although here the ranks are unchanged between populations, as the distributions of effect modifiers are similar). We specify `lower_better = FALSE`, since a higher outcome is better (higher chance of achieving PASI 75).

    (pso_ranks_FE <- posterior_ranks(pso_fit_FE, lower_better = FALSE))
    #> ---------------------------------------------------------------- Study: FIXTURE ----
    #>
    #> Covariate values:
    #>  durnpso prevsys  bsa weight  psa
    #>     1.65    0.64 0.34   8.32 0.15
    #>
    #>                        mean   sd 2.5% 25% 50% 75% 97.5% Bulk_ESS Tail_ESS Rhat
    #> rank[FIXTURE: PBO]     6.00 0.00    6   6   6   6     6       NA       NA   NA
    #> rank[FIXTURE: ETN]     5.00 0.00    5   5   5   5     5       NA       NA   NA
    #> rank[FIXTURE: IXE_Q2W] 1.00 0.02    1   1   1   1     1     4016     4016    1
    #> rank[FIXTURE: IXE_Q4W] 2.22 0.42    2   2   2   2     3     4083       NA    1
    #> rank[FIXTURE: SEC_150] 4.00 0.04    4   4   4   4     4     3160       NA    1
    #> rank[FIXTURE: SEC_300] 2.78 0.42    2   3   3   3     3     4203     3160    1
    #>
    #> -------------------------------------------------------------- Study: UNCOVER-1 ----
    #>
    #> Covariate values:
    #>  durnpso prevsys  bsa weight  psa
    #>        2    0.73 0.28   9.24 0.28
    #>
    #>                          mean   sd 2.5% 25% 50% 75% 97.5% Bulk_ESS Tail_ESS Rhat
    #> rank[UNCOVER-1: PBO]     6.00 0.00    6   6   6   6     6       NA       NA   NA
    #> rank[UNCOVER-1: ETN]     5.00 0.00    5   5   5   5     5       NA       NA   NA
    #> rank[UNCOVER-1: IXE_Q2W] 1.00 0.02    1   1   1   1     1     4016     4016    1
    #> rank[UNCOVER-1: IXE_Q4W] 2.22 0.42    2   2   2   2     3     4083       NA    1
    #> rank[UNCOVER-1: SEC_150] 4.00 0.04    4   4   4   4     4     3160       NA    1
    #> rank[UNCOVER-1: SEC_300] 2.78 0.42    2   3   3   3     3     4203     3160    1
    #>
    #> -------------------------------------------------------------- Study: UNCOVER-2 ----
    #>
    #> Covariate values:
    #>  durnpso prevsys  bsa weight  psa
    #>     1.87    0.64 0.27   9.17 0.24
    #>
    #>                          mean   sd 2.5% 25% 50% 75% 97.5% Bulk_ESS Tail_ESS Rhat
    #> rank[UNCOVER-2: PBO]     6.00 0.00    6   6   6   6     6       NA       NA   NA
    #> rank[UNCOVER-2: ETN]     5.00 0.00    5   5   5   5     5       NA       NA   NA
    #> rank[UNCOVER-2: IXE_Q2W] 1.00 0.02    1   1   1   1     1     4016     4016    1
    #> rank[UNCOVER-2: IXE_Q4W] 2.22 0.42    2   2   2   2     3     4083       NA    1
    #> rank[UNCOVER-2: SEC_150] 4.00 0.04    4   4   4   4     4     3160       NA    1
    #> rank[UNCOVER-2: SEC_300] 2.78 0.42    2   3   3   3     3     4203     3160    1
    #>
    #> -------------------------------------------------------------- Study: UNCOVER-3 ----
    #>
    #> Covariate values:
    #>  durnpso prevsys  bsa weight psa
    #>     1.78    0.59 0.28   9.01 0.2
    #>
    #>                          mean   sd 2.5% 25% 50% 75% 97.5% Bulk_ESS Tail_ESS Rhat
    #> rank[UNCOVER-3: PBO]     6.00 0.00    6   6   6   6     6       NA       NA   NA
    #> rank[UNCOVER-3: ETN]     5.00 0.00    5   5   5   5     5       NA       NA   NA
    #> rank[UNCOVER-3: IXE_Q2W] 1.00 0.02    1   1   1   1     1     4016     4016    1
    #> rank[UNCOVER-3: IXE_Q4W] 2.22 0.42    2   2   2   2     3     4083       NA    1
    #> rank[UNCOVER-3: SEC_150] 4.00 0.04    4   4   4   4     4     3160       NA    1
    #> rank[UNCOVER-3: SEC_300] 2.78 0.42    2   3   3   3     3     4203     3160    1
    plot(pso_ranks_FE)

    (pso_rankprobs_FE <- posterior_rank_probs(pso_fit_FE, lower_better = FALSE))
    #> ---------------------------------------------------------------- Study: FIXTURE ----
    #>
    #> Covariate values:
    #>  durnpso prevsys  bsa weight  psa
    #>     1.65    0.64 0.34   8.32 0.15
    #>
    #>                     p_rank[1] p_rank[2] p_rank[3] p_rank[4] p_rank[5] p_rank[6]
    #> d[FIXTURE: PBO]             0      0.00      0.00         0         0         1
    #> d[FIXTURE: ETN]             0      0.00      0.00         0         1         0
    #> d[FIXTURE: IXE_Q2W]         1      0.00      0.00         0         0         0
    #> d[FIXTURE: IXE_Q4W]         0      0.78      0.22         0         0         0
    #> d[FIXTURE: SEC_150]         0      0.00      0.00         1         0         0
    #> d[FIXTURE: SEC_300]         0      0.22      0.77         0         0         0
    #>
    #> -------------------------------------------------------------- Study: UNCOVER-1 ----
    #>
    #> Covariate values:
    #>  durnpso prevsys  bsa weight  psa
    #>        2    0.73 0.28   9.24 0.28
    #>
    #>                       p_rank[1] p_rank[2] p_rank[3] p_rank[4] p_rank[5] p_rank[6]
    #> d[UNCOVER-1: PBO]             0      0.00      0.00         0         0         1
    #> d[UNCOVER-1: ETN]             0      0.00      0.00         0         1         0
    #> d[UNCOVER-1: IXE_Q2W]         1      0.00      0.00         0         0         0
    #> d[UNCOVER-1: IXE_Q4W]         0      0.78      0.22         0         0         0
    #> d[UNCOVER-1: SEC_150]         0      0.00      0.00         1         0         0
    #> d[UNCOVER-1: SEC_300]         0      0.22      0.77         0         0         0
    #>
    #> -------------------------------------------------------------- Study: UNCOVER-2 ----
    #>
    #> Covariate values:
    #>  durnpso prevsys  bsa weight  psa
    #>     1.87    0.64 0.27   9.17 0.24
    #>
    #>                       p_rank[1] p_rank[2] p_rank[3] p_rank[4] p_rank[5] p_rank[6]
    #> d[UNCOVER-2: PBO]             0      0.00      0.00         0         0         1
    #> d[UNCOVER-2: ETN]             0      0.00      0.00         0         1         0
    #> d[UNCOVER-2: IXE_Q2W]         1      0.00      0.00         0         0         0
    #> d[UNCOVER-2: IXE_Q4W]         0      0.78      0.22         0         0         0
    #> d[UNCOVER-2: SEC_150]         0      0.00      0.00         1         0         0
    #> d[UNCOVER-2: SEC_300]         0      0.22      0.77         0         0         0
    #>
    #> -------------------------------------------------------------- Study: UNCOVER-3 ----
    #>
    #> Covariate values:
    #>  durnpso prevsys  bsa weight psa
    #>     1.78    0.59 0.28   9.01 0.2
    #>
    #>                       p_rank[1] p_rank[2] p_rank[3] p_rank[4] p_rank[5] p_rank[6]
    #> d[UNCOVER-3: PBO]             0      0.00      0.00         0         0         1
    #> d[UNCOVER-3: ETN]             0      0.00      0.00         0         1         0
    #> d[UNCOVER-3: IXE_Q2W]         1      0.00      0.00         0         0         0
    #> d[UNCOVER-3: IXE_Q4W]         0      0.78      0.22         0         0         0
    #> d[UNCOVER-3: SEC_150]         0      0.00      0.00         1         0         0
    #> d[UNCOVER-3: SEC_300]         0      0.22      0.77         0         0         0
    plot(pso_rankprobs_FE)

    (pso_cumrankprobs_FE <- posterior_rank_probs(pso_fit_FE, lower_better = FALSE, cumulative = TRUE))
    #> ---------------------------------------------------------------- Study: FIXTURE ----
    #>
    #> Covariate values:
    #>  durnpso prevsys  bsa weight  psa
    #>     1.65    0.64 0.34   8.32 0.15
    #>
    #>                     p_rank[1] p_rank[2] p_rank[3] p_rank[4] p_rank[5] p_rank[6]
    #> d[FIXTURE: PBO]             0      0.00         0         0         0         1
    #> d[FIXTURE: ETN]             0      0.00         0         0         1         1
    #> d[FIXTURE: IXE_Q2W]         1      1.00         1         1         1         1
    #> d[FIXTURE: IXE_Q4W]         0      0.78         1         1         1         1
    #> d[FIXTURE: SEC_150]         0      0.00         0         1         1         1
    #> d[FIXTURE: SEC_300]         0      0.22         1         1         1         1
    #>
    #> -------------------------------------------------------------- Study: UNCOVER-1 ----
    #>
    #> Covariate values:
    #>  durnpso prevsys  bsa weight  psa
    #>        2    0.73 0.28   9.24 0.28
    #>
    #>                       p_rank[1] p_rank[2] p_rank[3] p_rank[4] p_rank[5] p_rank[6]
    #> d[UNCOVER-1: PBO]             0      0.00         0         0         0         1
    #> d[UNCOVER-1: ETN]             0      0.00         0         0         1         1
    #> d[UNCOVER-1: IXE_Q2W]         1      1.00         1         1         1         1
    #> d[UNCOVER-1: IXE_Q4W]         0      0.78         1         1         1         1
    #> d[UNCOVER-1: SEC_150]         0      0.00         0         1         1         1
    #> d[UNCOVER-1: SEC_300]         0      0.22         1         1         1         1
    #>
    #> -------------------------------------------------------------- Study: UNCOVER-2 ----
    #>
    #> Covariate values:
    #>  durnpso prevsys  bsa weight  psa
    #>     1.87    0.64 0.27   9.17 0.24
    #>
    #>                       p_rank[1] p_rank[2] p_rank[3] p_rank[4] p_rank[5] p_rank[6]
    #> d[UNCOVER-2: PBO]             0      0.00         0         0         0         1
    #> d[UNCOVER-2: ETN]             0      0.00         0         0         1         1
    #> d[UNCOVER-2: IXE_Q2W]         1      1.00         1         1         1         1
    #> d[UNCOVER-2: IXE_Q4W]         0      0.78         1         1         1         1
    #> d[UNCOVER-2: SEC_150]         0      0.00         0         1         1         1
    #> d[UNCOVER-2: SEC_300]         0      0.22         1         1         1         1
    #>
    #> -------------------------------------------------------------- Study: UNCOVER-3 ----
    #>
    #> Covariate values:
    #>  durnpso prevsys  bsa weight psa
    #>     1.78    0.59 0.28   9.01 0.2
    #>
    #>                       p_rank[1] p_rank[2] p_rank[3] p_rank[4] p_rank[5] p_rank[6]
    #> d[UNCOVER-3: PBO]             0      0.00         0         0         0         1
    #> d[UNCOVER-3: ETN]             0      0.00         0         0         1         1
    #> d[UNCOVER-3: IXE_Q2W]         1      1.00         1         1         1         1
    #> d[UNCOVER-3: IXE_Q4W]         0      0.78         1         1         1         1
    #> d[UNCOVER-3: SEC_150]         0      0.00         0         1         1         1
    #> d[UNCOVER-3: SEC_300]         0      0.22         1         1         1         1
    plot(pso_cumrankprobs_FE)

All of the above estimates (relative effects, predictions, rankings) can also be produced for a specific target population or populations by providing a suitable `newdata` argument to for function (and a `baseline` distribution for `predict()`).

To produce population-adjusted relative effects (and corresponding rankings) for a chosen target population, we require only the mean covariate values in that population. For example, `newdata` could provide the following mean covariate values:

    new_agd_means <- tibble(
      bsa = 0.6,
      prevsys = 0.1,
      psa = 0.2,
      weight = 10,
      durnpso = 3)

Population-adjusted relative effects in this target population are then calculated using the `relative_effects()` function, and can be plotted with the corresponding `plot()` method:

    (pso_releff_FE_new <- relative_effects(pso_fit_FE, newdata = new_agd_means))
    #> ------------------------------------------------------------------ Study: New 1 ----
    #>
    #> Covariate values:
    #>  durnpso prevsys bsa weight psa
    #>        3     0.1 0.6     10 0.2
    #>
    #>                   mean   sd 2.5%  25%  50%  75% 97.5% Bulk_ESS Tail_ESS Rhat
    #> d[New 1: ETN]     1.25 0.23 0.80 1.10 1.25 1.40  1.71     5844     3065    1
    #> d[New 1: IXE_Q2W] 2.89 0.22 2.47 2.74 2.88 3.04  3.33     6062     3114    1
    #> d[New 1: IXE_Q4W] 2.48 0.22 2.06 2.32 2.47 2.62  2.92     6043     2967    1
    #> d[New 1: SEC_150] 2.08 0.23 1.65 1.93 2.08 2.23  2.53     6887     3562    1
    #> d[New 1: SEC_300] 2.39 0.23 1.96 2.23 2.38 2.54  2.85     7139     3326    1
    plot(pso_releff_FE_new, ref_line = 0)

For absolute predictions, we require information about the full covariate distribution in the target population, not just the mean values. If IPD are available for the target population, `newdata` is simply a data frame of the IPD. If AgD are available for the target population, `newdata` must be a data frame with added integration points created using the `add_integration()` function.

For example, suppose the aggregate target population introduced above had the following covariate means and standard deviations (for continuous covariates) or proportions (for discrete covariates):

    new_agd_int <- tibble(
      bsa_mean = 0.6,
      bsa_sd = 0.3,
      prevsys = 0.1,
      psa = 0.2,
      weight_mean = 10,
      weight_sd = 1,
      durnpso_mean = 3,
      durnpso_sd = 1
    )

We add integration points to this data frame in a similar manner to before. Again, we need to supply a correlation matrix for the joint covariate distribution; we use the same weighted mean correlation matrix computed earlier from the IPD in the network, which is stored in the network object as `int_cor`.

    new_agd_int <- add_integration(new_agd_int,
      durnpso = distr(qgamma, mean = durnpso_mean, sd = durnpso_sd),
      prevsys = distr(qbern, prob = prevsys),
      bsa = distr(qlogitnorm, mean = bsa_mean, sd = bsa_sd),
      weight = distr(qgamma, mean = weight_mean, sd = weight_sd),
      psa = distr(qbern, prob = psa),
      cor = pso_net$int_cor,
      n_int = 1000)

Predicted probabilities of achieving PASI 75 in this target population, given a \(N(-1.75, 0.08^2)\) distribution on the baseline probit-probability of response on Placebo (at the reference levels of the covariates), are then produced using the `predict()` method:

    (pso_pred_FE_new <- predict(pso_fit_FE,
                                type = "response",
                                newdata = new_agd_int,
                                baseline = distr(qnorm, -1.75, 0.08)))
    #> ------------------------------------------------------------------ Study: New 1 ----
    #>
    #>                      mean   sd 2.5%  25%  50%  75% 97.5% Bulk_ESS Tail_ESS Rhat
    #> pred[New 1: PBO]     0.06 0.03 0.02 0.04 0.06 0.07  0.12     5026     3172    1
    #> pred[New 1: ETN]     0.37 0.06 0.26 0.33 0.37 0.41  0.49     5967     3444    1
    #> pred[New 1: IXE_Q2W] 0.90 0.03 0.84 0.88 0.90 0.92  0.94     5185     3657    1
    #> pred[New 1: IXE_Q4W] 0.81 0.04 0.72 0.78 0.81 0.83  0.87     5250     3689    1
    #> pred[New 1: SEC_150] 0.68 0.05 0.57 0.65 0.68 0.72  0.78     5367     3540    1
    #> pred[New 1: SEC_300] 0.78 0.05 0.68 0.75 0.78 0.81  0.86     5880     3578    1
    plot(pso_pred_FE_new, ref_line = c(0, 1))

# References

Caflisch, R. E. 1998. “Monte Carlo and Quasi-Monte Carlo Methods.” Acta Numerica 7: 1–49. https://doi.org/10.1017/S0962492900002804.
Gordon, K. B., A. Blauvelt, K. A. Papp, R. G. Langley, T. Luger, M. Ohtsuki, K. Reich, et al. 2016. “Phase 3 Trials of Ixekizumab in Moderate-to-Severe Plaque Psoriasis.” New England Journal of Medicine 375 (4): 345–56. https://doi.org/10.1056/nejmoa1512711.
Griffiths, C. E. M., K. Reich, M. Lebwohl, P. van de Kerkhof, C. Paul, A. Menter, G. S. Cameron, et al. 2015. “Comparison of Ixekizumab with Etanercept or Placebo in Moderate-to-Severe Psoriasis (UNCOVER-2 and UNCOVER-3): Results from Two Phase 3 Randomised Trials.” The Lancet 386 (9993): 541–51. https://doi.org/10.1016/s0140-6736(15)60125-8.
Langley, R. G., B. E. Elewski, M. Lebwohl, K. Reich, C. E. M. Griffiths, K. Papp, L. Puig, et al. 2014. “Secukinumab in Plaque Psoriasis — Results of Two Phase 3 Trials.” New England Journal of Medicine 371 (4): 326–38. https://doi.org/10.1056/nejmoa1314258.
Niederreiter, H. 1978. “Quasi-Monte Carlo Methods and Pseudo-Random Numbers.” Bulletin of the American Mathematical Society 84 (6): 957–1041. https://doi.org/10.1090/S0002-9904-1978-14532-7.
Phillippo, D. M. 2019. “Calibration of Treatment Effects in Network Meta-Analysis Using Individual Patient Data.” PhD thesis, University of Bristol.
Phillippo, D. M., S. Dias, A. E. Ades, M. Belger, A. Brnabic, A. Schacht, D. Saure, Z. Kadziola, and N. J. Welton. 2020. “Multilevel Network Meta-Regression for Population-Adjusted Treatment Comparisons.” Journal of the Royal Statistical Society: Series A (Statistics in Society) 183 (3): 1189–1210. https://doi.org/10.1111/rssa.12579.

1. The convergence rate of QMC is typically \(\mathcal{O}(1/n)\), whereas the expected convergence rate of standard MC is \(\mathcal{O}(1/n^\frac{1}{2})\) (Caflisch 1998; Niederreiter 1978).↩︎