Example: Plaque psoriasis ML-NMR

library(multinma)
library(dplyr)      # dplyr and tidyr for data manipulation
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.6496080 0.2917665 8.9328237 0.2167777 
#> 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.04    0.00 0.06    -0.08     0.00     0.05     0.09
#> beta[prevsys]                          -0.13    0.00 0.16    -0.45    -0.23    -0.13    -0.03
#> beta[bsa]                              -0.06    0.01 0.44    -0.98    -0.36    -0.04     0.24
#> beta[weight]                            0.04    0.00 0.03    -0.02     0.02     0.04     0.06
#> beta[psa]                              -0.07    0.00 0.16    -0.40    -0.19    -0.07     0.04
#> beta[durnpso:.trtclassTNFa blocker]    -0.03    0.00 0.07    -0.17    -0.08    -0.03     0.02
#> beta[durnpso:.trtclassIL blocker]      -0.01    0.00 0.07    -0.15    -0.06    -0.01     0.03
#> beta[prevsys:.trtclassTNFa blocker]     0.18    0.00 0.18    -0.17     0.06     0.18     0.31
#> beta[prevsys:.trtclassIL blocker]       0.06    0.00 0.18    -0.27    -0.06     0.06     0.18
#> beta[bsa:.trtclassTNFa blocker]         0.05    0.01 0.51    -0.93    -0.29     0.04     0.38
#> beta[bsa:.trtclassIL blocker]           0.29    0.01 0.50    -0.64    -0.05     0.26     0.61
#> beta[weight:.trtclassTNFa blocker]     -0.17    0.00 0.03    -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.20    -0.47    -0.19    -0.06     0.08
#> beta[psa:.trtclassIL blocker]           0.00    0.00 0.19    -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.96    0.00 0.08     2.79     2.90     2.96     3.01
#> d[IXE_Q4W]                              2.54    0.00 0.08     2.39     2.49     2.54     2.60
#> d[SEC_150]                              2.15    0.00 0.12     1.91     2.07     2.15     2.23
#> d[SEC_300]                              2.45    0.00 0.12     2.22     2.37     2.45     2.53
#> lp__                                -1576.40    0.09 3.54 -1584.12 -1578.65 -1576.15 -1573.83
#>                                        97.5% n_eff Rhat
#> beta[durnpso]                           0.16  5316 1.00
#> beta[prevsys]                           0.18  6873 1.00
#> beta[bsa]                               0.75  5554 1.00
#> beta[weight]                            0.09  5200 1.00
#> beta[psa]                               0.25  5914 1.00
#> beta[durnpso:.trtclassTNFa blocker]     0.11  5606 1.00
#> beta[durnpso:.trtclassIL blocker]       0.12  5777 1.00
#> beta[prevsys:.trtclassTNFa blocker]     0.55  6108 1.00
#> beta[prevsys:.trtclassIL blocker]       0.41  7685 1.00
#> beta[bsa:.trtclassTNFa blocker]         1.09  5958 1.00
#> beta[bsa:.trtclassIL blocker]           1.28  7103 1.00
#> beta[weight:.trtclassTNFa blocker]     -0.10  5812 1.00
#> beta[weight:.trtclassIL blocker]       -0.04  6067 1.00
#> beta[psa:.trtclassTNFa blocker]         0.35  6882 1.00
#> beta[psa:.trtclassIL blocker]           0.38  6710 1.00
#> d[ETN]                                  1.71  5403 1.00
#> d[IXE_Q2W]                              3.12  5068 1.00
#> d[IXE_Q4W]                              2.70  5665 1.00
#> d[SEC_150]                              2.38  5511 1.00
#> d[SEC_300]                              2.68  8290 1.00
#> lp__                                -1570.36  1478 1.01
#> 
#> Samples were drawn using NUTS(diag_e) at Tue Jan 18 09:45:05 2022.
#> 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 16 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.6496080 0.2917665 8.9328237 0.2167777 
#> 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.43    -0.23    -0.12    -0.01
#> beta[bsa]                              -0.10    0.01 0.44    -1.00    -0.39    -0.09     0.21
#> beta[weight]                            0.04    0.00 0.03    -0.01     0.02     0.04     0.06
#> beta[psa]                              -0.06    0.00 0.17    -0.41    -0.17    -0.06     0.05
#> beta[durnpso:.trtclassTNFa blocker]    -0.03    0.00 0.08    -0.18    -0.09    -0.03     0.02
#> beta[durnpso:.trtclassIL blocker]      -0.02    0.00 0.07    -0.15    -0.07    -0.02     0.03
#> beta[prevsys:.trtclassTNFa blocker]     0.18    0.00 0.18    -0.18     0.06     0.18     0.30
#> beta[prevsys:.trtclassIL blocker]       0.05    0.00 0.18    -0.30    -0.07     0.05     0.17
#> beta[bsa:.trtclassTNFa blocker]         0.08    0.01 0.52    -0.91    -0.28     0.09     0.44
#> beta[bsa:.trtclassIL blocker]           0.33    0.01 0.48    -0.58     0.01     0.32     0.65
#> 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.13    -0.10    -0.08
#> beta[psa:.trtclassTNFa blocker]        -0.07    0.00 0.20    -0.48    -0.20    -0.06     0.07
#> beta[psa:.trtclassIL blocker]          -0.01    0.00 0.18    -0.37    -0.13    -0.01     0.11
#> d[ETN]                                  1.56    0.00 0.16     1.26     1.48     1.56     1.64
#> d[IXE_Q2W]                              2.98    0.00 0.16     2.67     2.89     2.97     3.06
#> d[IXE_Q4W]                              2.56    0.00 0.16     2.24     2.47     2.56     2.65
#> d[SEC_150]                              2.14    0.01 0.25     1.62     2.01     2.13     2.27
#> d[SEC_300]                              2.44    0.00 0.24     1.94     2.30     2.44     2.57
#> lp__                                -1580.22    0.16 4.89 -1590.93 -1583.22 -1579.85 -1576.76
#> tau                                     0.20    0.01 0.13     0.02     0.10     0.17     0.26
#>                                        97.5% n_eff Rhat
#> beta[durnpso]                           0.17  4545 1.00
#> beta[prevsys]                           0.19  4618 1.00
#> beta[bsa]                               0.75  3808 1.00
#> beta[weight]                            0.10  5083 1.00
#> beta[psa]                               0.26  3919 1.00
#> beta[durnpso:.trtclassTNFa blocker]     0.12  4768 1.00
#> beta[durnpso:.trtclassIL blocker]       0.12  4879 1.00
#> beta[prevsys:.trtclassTNFa blocker]     0.54  4830 1.00
#> beta[prevsys:.trtclassIL blocker]       0.38  4587 1.00
#> beta[bsa:.trtclassTNFa blocker]         1.12  3955 1.00
#> beta[bsa:.trtclassIL blocker]           1.30  4441 1.00
#> beta[weight:.trtclassTNFa blocker]     -0.10  5216 1.00
#> beta[weight:.trtclassIL blocker]       -0.04  5021 1.00
#> beta[psa:.trtclassTNFa blocker]         0.33  3541 1.00
#> beta[psa:.trtclassIL blocker]           0.37  4126 1.00
#> d[ETN]                                  1.91  1684 1.00
#> d[IXE_Q2W]                              3.30  2103 1.00
#> d[IXE_Q4W]                              2.91  1155 1.00
#> d[SEC_150]                              2.62  1813 1.00
#> d[SEC_300]                              2.93  2353 1.00
#> lp__                                -1571.58   965 1.00
#> tau                                     0.54   570 1.01
#> 
#> Samples were drawn using NUTS(diag_e) at Tue Jan 18 09:54:52 2022.
#> 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.2
#>               DIC: 3153.8
(pso_dic_RE <- dic(pso_fit_RE))
#> Residual deviance: 3123.4 (on 3858 data points)
#>                pD: 28
#>               DIC: 3151.4

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.85     4253     3249    1
#> d[FIXTURE: IXE_Q2W] 3.03 0.10 2.84 2.96 3.03 3.10  3.23     4704     3054    1
#> d[FIXTURE: IXE_Q4W] 2.62 0.09 2.44 2.56 2.62 2.68  2.81     4982     3134    1
#> d[FIXTURE: SEC_150] 2.22 0.12 1.98 2.14 2.22 2.31  2.46     4919     3321    1
#> d[FIXTURE: SEC_300] 2.53 0.12 2.29 2.44 2.53 2.62  2.76     6773     3490    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.68     6669     3203    1
#> d[UNCOVER-1: IXE_Q2W] 2.92 0.08 2.76 2.87 2.92 2.98  3.09     5499     3396    1
#> d[UNCOVER-1: IXE_Q4W] 2.51 0.08 2.36 2.45 2.51 2.57  2.67     6437     3202    1
#> d[UNCOVER-1: SEC_150] 2.11 0.12 1.87 2.03 2.11 2.19  2.36     6546     2935    1
#> d[UNCOVER-1: SEC_300] 2.42 0.12 2.18 2.33 2.42 2.50  2.66     8678     2847    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.56  1.67     6310     2913    1
#> d[UNCOVER-2: IXE_Q2W] 2.92 0.08 2.76 2.87 2.92 2.98  3.09     5580     3393    1
#> d[UNCOVER-2: IXE_Q4W] 2.51 0.08 2.36 2.46 2.51 2.57  2.66     6391     3361    1
#> d[UNCOVER-2: SEC_150] 2.11 0.12 1.87 2.03 2.11 2.19  2.35     6158     3336    1
#> d[UNCOVER-2: SEC_300] 2.42 0.12 2.18 2.34 2.42 2.50  2.65     9278     3183    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.48 1.53 1.58  1.69     5599     3293    1
#> d[UNCOVER-3: IXE_Q2W] 2.94 0.08 2.77 2.88 2.94 3.00  3.11     5359     3660    1
#> d[UNCOVER-3: IXE_Q4W] 2.53 0.08 2.38 2.47 2.53 2.58  2.68     5988     3389    1
#> d[UNCOVER-3: SEC_150] 2.13 0.12 1.89 2.05 2.13 2.21  2.37     5775     3197    1
#> d[UNCOVER-3: SEC_300] 2.44 0.12 2.20 2.35 2.44 2.52  2.67     8638     3015    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     5154     3518    1
#> pred[FIXTURE: ETN]     0.46 0.03 0.41 0.44 0.46 0.47  0.51     6863     2897    1
#> pred[FIXTURE: IXE_Q2W] 0.89 0.02 0.85 0.88 0.89 0.90  0.92     7019     3061    1
#> pred[FIXTURE: IXE_Q4W] 0.80 0.03 0.74 0.78 0.80 0.81  0.84     7871     3134    1
#> pred[FIXTURE: SEC_150] 0.67 0.03 0.62 0.65 0.67 0.69  0.72     9539     3084    1
#> pred[FIXTURE: SEC_300] 0.77 0.02 0.72 0.76 0.77 0.79  0.81     9461     2592    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     6853     3363    1
#> pred[UNCOVER-1: ETN]     0.46 0.03 0.40 0.44 0.46 0.48  0.52    10415     3403    1
#> pred[UNCOVER-1: IXE_Q2W] 0.90 0.01 0.88 0.89 0.90 0.91  0.92     9008     2835    1
#> pred[UNCOVER-1: IXE_Q4W] 0.81 0.02 0.78 0.80 0.81 0.82  0.84     9584     3104    1
#> pred[UNCOVER-1: SEC_150] 0.69 0.04 0.60 0.66 0.69 0.72  0.77     6502     2747    1
#> pred[UNCOVER-1: SEC_300] 0.78 0.04 0.71 0.76 0.79 0.81  0.85     8172     3087    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     6395     3427    1
#> pred[UNCOVER-2: ETN]     0.42 0.02 0.38 0.40 0.42 0.43  0.46     8999     3406    1
#> pred[UNCOVER-2: IXE_Q2W] 0.88 0.01 0.86 0.87 0.88 0.89  0.91     8142     2704    1
#> pred[UNCOVER-2: IXE_Q4W] 0.78 0.02 0.75 0.77 0.78 0.79  0.81     8804     3251    1
#> pred[UNCOVER-2: SEC_150] 0.65 0.04 0.57 0.62 0.65 0.68  0.73     6696     3193    1
#> pred[UNCOVER-2: SEC_300] 0.75 0.04 0.68 0.73 0.76 0.78  0.82     8373     3321    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     6097     3318    1
#> pred[UNCOVER-3: ETN]     0.53 0.02 0.49 0.51 0.53 0.54  0.57     8825     3283    1
#> pred[UNCOVER-3: IXE_Q2W] 0.93 0.01 0.91 0.92 0.93 0.93  0.94     7193     3102    1
#> pred[UNCOVER-3: IXE_Q4W] 0.85 0.01 0.83 0.85 0.85 0.86  0.88     8614     3056    1
#> pred[UNCOVER-3: SEC_150] 0.75 0.04 0.67 0.72 0.75 0.77  0.81     6138     3269    1
#> pred[UNCOVER-3: SEC_300] 0.83 0.03 0.77 0.81 0.83 0.85  0.88     8289     3225    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.00    1   1   1   1     1       NA       NA   NA
#> rank[FIXTURE: IXE_Q4W] 2.24 0.43    2   2   2   2     3     4112     4016    1
#> rank[FIXTURE: SEC_150] 4.00 0.07    4   4   4   4     4     2784       NA    1
#> rank[FIXTURE: SEC_300] 2.77 0.43    2   3   3   3     3     4160     2662    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.00    1   1   1   1     1       NA       NA   NA
#> rank[UNCOVER-1: IXE_Q4W] 2.24 0.43    2   2   2   2     3     4112     4016    1
#> rank[UNCOVER-1: SEC_150] 4.00 0.07    4   4   4   4     4     2784       NA    1
#> rank[UNCOVER-1: SEC_300] 2.77 0.43    2   3   3   3     3     4160     2662    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.00    1   1   1   1     1       NA       NA   NA
#> rank[UNCOVER-2: IXE_Q4W] 2.24 0.43    2   2   2   2     3     4112     4016    1
#> rank[UNCOVER-2: SEC_150] 4.00 0.07    4   4   4   4     4     2784       NA    1
#> rank[UNCOVER-2: SEC_300] 2.77 0.43    2   3   3   3     3     4160     2662    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.00    1   1   1   1     1       NA       NA   NA
#> rank[UNCOVER-3: IXE_Q4W] 2.24 0.43    2   2   2   2     3     4112     4016    1
#> rank[UNCOVER-3: SEC_150] 4.00 0.07    4   4   4   4     4     2784       NA    1
#> rank[UNCOVER-3: SEC_300] 2.77 0.43    2   3   3   3     3     4160     2662    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.76      0.23         0         0         0
#> d[FIXTURE: SEC_150]         0      0.00      0.00         1         0         0
#> d[FIXTURE: SEC_300]         0      0.23      0.76         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.76      0.23         0         0         0
#> d[UNCOVER-1: SEC_150]         0      0.00      0.00         1         0         0
#> d[UNCOVER-1: SEC_300]         0      0.23      0.76         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.76      0.23         0         0         0
#> d[UNCOVER-2: SEC_150]         0      0.00      0.00         1         0         0
#> d[UNCOVER-2: SEC_300]         0      0.23      0.76         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.76      0.23         0         0         0
#> d[UNCOVER-3: SEC_150]         0      0.00      0.00         1         0         0
#> d[UNCOVER-3: SEC_300]         0      0.23      0.76         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.76         1         1         1         1
#> d[FIXTURE: SEC_150]         0      0.00         0         1         1         1
#> d[FIXTURE: SEC_300]         0      0.23         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.76         1         1         1         1
#> d[UNCOVER-1: SEC_150]         0      0.00         0         1         1         1
#> d[UNCOVER-1: SEC_300]         0      0.23         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.76         1         1         1         1
#> d[UNCOVER-2: SEC_150]         0      0.00         0         1         1         1
#> d[UNCOVER-2: SEC_300]         0      0.23         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.76         1         1         1         1
#> d[UNCOVER-3: SEC_150]         0      0.00         0         1         1         1
#> d[UNCOVER-3: SEC_300]         0      0.23         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.22 0.82 1.11 1.26 1.41  1.69     5838     2848    1
#> d[New 1: IXE_Q2W] 2.89 0.22 2.45 2.74 2.89 3.04  3.34     7185     3196    1
#> d[New 1: IXE_Q4W] 2.48 0.22 2.05 2.34 2.47 2.62  2.92     7206     3097    1
#> d[New 1: SEC_150] 2.08 0.22 1.65 1.93 2.08 2.23  2.54     6054     2842    1
#> d[New 1: SEC_300] 2.39 0.23 1.96 2.23 2.39 2.53  2.85     7254     2871    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.02 0.03 0.04 0.06 0.07  0.12     5686     3671    1
#> pred[New 1: ETN]     0.37 0.06 0.26 0.33 0.37 0.41  0.49     6683     3562    1
#> pred[New 1: IXE_Q2W] 0.90 0.03 0.84 0.88 0.90 0.92  0.94     5943     3386    1
#> pred[New 1: IXE_Q4W] 0.81 0.04 0.73 0.78 0.81 0.83  0.87     6026     3578    1
#> pred[New 1: SEC_150] 0.68 0.06 0.57 0.64 0.68 0.72  0.79     5629     3392    1
#> pred[New 1: SEC_300] 0.78 0.05 0.68 0.75 0.78 0.81  0.86     6655     3595    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).↩︎