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) and Phillippo et al. (2022; see also Phillippo 2019).

In the first analysis (Phillippo et al. 2020), we consider a network of four studies with a binary outcome (success/failure to achieve a 75% reduction on the psoriasis area and severity index, PASI 75).

In the second analysis (Phillippo et al. 2022), we extend this network with a further five studies and demonstrate how the key assumptions of population adjustment can be assessed in this larger network. We also demonstrate how to produce estimates for three external target populations, and fit a multinomial model to incorporate ordered categorical outcomes (PASI 75, PASI 90, and PASI 100).

Initial analysis

We start by recreating the analysis presented by Phillippo et al. (2020). 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     Treatment arms                  
#>  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   Treatment arms                  
#>  FIXTURE 4: PBO | ETN | 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.¹ 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.07     0.00     0.05     0.09
#> beta[prevsys]                          -0.13    0.00 0.16    -0.44    -0.24    -0.13    -0.02
#> beta[bsa]                              -0.06    0.01 0.44    -0.94    -0.36    -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.42    -0.19    -0.08     0.04
#> 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.14    -0.06    -0.01     0.03
#> beta[prevsys:.trtclassTNFa blocker]     0.18    0.00 0.19    -0.18     0.05     0.18     0.31
#> beta[prevsys:.trtclassIL blocker]       0.06    0.00 0.17    -0.29    -0.06     0.06     0.18
#> beta[bsa:.trtclassTNFa blocker]         0.05    0.01 0.51    -0.93    -0.30     0.03     0.39
#> beta[bsa:.trtclassIL blocker]           0.29    0.01 0.48    -0.63    -0.04     0.28     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.05    0.00 0.20    -0.45    -0.19    -0.06     0.09
#> beta[psa:.trtclassIL blocker]           0.01    0.00 0.19    -0.36    -0.12     0.01     0.13
#> d[ETN]                                  1.55    0.00 0.08     1.39     1.49     1.55     1.60
#> d[IXE_Q2W]                              2.95    0.00 0.09     2.79     2.90     2.95     3.01
#> d[IXE_Q4W]                              2.54    0.00 0.08     2.39     2.49     2.54     2.60
#> d[SEC_150]                              2.14    0.00 0.11     1.91     2.07     2.14     2.22
#> d[SEC_300]                              2.45    0.00 0.12     2.22     2.37     2.45     2.53
#> lp__                                -1576.33    0.08 3.42 -1583.99 -1578.45 -1575.95 -1573.86
#>                                        97.5% n_eff Rhat
#> beta[durnpso]                           0.16  5352    1
#> beta[prevsys]                           0.19  6199    1
#> beta[bsa]                               0.77  4676    1
#> beta[weight]                            0.10  6123    1
#> beta[psa]                               0.23  5817    1
#> beta[durnpso:.trtclassTNFa blocker]     0.12  5105    1
#> beta[durnpso:.trtclassIL blocker]       0.12  6136    1
#> beta[prevsys:.trtclassTNFa blocker]     0.54  6347    1
#> beta[prevsys:.trtclassIL blocker]       0.39  7872    1
#> beta[bsa:.trtclassTNFa blocker]         1.05  5101    1
#> beta[bsa:.trtclassIL blocker]           1.23  5477    1
#> beta[weight:.trtclassTNFa blocker]     -0.10  6300    1
#> beta[weight:.trtclassIL blocker]       -0.04  7778    1
#> beta[psa:.trtclassTNFa blocker]         0.35  6504    1
#> beta[psa:.trtclassIL blocker]           0.38  6998    1
#> d[ETN]                                  1.70  4713    1
#> d[IXE_Q2W]                              3.13  5290    1
#> d[IXE_Q4W]                              2.70  5589    1
#> d[SEC_150]                              2.37  5583    1
#> d[SEC_300]                              2.68  6434    1
#> lp__                                -1570.58  1722    1
#> 
#> Samples were drawn using NUTS(diag_e) at Mon Aug 29 16:48:26 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 17 divergent transitions after warmup. See
#> https://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.00     0.05     0.09
#> beta[prevsys]                          -0.12    0.00 0.16    -0.44    -0.23    -0.12    -0.01
#> beta[bsa]                              -0.10    0.01 0.46    -1.01    -0.42    -0.08     0.22
#> 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.40    -0.18    -0.06     0.05
#> beta[durnpso:.trtclassTNFa blocker]    -0.03    0.00 0.07    -0.19    -0.08    -0.03     0.02
#> beta[durnpso:.trtclassIL blocker]      -0.02    0.00 0.07    -0.14    -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.05    0.00 0.18    -0.31    -0.08     0.05     0.17
#> beta[bsa:.trtclassTNFa blocker]         0.09    0.01 0.53    -0.94    -0.28     0.07     0.45
#> beta[bsa:.trtclassIL blocker]           0.34    0.01 0.50    -0.62    -0.01     0.32     0.69
#> beta[weight:.trtclassTNFa blocker]     -0.17    0.00 0.04    -0.24    -0.19    -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.21    -0.46    -0.21    -0.08     0.07
#> beta[psa:.trtclassIL blocker]          -0.01    0.00 0.19    -0.37    -0.14    -0.01     0.11
#> d[ETN]                                  1.56    0.00 0.15     1.28     1.47     1.56     1.65
#> d[IXE_Q2W]                              2.98    0.00 0.15     2.69     2.88     2.97     3.06
#> d[IXE_Q4W]                              2.57    0.00 0.16     2.28     2.47     2.56     2.65
#> d[SEC_150]                              2.14    0.01 0.23     1.67     2.01     2.14     2.27
#> d[SEC_300]                              2.44    0.01 0.23     1.98     2.31     2.44     2.58
#> lp__                                -1580.50    0.17 4.86 -1590.80 -1583.66 -1580.20 -1577.06
#> tau                                     0.19    0.01 0.12     0.02     0.10     0.17     0.25
#>                                        97.5% n_eff Rhat
#> beta[durnpso]                           0.17  3933 1.00
#> beta[prevsys]                           0.20  4787 1.00
#> beta[bsa]                               0.74  3828 1.00
#> beta[weight]                            0.10  4359 1.00
#> beta[psa]                               0.26  3843 1.00
#> beta[durnpso:.trtclassTNFa blocker]     0.11  4361 1.00
#> beta[durnpso:.trtclassIL blocker]       0.12  4262 1.00
#> beta[prevsys:.trtclassTNFa blocker]     0.54  5022 1.00
#> beta[prevsys:.trtclassIL blocker]       0.39  6186 1.00
#> beta[bsa:.trtclassTNFa blocker]         1.15  3932 1.00
#> beta[bsa:.trtclassIL blocker]           1.34  4581 1.00
#> beta[weight:.trtclassTNFa blocker]     -0.10  4496 1.00
#> beta[weight:.trtclassIL blocker]       -0.04  4718 1.00
#> beta[psa:.trtclassTNFa blocker]         0.34  3944 1.00
#> beta[psa:.trtclassIL blocker]           0.35  4874 1.00
#> d[ETN]                                  1.88  1711 1.00
#> d[IXE_Q2W]                              3.31  1544 1.00
#> d[IXE_Q4W]                              2.93   995 1.01
#> d[SEC_150]                              2.59  1612 1.00
#> d[SEC_300]                              2.91  1481 1.00
#> lp__                                -1571.85   860 1.00
#> tau                                     0.49   521 1.00
#> 
#> Samples were drawn using NUTS(diag_e) at Mon Aug 29 16:56:04 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.5 (on 3858 data points)
#>                pD: 24.1
#>               DIC: 3153.6
(pso_dic_RE <- dic(pso_fit_RE))
#> Residual deviance: 3123.8 (on 3858 data points)
#>                pD: 28.5
#>               DIC: 3152.3

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

Producing relative effects and event probabilities

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.72  1.84     4000     3262    1
#> d[FIXTURE: IXE_Q2W] 3.03 0.10 2.84 2.96 3.03 3.09  3.23     4853     2888    1
#> d[FIXTURE: IXE_Q4W] 2.62 0.09 2.43 2.55 2.61 2.68  2.80     5046     3299    1
#> d[FIXTURE: SEC_150] 2.22 0.12 1.98 2.14 2.22 2.30  2.45     4940     3108    1
#> d[FIXTURE: SEC_300] 2.53 0.12 2.29 2.45 2.52 2.61  2.76     5675     3209    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.50 0.08 1.34 1.45 1.50 1.56  1.67     5925     3213 1.00
#> d[UNCOVER-1: IXE_Q2W] 2.92 0.08 2.76 2.86 2.92 2.98  3.09     5675     3282 1.01
#> d[UNCOVER-1: IXE_Q4W] 2.51 0.08 2.35 2.45 2.51 2.56  2.67     5959     3148 1.00
#> d[UNCOVER-1: SEC_150] 2.11 0.12 1.88 2.03 2.11 2.19  2.34     5856     3244 1.00
#> d[UNCOVER-1: SEC_300] 2.42 0.12 2.18 2.33 2.42 2.50  2.66     6771     3018 1.00
#> 
#> -------------------------------------------------------------- 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.50 0.08 1.35 1.45 1.50 1.56  1.66     5634     3188    1
#> d[UNCOVER-2: IXE_Q2W] 2.92 0.08 2.76 2.86 2.92 2.98  3.09     5703     3337    1
#> d[UNCOVER-2: IXE_Q4W] 2.51 0.08 2.36 2.45 2.51 2.56  2.66     6065     3029    1
#> d[UNCOVER-2: SEC_150] 2.11 0.12 1.87 2.03 2.11 2.19  2.34     5946     3111    1
#> d[UNCOVER-2: SEC_300] 2.42 0.12 2.19 2.34 2.42 2.49  2.66     6838     3048    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.58  1.68     5088     3070    1
#> d[UNCOVER-3: IXE_Q2W] 2.94 0.09 2.77 2.88 2.94 3.00  3.11     5536     3086    1
#> d[UNCOVER-3: IXE_Q4W] 2.53 0.08 2.37 2.47 2.52 2.58  2.68     5850     3315    1
#> d[UNCOVER-3: SEC_150] 2.13 0.11 1.89 2.05 2.13 2.21  2.35     5822     3164    1
#> d[UNCOVER-3: SEC_300] 2.43 0.12 2.21 2.36 2.44 2.51  2.67     6709     3149    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     4472     3341    1
#> pred[FIXTURE: ETN]     0.46 0.03 0.41 0.44 0.46 0.47  0.51     6925     3171    1
#> pred[FIXTURE: IXE_Q2W] 0.89 0.02 0.85 0.88 0.89 0.90  0.92     6743     3099    1
#> pred[FIXTURE: IXE_Q4W] 0.80 0.03 0.74 0.78 0.80 0.81  0.84     7415     3015    1
#> pred[FIXTURE: SEC_150] 0.67 0.03 0.62 0.65 0.67 0.69  0.72     9012     2567    1
#> pred[FIXTURE: SEC_300] 0.77 0.02 0.73 0.76 0.77 0.79  0.81     8914     3352    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.07     6078     3247    1
#> pred[UNCOVER-1: ETN]     0.46 0.03 0.41 0.44 0.46 0.48  0.52     9265     3219    1
#> pred[UNCOVER-1: IXE_Q2W] 0.90 0.01 0.88 0.89 0.90 0.91  0.92     7819     3198    1
#> pred[UNCOVER-1: IXE_Q4W] 0.81 0.02 0.78 0.80 0.81 0.82  0.84    10593     3414    1
#> pred[UNCOVER-1: SEC_150] 0.69 0.04 0.60 0.66 0.69 0.72  0.77     7116     3212    1
#> pred[UNCOVER-1: SEC_300] 0.78 0.03 0.71 0.76 0.79 0.81  0.85     7934     2987    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     6252     3015    1
#> pred[UNCOVER-2: ETN]     0.42 0.02 0.38 0.40 0.42 0.43  0.46     8076     3303    1
#> pred[UNCOVER-2: IXE_Q2W] 0.88 0.01 0.86 0.87 0.88 0.89  0.91     7216     3314    1
#> pred[UNCOVER-2: IXE_Q4W] 0.78 0.02 0.75 0.77 0.78 0.79  0.81     9341     2965    1
#> pred[UNCOVER-2: SEC_150] 0.65 0.04 0.56 0.62 0.65 0.68  0.73     6898     2925    1
#> pred[UNCOVER-2: SEC_300] 0.75 0.04 0.68 0.73 0.75 0.78  0.82     6887     2928    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     5962     3081    1
#> pred[UNCOVER-3: ETN]     0.53 0.02 0.49 0.51 0.53 0.54  0.57     9197     3154    1
#> pred[UNCOVER-3: IXE_Q2W] 0.93 0.01 0.91 0.92 0.93 0.93  0.94     7615     3408    1
#> pred[UNCOVER-3: IXE_Q4W] 0.85 0.01 0.83 0.85 0.85 0.86  0.88     8008     3068    1
#> pred[UNCOVER-3: SEC_150] 0.75 0.04 0.67 0.72 0.75 0.77  0.81     6728     2660    1
#> pred[UNCOVER-3: SEC_300] 0.83 0.03 0.77 0.81 0.83 0.85  0.88     7529     2983    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.22 0.41    2   2   2   2     3     4317     4018    1
#> rank[FIXTURE: SEC_150] 4.00 0.05    4   4   4   4     4     4033       NA    1
#> rank[FIXTURE: SEC_300] 2.78 0.42    2   3   3   3     3     4265     4029    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.22 0.41    2   2   2   2     3     4317     4018    1
#> rank[UNCOVER-1: SEC_150] 4.00 0.05    4   4   4   4     4     4033       NA    1
#> rank[UNCOVER-1: SEC_300] 2.78 0.42    2   3   3   3     3     4265     4029    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.22 0.41    2   2   2   2     3     4317     4018    1
#> rank[UNCOVER-2: SEC_150] 4.00 0.05    4   4   4   4     4     4033       NA    1
#> rank[UNCOVER-2: SEC_300] 2.78 0.42    2   3   3   3     3     4265     4029    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.22 0.41    2   2   2   2     3     4317     4018    1
#> rank[UNCOVER-3: SEC_150] 4.00 0.05    4   4   4   4     4     4033       NA    1
#> rank[UNCOVER-3: SEC_300] 2.78 0.42    2   3   3   3     3     4265     4029    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.78         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.78         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.78         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.78         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.81 1.09 1.25 1.40  1.72     5985     3283    1
#> d[New 1: IXE_Q2W] 2.89 0.22 2.48 2.74 2.89 3.03  3.33     6586     3380    1
#> d[New 1: IXE_Q4W] 2.48 0.22 2.06 2.33 2.48 2.62  2.91     6686     3309    1
#> d[New 1: SEC_150] 2.08 0.22 1.66 1.93 2.08 2.23  2.52     6643     3190    1
#> d[New 1: SEC_300] 2.39 0.22 1.96 2.23 2.38 2.53  2.84     6800     3009    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.08  0.12     5452     3179    1
#> pred[New 1: ETN]     0.37 0.06 0.26 0.33 0.36 0.41  0.49     6094     3553    1
#> pred[New 1: IXE_Q2W] 0.90 0.03 0.84 0.88 0.90 0.92  0.94     5384     3695    1
#> pred[New 1: IXE_Q4W] 0.81 0.04 0.72 0.78 0.81 0.83  0.87     5340     3204    1
#> pred[New 1: SEC_150] 0.68 0.06 0.57 0.64 0.68 0.72  0.78     5221     3012    1
#> pred[New 1: SEC_300] 0.78 0.05 0.68 0.75 0.78 0.81  0.86     5763     3147    1
plot(pso_pred_FE_new, ref_line = c(0, 1))

Extended analysis

We now extend the network to include a further five studies (four AgD and one IPD), recreating the analysis of Phillippo et al. (2022). This larger network allows us to assess the key assumptions underlying population adjustment.

Setup

Preparing the data

We begin, as before, with some data transformations for each of the covariates and set up a treatment class variable trtclass.

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

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

There are a very small number of individuals with missing values in the IPD, which we simply exclude from the analysis.

pso_ipd %>% 
  group_by(studyc) %>% 
  summarise(n_total = n(),
            n_missing = sum(!is_complete), 
            pct_missing = mean(!is_complete) * 100)
#> # A tibble: 4 × 4
#>   studyc    n_total n_missing pct_missing
#>   <chr>       <int>     <int>       <dbl>
#> 1 IXORA-S       260         0       0    
#> 2 UNCOVER-1    1296         0       0    
#> 3 UNCOVER-2    1221         2       0.164
#> 4 UNCOVER-3    1341         2       0.149

pso_ipd <- filter(pso_ipd, is_complete)

Creating the network

Next we set up the network. We set the IPD with set_ipd() and AgD (arm-based) with set_agd_arm(), and combine these together using combine_network(). We specify an ordered categorical (multinomial) outcome using the multi() helper function. The outcome data are in “inclusive” format, i.e. the lowest category is the sample size (or 1 for IPD), the second category counts those achieving PASI 75 or greater (\(\ge 75\%\) reduction in symptoms), the third counts those achieving PASI 90 or greater (\(\ge 90\%\) reduction), and the final category counts those achieving PASI 100 (\(100\%\) reduction).² We specify the treatment classes with trt_class = trtclass.

pso_net <- combine_network(
  set_ipd(pso_ipd,
    study = studyc,
    trt = trtc,
    r = multi(r0 = 1, 
              PASI75 = pasi75,
              PASI90 = pasi90,
              PASI100 = pasi100,
              type = "ordered", inclusive = TRUE),
    trt_class = trtclass),
  set_agd_arm(pso_agd,
    study = studyc,
    trt = trtc,
    r = multi(r0 = pasi75_n, 
              PASI75 = pasi75_r,
              PASI90 = pasi90_r,
              PASI100 = pasi100_r,
              type = "ordered", inclusive = TRUE),
    trt_class = trtclass)
)

pso_net
#> A network with 4 IPD studies, and 5 AgD studies (arm-based).
#> 
#> ------------------------------------------------------------------- IPD studies ---- 
#>  Study     Treatment arms                  
#>  IXORA-S   2: IXE_Q2W | UST                
#>  UNCOVER-1 3: PBO | IXE_Q2W | IXE_Q4W      
#>  UNCOVER-2 4: PBO | IXE_Q2W | IXE_Q4W | ETN
#>  UNCOVER-3 4: PBO | IXE_Q2W | IXE_Q4W | ETN
#> 
#>  Outcome type: ordered (4 categories)
#> ------------------------------------------------------- AgD studies (arm-based) ---- 
#>  Study    Treatment arms                  
#>  CLEAR    2: SEC_300 | UST                
#>  ERASURE  3: PBO | SEC_150 | SEC_300      
#>  FEATURE  3: PBO | SEC_150 | SEC_300      
#>  FIXTURE  4: PBO | ETN | SEC_150 | SEC_300
#>  JUNCTURE 3: PBO | SEC_150 | SEC_300      
#> 
#>  Outcome type: ordered (4 categories)
#> ------------------------------------------------------------------------------------
#> Total number of treatments: 7, in 4 classes
#> Total number of studies: 9
#> Reference treatment is: PBO
#> Network is connected

We create a network plot using the plot() function applied to the pso_net network object, choosing to scale the edges and nodes by the number of studies/sample size (weight_edges and weight_nodes = TRUE), colour the treatment nodes by class (show_trt_class = TRUE), and nudge the treatment names away from the nodes (nudge = 0.1). We further customise the plot using ggplot syntax to alter the colour scheme.

class_pal <- c("#D95F02", "#7570B3", "#E7298A", "#E6AB02")

plot(pso_net, weight_nodes = TRUE, weight_edges = TRUE, show_trt_class = TRUE, nudge = 0.1) +
  ggraph::scale_edge_colour_manual("Data", 
                                   values = c(AgD = "#113259", IPD = "#55A480"),
                                   guide = guide_legend(override.aes = list(edge_width = 2))) +
  scale_fill_manual("Treatment class", 
                    values = class_pal,
                    aesthetics = c("fill", "colour"),
                    guide = guide_legend(override.aes = list(size = 2)))
#> Scale for 'edge_colour' is already present. Adding another scale for 'edge_colour', which will
#> replace the existing scale.
#> Scale for 'fill' is already present. Adding another scale for 'fill', which will replace the
#> existing scale.
#> Warning: Duplicated override.aes is ignored.

Numerical integration for ML-NMR models

We add integration points to the AgD studies in the network using the add_integration() function, specifying the chosen marginal distribution for each covariate using the distr() function. As before, we specify Gamma distributions for weight and duration of psoriasis, a logit-Normal distribution for body surface area, and Bernoulli distributions for previous systemic treatment and psoriatic arthritis as binary covariates. Since we do not know the correlations between covariates in the AgD studies, we once again 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.

ML-NMR model

Using the nma() function, we fit a (fixed effect) ML-NMR model which includes main effects (prognostic terms) and covariate-treatment interactions (effect-modifying terms) for each of the five covariates. Ideally, we would fit independent interaction terms for each treatment; however, this requires either IPD or several AgD studies at a range of covariate values on each treatment. The data here are insufficient to fit independent interaction terms for each treatment, so we make the shared effect modifier assumption within each class of treatments (Phillippo et al. 2016) and specify common interaction terms within each treatment class (class_interactions = "common"). As before, we specify \(\mathrm{N}(0, 10^2)\) prior distributions on the study-specific intercepts, treatment effects, and regression parameters. However, since we now have an ordered multinomial likelihood we also need to specify priors for the differences between the latent cutoffs for each outcome category; we choose an improper flat prior \(\mathrm{U}(-\infty,\infty)\) which will automatically be truncated to meet the ordering constraints (prior_aux = flat()).

pso_fit_FE <- nma(pso_net, 
                  trt_effects = "fixed",
                  link = "probit", 
                  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_aux = flat(),
                  QR = TRUE,
                  init_r = 0.5)

#> Note: Setting "PBO" as the network reference treatment.

pso_fit_FE

#> A fixed effects ML-NMR with a ordered 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.8159830 0.6576489 0.2987820 8.9097263 0.2104826 
#> Inference for Stan model: ordered_multinomial.
#> 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.03    0.00 0.06    -0.08    -0.01     0.03     0.08
#> beta[prevsys]                              -0.17    0.00 0.16    -0.49    -0.28    -0.17    -0.06
#> beta[bsa]                                  -0.10    0.01 0.45    -1.01    -0.40    -0.08     0.21
#> beta[weight]                                0.04    0.00 0.03    -0.01     0.02     0.04     0.06
#> beta[psa]                                  -0.08    0.00 0.17    -0.42    -0.20    -0.08     0.03
#> beta[durnpso:.trtclassTNFa blocker]        -0.02    0.00 0.07    -0.16    -0.07    -0.02     0.03
#> beta[durnpso:.trtclassIL-12/23 blocker]    -0.06    0.00 0.10    -0.26    -0.13    -0.07     0.00
#> beta[durnpso:.trtclassIL-17 blocker]       -0.02    0.00 0.06    -0.14    -0.07    -0.02     0.02
#> beta[prevsys:.trtclassTNFa blocker]         0.19    0.00 0.19    -0.18     0.07     0.19     0.32
#> beta[prevsys:.trtclassIL-12/23 blocker]     0.45    0.01 0.34    -0.24     0.24     0.47     0.68
#> beta[prevsys:.trtclassIL-17 blocker]        0.16    0.00 0.17    -0.17     0.05     0.16     0.28
#> beta[bsa:.trtclassTNFa blocker]             0.24    0.01 0.51    -0.73    -0.11     0.22     0.57
#> beta[bsa:.trtclassIL-12/23 blocker]         0.61    0.01 0.68    -0.71     0.14     0.61     1.07
#> beta[bsa:.trtclassIL-17 blocker]            0.27    0.01 0.47    -0.63    -0.06     0.26     0.58
#> beta[weight:.trtclassTNFa blocker]         -0.16    0.00 0.03    -0.23    -0.18    -0.16    -0.14
#> beta[weight:.trtclassIL-12/23 blocker]     -0.09    0.00 0.05    -0.18    -0.12    -0.09    -0.06
#> beta[weight:.trtclassIL-17 blocker]        -0.13    0.00 0.03    -0.19    -0.15    -0.13    -0.11
#> beta[psa:.trtclassTNFa blocker]            -0.05    0.00 0.20    -0.43    -0.18    -0.04     0.08
#> beta[psa:.trtclassIL-12/23 blocker]         0.12    0.01 0.33    -0.51    -0.11     0.12     0.36
#> beta[psa:.trtclassIL-17 blocker]            0.10    0.00 0.18    -0.24    -0.02     0.10     0.21
#> d[ETN]                                      1.58    0.00 0.07     1.44     1.53     1.58     1.63
#> d[IXE_Q2W]                                  2.91    0.00 0.07     2.76     2.86     2.91     2.96
#> d[IXE_Q4W]                                  2.69    0.00 0.08     2.54     2.64     2.69     2.74
#> d[SEC_150]                                  2.19    0.00 0.08     2.03     2.13     2.19     2.24
#> d[SEC_300]                                  2.60    0.00 0.08     2.45     2.54     2.60     2.65
#> d[UST]                                      2.13    0.00 0.11     1.91     2.06     2.13     2.20
#> lp__                                    -7640.27    0.11 4.31 -7649.63 -7642.96 -7639.97 -7637.16
#> cc[PASI75]                                  0.00     NaN 0.00     0.00     0.00     0.00     0.00
#> cc[PASI90]                                  0.69    0.00 0.02     0.65     0.68     0.69     0.70
#> cc[PASI100]                                 1.53    0.00 0.02     1.49     1.52     1.53     1.55
#>                                            97.5% n_eff Rhat
#> beta[durnpso]                               0.15  2512    1
#> beta[prevsys]                               0.14  2701    1
#> beta[bsa]                                   0.77  2254    1
#> beta[weight]                                0.10  2331    1
#> beta[psa]                                   0.24  3019    1
#> beta[durnpso:.trtclassTNFa blocker]         0.12  2678    1
#> beta[durnpso:.trtclassIL-12/23 blocker]     0.14  3490    1
#> beta[durnpso:.trtclassIL-17 blocker]        0.10  2899    1
#> beta[prevsys:.trtclassTNFa blocker]         0.55  2824    1
#> beta[prevsys:.trtclassIL-12/23 blocker]     1.11  4293    1
#> beta[prevsys:.trtclassIL-17 blocker]        0.48  3128    1
#> beta[bsa:.trtclassTNFa blocker]             1.27  2433    1
#> beta[bsa:.trtclassIL-12/23 blocker]         1.93  2921    1
#> beta[bsa:.trtclassIL-17 blocker]            1.25  2728    1
#> beta[weight:.trtclassTNFa blocker]         -0.10  2590    1
#> beta[weight:.trtclassIL-12/23 blocker]      0.01  3473    1
#> beta[weight:.trtclassIL-17 blocker]        -0.07  2731    1
#> beta[psa:.trtclassTNFa blocker]             0.33  2990    1
#> beta[psa:.trtclassIL-12/23 blocker]         0.77  3967    1
#> beta[psa:.trtclassIL-17 blocker]            0.46  3516    1
#> d[ETN]                                      1.72  1969    1
#> d[IXE_Q2W]                                  3.06  2104    1
#> d[IXE_Q4W]                                  2.84  2327    1
#> d[SEC_150]                                  2.36  2267    1
#> d[SEC_300]                                  2.76  2396    1
#> d[UST]                                      2.34  3215    1
#> lp__                                    -7632.94  1664    1
#> cc[PASI75]                                  0.00   NaN  NaN
#> cc[PASI90]                                  0.72  3323    1
#> cc[PASI100]                                 1.58  2862    1
#> 
#> Samples were drawn using NUTS(diag_e) at Sun Aug 28 12:56:36 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).

Assessing assumptions

In the first analysis, the small network made assessing assumptions difficult. With this larger network (although still only nine studies) we have greater opportunity to assess the key assumptions.

The key assumption made by ML-NMR (and indeed all population adjustment methods in connected networks) is the conditional constancy of relative effects assumption (Phillippo et al. 2016). This means that there are no unobserved effect modifiers, so that the relative effects are constant given the included effect-modifying covariates. This assumption implies that there is no residual heterogeneity or inconsistency, which can be assessed using standard network meta-analysis techniques. We assess residual heterogeneity using a random effects model, and residual inconsistency using an unrelated mean effects (UME) model.

Assessing residual heterogeneity with a random effects ML-NMR

First, we fit a random effects model to assess residual heterogeneity. The call to the nma() function is identical to the fixed effect model above, except that now we specify trt_effects = "random" and need to provide a prior for the between-study heterogeneity (we choose a \(\textrm{half-N}(0, 2.5^2)\) prior with prior_het = half_normal(scale = 2.5).

pso_fit_RE <- nma(pso_net, 
                  trt_effects = "random",
                  link = "probit", 
                  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_aux = flat(),
                  prior_het = half_normal(scale = 2.5),
                  QR = TRUE,
                  init_r = 0.5)

#> Note: Setting "PBO" as the network reference treatment.
#> Warning: There were 1 divergent transitions after warmup. See
#> https://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

pso_fit_RE

#> A random effects ML-NMR with a ordered 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.8159830 0.6576489 0.2987820 8.9097263 0.2104826 
#> Inference for Stan model: ordered_multinomial.
#> 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.01     0.04     0.08
#> beta[prevsys]                              -0.16    0.00 0.16    -0.46    -0.26    -0.15    -0.06
#> beta[bsa]                                  -0.14    0.01 0.46    -1.09    -0.45    -0.13     0.17
#> beta[weight]                                0.05    0.00 0.03    -0.01     0.03     0.05     0.07
#> beta[psa]                                  -0.07    0.00 0.17    -0.40    -0.18    -0.06     0.04
#> beta[durnpso:.trtclassTNFa blocker]        -0.02    0.00 0.07    -0.17    -0.07    -0.02     0.03
#> beta[durnpso:.trtclassIL-12/23 blocker]    -0.07    0.00 0.10    -0.27    -0.14    -0.07     0.01
#> beta[durnpso:.trtclassIL-17 blocker]       -0.03    0.00 0.06    -0.15    -0.07    -0.02     0.02
#> beta[prevsys:.trtclassTNFa blocker]         0.19    0.00 0.18    -0.17     0.07     0.19     0.31
#> beta[prevsys:.trtclassIL-12/23 blocker]     0.43    0.00 0.35    -0.27     0.20     0.44     0.67
#> beta[prevsys:.trtclassIL-17 blocker]        0.15    0.00 0.16    -0.18     0.04     0.15     0.26
#> beta[bsa:.trtclassTNFa blocker]             0.27    0.01 0.53    -0.73    -0.08     0.25     0.61
#> beta[bsa:.trtclassIL-12/23 blocker]         0.66    0.01 0.66    -0.65     0.21     0.67     1.09
#> beta[bsa:.trtclassIL-17 blocker]            0.32    0.01 0.48    -0.58    -0.01     0.31     0.63
#> beta[weight:.trtclassTNFa blocker]         -0.16    0.00 0.03    -0.23    -0.19    -0.16    -0.14
#> beta[weight:.trtclassIL-12/23 blocker]     -0.09    0.00 0.05    -0.18    -0.12    -0.09    -0.06
#> beta[weight:.trtclassIL-17 blocker]        -0.13    0.00 0.03    -0.19    -0.15    -0.13    -0.11
#> beta[psa:.trtclassTNFa blocker]            -0.06    0.00 0.20    -0.46    -0.20    -0.06     0.07
#> beta[psa:.trtclassIL-12/23 blocker]         0.11    0.00 0.33    -0.53    -0.10     0.11     0.34
#> beta[psa:.trtclassIL-17 blocker]            0.08    0.00 0.18    -0.27    -0.04     0.08     0.20
#> d[ETN]                                      1.59    0.00 0.11     1.38     1.52     1.59     1.66
#> d[IXE_Q2W]                                  2.93    0.00 0.11     2.73     2.86     2.93     3.00
#> d[IXE_Q4W]                                  2.71    0.00 0.12     2.49     2.64     2.71     2.78
#> d[SEC_150]                                  2.21    0.00 0.12     1.99     2.13     2.21     2.28
#> d[SEC_300]                                  2.64    0.00 0.12     2.42     2.56     2.63     2.71
#> d[UST]                                      2.17    0.00 0.17     1.85     2.06     2.16     2.27
#> lp__                                    -7646.77    0.21 6.32 -7660.09 -7650.88 -7646.44 -7642.40
#> tau                                         0.13    0.00 0.07     0.02     0.09     0.13     0.17
#> cc[PASI75]                                  0.00     NaN 0.00     0.00     0.00     0.00     0.00
#> cc[PASI90]                                  0.69    0.00 0.02     0.66     0.68     0.69     0.70
#> cc[PASI100]                                 1.53    0.00 0.02     1.49     1.52     1.53     1.55
#>                                            97.5% n_eff Rhat
#> beta[durnpso]                               0.16  3817    1
#> beta[prevsys]                               0.16  4212    1
#> beta[bsa]                                   0.73  3325    1
#> beta[weight]                                0.10  3585    1
#> beta[psa]                                   0.27  4167    1
#> beta[durnpso:.trtclassTNFa blocker]         0.12  3797    1
#> beta[durnpso:.trtclassIL-12/23 blocker]     0.13  5280    1
#> beta[durnpso:.trtclassIL-17 blocker]        0.10  4181    1
#> beta[prevsys:.trtclassTNFa blocker]         0.54  4064    1
#> beta[prevsys:.trtclassIL-12/23 blocker]     1.06  5604    1
#> beta[prevsys:.trtclassIL-17 blocker]        0.46  4509    1
#> beta[bsa:.trtclassTNFa blocker]             1.31  3573    1
#> beta[bsa:.trtclassIL-12/23 blocker]         1.98  4347    1
#> beta[bsa:.trtclassIL-17 blocker]            1.29  3634    1
#> beta[weight:.trtclassTNFa blocker]         -0.10  3833    1
#> beta[weight:.trtclassIL-12/23 blocker]      0.00  5021    1
#> beta[weight:.trtclassIL-17 blocker]        -0.07  4493    1
#> beta[psa:.trtclassTNFa blocker]             0.35  4225    1
#> beta[psa:.trtclassIL-12/23 blocker]         0.75  6190    1
#> beta[psa:.trtclassIL-17 blocker]            0.43  4403    1
#> d[ETN]                                      1.81  2805    1
#> d[IXE_Q2W]                                  3.15  3031    1
#> d[IXE_Q4W]                                  2.95  2966    1
#> d[SEC_150]                                  2.47  2726    1
#> d[SEC_300]                                  2.88  2575    1
#> d[UST]                                      2.52  2667    1
#> lp__                                    -7635.39   899    1
#> tau                                         0.29   620    1
#> cc[PASI75]                                  0.00   NaN  NaN
#> cc[PASI90]                                  0.72  5665    1
#> cc[PASI100]                                 1.58  6212    1
#> 
#> Samples were drawn using NUTS(diag_e) at Sun Aug 28 14:09:25 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).

The estimated between-study heterogeneity standard deviation tau is small compared to the relative treatment effects. We compare the model fit using the DIC:

(pso_dic_FE <- dic(pso_fit_FE))

#> Residual deviance: 8811.4 (on 12387 data points)
#>                pD: 36
#>               DIC: 8847.4

(pso_dic_RE <- dic(pso_fit_RE))

#> Residual deviance: 8800.1 (on 12387 data points)
#>                pD: 42.3
#>               DIC: 8842.4

The DIC is lower for the RE model, indicating that there may be residual heterogeneity in the network and that the conditional constancy of relative effects assumption may be invalid—there may be additional effect modifiers that we have not accounted for. This result is different to the actual analysis reported by Phillippo et al. (2022), since here we are using synthetic IPD that have been simulated to closely resemble the original IPD. In the actual analysis the DIC was similar between the FE and RE models, so we might choose the more parsimonious FE model based on DIC alone, and there was no evidence for residual heterogeneity in this network.

Assessing residual inconsistency with an unrelated mean effects ML-NMR

We assess residual inconsistency using an unrelated mean effects model (Dias et al. 2011). Again, the call to the nma() function is identical, except this time we specify consistency = "ume". Node-splitting is also a possibility (with consistency = "nodesplit"), but this takes substantially longer since the model is re-run for each node-split comparison. Here we will proceed as in the analysis of Phillippo et al. (2022) and fit a fixed effect UME model (since there was no evidence for heterogeneity in the actual analysis); however, in our recreated analysis using synthetic IPD there was evidence of heterogeneity and we should really fit a random effects UME model instead.

pso_fit_UME <- nma(pso_net, 
                   trt_effects = "fixed",
                   consistency = "ume",
                   link = "probit", 
                   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_aux = flat(),
                   QR = TRUE,
                   init_r = 0.5)

#> Note: Setting "PBO" as the network reference treatment.

pso_fit_UME

#> A fixed effects ML-NMR with a ordered likelihood (probit link).
#> An inconsistency model ('ume') was fitted.
#> Regression model: ~(durnpso + prevsys + bsa + weight + psa) * .trt.
#> Centred covariates at the following overall mean values:
#>   durnpso   prevsys       bsa    weight       psa 
#> 1.8159830 0.6576489 0.2987820 8.9097263 0.2104826 
#> Inference for Stan model: ordered_multinomial.
#> 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.03    0.00 0.06    -0.09    -0.01     0.03     0.07
#> beta[prevsys]                              -0.17    0.00 0.16    -0.47    -0.27    -0.17    -0.06
#> beta[bsa]                                  -0.09    0.01 0.44    -0.98    -0.38    -0.08     0.21
#> beta[weight]                                0.04    0.00 0.03    -0.01     0.02     0.04     0.06
#> beta[psa]                                  -0.08    0.00 0.17    -0.42    -0.20    -0.08     0.03
#> beta[durnpso:.trtclassTNFa blocker]        -0.02    0.00 0.07    -0.16    -0.07    -0.02     0.03
#> beta[durnpso:.trtclassIL-12/23 blocker]    -0.06    0.00 0.10    -0.26    -0.13    -0.06     0.01
#> beta[durnpso:.trtclassIL-17 blocker]       -0.02    0.00 0.06    -0.14    -0.06    -0.02     0.02
#> beta[prevsys:.trtclassTNFa blocker]         0.19    0.00 0.18    -0.16     0.07     0.19     0.31
#> beta[prevsys:.trtclassIL-12/23 blocker]     0.45    0.01 0.35    -0.24     0.21     0.46     0.70
#> beta[prevsys:.trtclassIL-17 blocker]        0.16    0.00 0.16    -0.16     0.05     0.16     0.27
#> beta[bsa:.trtclassTNFa blocker]             0.22    0.01 0.50    -0.76    -0.12     0.20     0.56
#> beta[bsa:.trtclassIL-12/23 blocker]         0.59    0.01 0.67    -0.72     0.13     0.60     1.05
#> beta[bsa:.trtclassIL-17 blocker]            0.26    0.01 0.46    -0.64    -0.05     0.25     0.56
#> beta[weight:.trtclassTNFa blocker]         -0.16    0.00 0.03    -0.23    -0.18    -0.16    -0.14
#> beta[weight:.trtclassIL-12/23 blocker]     -0.09    0.00 0.05    -0.18    -0.12    -0.09    -0.06
#> beta[weight:.trtclassIL-17 blocker]        -0.13    0.00 0.03    -0.19    -0.15    -0.13    -0.11
#> beta[psa:.trtclassTNFa blocker]            -0.05    0.00 0.20    -0.44    -0.18    -0.05     0.09
#> beta[psa:.trtclassIL-12/23 blocker]         0.12    0.01 0.35    -0.54    -0.12     0.12     0.36
#> beta[psa:.trtclassIL-17 blocker]            0.10    0.00 0.18    -0.26    -0.03     0.09     0.22
#> d[ETN vs. PBO]                              1.58    0.00 0.07     1.44     1.53     1.58     1.63
#> d[IXE_Q2W vs. PBO]                          2.91    0.00 0.07     2.77     2.86     2.91     2.96
#> d[IXE_Q4W vs. PBO]                          2.69    0.00 0.07     2.54     2.64     2.69     2.74
#> d[SEC_150 vs. PBO]                          2.19    0.00 0.08     2.02     2.13     2.19     2.25
#> d[SEC_300 vs. PBO]                          2.60    0.00 0.08     2.44     2.54     2.60     2.65
#> d[UST vs. IXE_Q2W]                         -0.78    0.00 0.16    -1.09    -0.89    -0.78    -0.68
#> d[UST vs. SEC_300]                         -0.47    0.00 0.09    -0.65    -0.53    -0.47    -0.41
#> lp__                                    -7640.40    0.10 4.39 -7650.22 -7643.16 -7639.97 -7637.34
#> cc[PASI75]                                  0.00     NaN 0.00     0.00     0.00     0.00     0.00
#> cc[PASI90]                                  0.69    0.00 0.02     0.65     0.68     0.69     0.70
#> cc[PASI100]                                 1.53    0.00 0.02     1.49     1.52     1.53     1.55
#>                                            97.5% n_eff Rhat
#> beta[durnpso]                               0.15  2591    1
#> beta[prevsys]                               0.14  2698    1
#> beta[bsa]                                   0.74  2778    1
#> beta[weight]                                0.10  2854    1
#> beta[psa]                                   0.25  2997    1
#> beta[durnpso:.trtclassTNFa blocker]         0.12  2850    1
#> beta[durnpso:.trtclassIL-12/23 blocker]     0.13  3732    1
#> beta[durnpso:.trtclassIL-17 blocker]        0.11  3146    1
#> beta[prevsys:.trtclassTNFa blocker]         0.54  2803    1
#> beta[prevsys:.trtclassIL-12/23 blocker]     1.09  4034    1
#> beta[prevsys:.trtclassIL-17 blocker]        0.48  3191    1
#> beta[bsa:.trtclassTNFa blocker]             1.20  2846    1
#> beta[bsa:.trtclassIL-12/23 blocker]         1.91  3889    1
#> beta[bsa:.trtclassIL-17 blocker]            1.19  3282    1
#> beta[weight:.trtclassTNFa blocker]         -0.09  3169    1
#> beta[weight:.trtclassIL-12/23 blocker]      0.00  4218    1
#> beta[weight:.trtclassIL-17 blocker]        -0.07  3421    1
#> beta[psa:.trtclassTNFa blocker]             0.35  3036    1
#> beta[psa:.trtclassIL-12/23 blocker]         0.81  4680    1
#> beta[psa:.trtclassIL-17 blocker]            0.45  3454    1
#> d[ETN vs. PBO]                              1.73  2815    1
#> d[IXE_Q2W vs. PBO]                          3.05  2892    1
#> d[IXE_Q4W vs. PBO]                          2.84  3286    1
#> d[SEC_150 vs. PBO]                          2.36  2992    1
#> d[SEC_300 vs. PBO]                          2.76  3203    1
#> d[UST vs. IXE_Q2W]                         -0.47  4990    1
#> d[UST vs. SEC_300]                         -0.29  6387    1
#> lp__                                    -7633.00  1771    1
#> cc[PASI75]                                  0.00   NaN  NaN
#> cc[PASI90]                                  0.72  3712    1
#> cc[PASI100]                                 1.58  3257    1
#> 
#> Samples were drawn using NUTS(diag_e) at Sun Aug 28 14:26:48 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).

We compare model fit to the FE ML-NMR model using the DIC.

pso_dic_FE

#> Residual deviance: 8811.4 (on 12387 data points)
#>                pD: 36
#>               DIC: 8847.4

(pso_dic_UME <- dic(pso_fit_UME))

#> Residual deviance: 8811.7 (on 12387 data points)
#>                pD: 36.3
#>               DIC: 8848

The DIC values are similar between the FE model (assuming consistency) and the UME (inconsistency) model, which suggests no evidence for inconsistency overall. It is also important to compare the residual deviance contributions under each model to see whether there are any points that fit better under the UME model, as this can also indicate inconsistency. Using the plot() function produces a “dev-dev” plot of the residual deviance contributions under either model.

plot(pso_dic_FE, pso_dic_UME, show_uncertainty = FALSE) +
  xlab("Residual deviance - consistency model") +
  ylab("Residual deviance - inconsistency (UME) model")

plot of chunk pso-full-dev-dev

All points lie on the line of equality, so there is no evidence of inconsistency. If random effects models had been fitted then the heterogeneity estimates should also be compared as a drop in tau for the UME model can also indicate inconsistency.

Relaxing the shared effect modifier assumption

The treatment classes in the network are as follows:

data.frame(classes = pso_net$classes, treatments = pso_net$treatments)
#>            classes treatments
#> 1          Placebo        PBO
#> 2     TNFa blocker        ETN
#> 3    IL-17 blocker    IXE_Q2W
#> 4    IL-17 blocker    IXE_Q4W
#> 5    IL-17 blocker    SEC_150
#> 6    IL-17 blocker    SEC_300
#> 7 IL-12/23 blocker        UST

We fitted common interaction terms within each treatment class, under the shared effect modifier assumption, in order to make the model estimable with the available data. Note that only the interleukin-17 blocker class has more than one treatment; etanercept and ustekinumab are in classes of their own and so are unaffected by specifying class_interactions = "common". To assess this assumption we cannot simply fit independent interaction terms for all treatments and all effect modifiers at once as we do not have sufficient data. Instead, we relax this assumption one covariate at a time, estimating independent interactions for one covariate whilst keeping the shared effect modifier assumption (common interactions within each treatment class) for the other covariates.

To specify these relaxed models, we need to somehow mix class_interactions = "common" and class_interactions = "independent" for different covariates. The way we do this is with the .trt and .trtclass specials when specifying the regression model. To see how this works, first note that the model making the shared effect modifiers assumption

regression = ~(durnpso + prevsys + bsa + weight + psa)*.trt,
class_interactions = "common"

can be written equivalently using the .trtclass special as

regression = ~(durnpso + prevsys + bsa + weight + psa)*.trtclass

The .trtclass special is essentially a factor variable containing the treatment classes, and is available any time treatment classes have been specified in the network; this regression formula therefore has a single interaction term for each covariate within each treatment class (the same result as specifying class_interactions = "common" above). Finally, to fit independent interactions for a single covariate, say durnpso, we split these out using the .trt special with class_interactions = "independent" (i.e. telling the model not to combine interactions for .trt within classes):

regression = ~(prevsys + bsa + weight + psa)*.trtclass + durnpso*.trt,
class_interactions = "independent"

Since we are fitting several of these models, let us set up a list of model specifications and iterate over these.

noSEM_mods <- list(
  durnpso = ~(prevsys + bsa + weight + psa)*.trtclass + durnpso*.trt,
  prevsys = ~(durnpso + bsa + weight + psa)*.trtclass + prevsys*.trt,
  bsa = ~(durnpso + prevsys + weight + psa)*.trtclass + bsa*.trt,
  weight = ~(durnpso + prevsys + bsa + psa)*.trtclass + weight*.trt,
  psa = ~(durnpso + prevsys + bsa + weight)*.trtclass + psa*.trt
  )

noSEM_fits <- noSEM_mods

for (m in 1:length(noSEM_mods)) {
  cat("Fitting model with independent interactions for", names(noSEM_mods)[m], "\n")
  
  noSEM_fits[[m]] <- 
    nma(pso_net, 
        trt_effects = "fixed",
        link = "probit", 
        regression = noSEM_mods[[m]],
        class_interactions = "independent",
        prior_intercept = normal(scale = 10),
        prior_trt = normal(scale = 10),
        prior_reg = normal(scale = 10),
        prior_aux = flat(),
        QR = TRUE,
        init_r = 0.5,
        # Using save_warmup = FALSE reduces memory footprint when 
        # fitting many models in one session
        save_warmup = FALSE)
}

#> Fitting model with independent interactions for durnpso
#> Note: Setting "PBO" as the network reference treatment.
#> Fitting model with independent interactions for prevsys
#> Note: Setting "PBO" as the network reference treatment.
#> Fitting model with independent interactions for bsa
#> Note: Setting "PBO" as the network reference treatment.
#> Fitting model with independent interactions for weight
#> Note: Setting "PBO" as the network reference treatment.
#> Fitting model with independent interactions for psa
#> Note: Setting "PBO" as the network reference treatment.

Comparing model fit using the DIC

pso_dic_FE

#> Residual deviance: 8811.4 (on 12387 data points)
#>                pD: 36
#>               DIC: 8847.4

lapply(noSEM_fits, dic)

#> $durnpso
#> Residual deviance: 8812.5 (on 12387 data points)
#>                pD: 37.7
#>               DIC: 8850.3
#> 
#> $prevsys
#> Residual deviance: 8813.4 (on 12387 data points)
#>                pD: 37.8
#>               DIC: 8851.2
#> 
#> $bsa
#> Residual deviance: 8813 (on 12387 data points)
#>                pD: 37.8
#>               DIC: 8850.8
#> 
#> $weight
#> Residual deviance: 8807.5 (on 12387 data points)
#>                pD: 37.9
#>               DIC: 8845.4
#> 
#> $psa
#> Residual deviance: 8812.2 (on 12387 data points)
#>                pD: 37.8
#>               DIC: 8850

All of the models have similar or higher DIC to the original model making the shared effect modifier assumption for all covariates, with the only exception being the model with independent interactions for weight which has slightly lower DIC.

We also visually examine the differences between the estimated interaction terms under the original model (shared effect modifier assumption for all covariates) and the relaxed models (independent interactions, one covariate at a time).

library(purrr)
library(stringr)
library(forcats)

# Extract draws from relaxed models
imap_dfr(noSEM_fits,
        ~as_tibble(as.matrix(.x, pars = "beta")) %>% 
           pivot_longer(cols = everything(), names_to = "parameter", values_to = "value") %>% 
           filter(str_detect(parameter, paste0("(IXE|SEC).+:", .y))) %>% 
           mutate(model = .y)) %>% 
  
  # Add in draws from the original model
  bind_rows(
    as_tibble(as.matrix(pso_fit_FE, pars = "beta")) %>% 
    pivot_longer(cols = everything(), names_to = "parameter", values_to = "value") %>% 
    filter(str_detect(parameter, ":.+IL\\-17 blocker")) %>% 
    mutate(model = "all")
  ) %>% 
  
  mutate(
    # Rescale BSA to per 10% 
    value = if_else(str_detect(parameter, "bsa"), value / 10, value),
    # Create labels
    covariate = str_extract(parameter, "durnpso|prevsys|bsa|weight|psa"),
    covariatef = recode_factor(covariate,
                               durnpso = "Duration of psoriasis, per 10 years",
                               prevsys = "Previous systemic use",
                               bsa = "Body surface area, per 10%",
                               weight = "Weight, per 10 kg",
                               psa = "Psoriatic arthritis"),
    treatment = str_remove(str_extract(parameter, "\\.trt(class)?.+?(?=[\\]:])"),
                           "\\.trt(class)?"),
    Interactions = fct_collapse(factor(model), 
                                Common = "all", 
                                other_level = "Independent")) %>% 
  
# Plot
ggplot(aes(x = value, y = fct_rev(treatment), colour = Interactions, fill = Interactions)) +
  geom_vline(xintercept = 0, colour = "grey70") +
  ggdist::stat_halfeye(normalize = "panels", slab_alpha = 0.3, .width = c(0, 0.95)) +
  facet_wrap("covariatef", scales = "free") +
  xlab("Interaction effect (SMD)") + 
  ylab("Treatment / Class") +
  scale_colour_manual(values = c(Common = "#7B3294", Independent = "#91D388"),
                      aesthetics = c("colour", "fill")) +
  theme_multinma() +
  theme(legend.position = c(0.85, 0.2))

plot of chunk pso-full-relaxing-SEM

The independent interaction estimates are very similar to the common interaction estimates, but with much more uncertainty—particularly for the secukinumab regimens which are estimated only from aggregate data. The only exception is for weight, where there is some suggestion that this covariate may interact differently with the secukinumab treatment regimens to the ixekizumab regimens. However, the credible intervals for the secukinumab interactions are wide and overlap those for the ixekizumab regimens and the common interaction. Overall, there is some weak evidence that the shared effect modifier assumption (for the class of interleukin-17 blockers) may be invalid for weight. Since we are fitting multiple models here we should be mindful of multiple testing and the possibility that such differences have occurred by chance. On the other hand, this approach is likely to have low power to detect violations of the shared effect modifier assumption, particularly when the data are lacking. In this case, results from the model relaxing the shared effect modifier assumption for weight are very similar to the original model (see Phillippo et al. 2022).

Producing relative effects and event probabilities

Study populations included in the network

Population-average treatment effects can be produced for all the study populations represented in the network using the relative_effects() function.

(pso_releff_FE <- relative_effects(pso_fit_FE))

These relative effects can then be plotted using the plot() function.

plot(pso_releff_FE, ref_line = 0)

plot of chunk pso-full-releff

Similarly, average response probabilities on each treatment, in each study population, at each PASI cutoff can be produced using the predict() function. We specify type = "response" to produce predicted probabilities (rather than probit-probabilities).

(pso_pred_FE <- predict(pso_fit_FE, type = "response"))

Again, these can be plotted using the plot() function.

plot(pso_pred_FE, ref_line = c(0, 1))

plot of chunk pso-full-pred

External target populations

For the purposes of decision-making it is crucial that population-average estimates are produced for the decision target population of interest. The decision target population may not be represented by any of the study populations in the network, indeed it is likely best represented by an external registry or cohort study, or perhaps expert knowledge (Phillippo et al. 2016).

As an example, Phillippo et al. (2022) produce estimates for three external target populations represented by the PsoBest registry (Reich et al. 2015; Augustin et al. 2014), and the PROSPECT (Thaçi et al. 2019) and Chiricozzi 2019 (Chiricozzi et al. 2019) cohort studies. First of all, we need the covariate means and standard deviations in each of these populations:

new_agd_means <- tibble::tribble(
             ~study, ~covariate,  ~mean,   ~sd,
          "PsoBest",      "bsa",     24,  20.5,
          "PsoBest",  "durnpso",   18.2,  14.1,
          "PsoBest",  "prevsys",   0.54,    NA,
          "PsoBest",      "psa",  0.207,    NA,
          "PsoBest",   "weight",     85,  19.1,
         "PROSPECT",      "bsa",   18.7,  18.4,
         "PROSPECT",  "durnpso",   19.6,  13.5,
         "PROSPECT",  "prevsys", 0.9095,    NA,
         "PROSPECT",      "psa",  0.202,    NA,
         "PROSPECT",   "weight",   87.5,  20.3,
  "Chiricozzi 2019",      "bsa",     23, 16.79,
  "Chiricozzi 2019",  "durnpso",  16.93, 10.82,
  "Chiricozzi 2019",  "prevsys", 0.9061,    NA,
  "Chiricozzi 2019",      "psa", 0.2152,    NA,
  "Chiricozzi 2019",   "weight",   78.3, 15.87
  ) %>%
  # Tidy up
  pivot_wider(id_cols = study, 
              names_from = covariate, 
              values_from = c(mean, sd),
              names_glue = "{covariate}_{.value}") %>% 
  # Rescale as per analysis
  transmute(study,
            bsa_mean = bsa_mean / 100, 
            bsa_sd = bsa_sd / 100,
            weight_mean = weight_mean / 10,
            weight_sd = weight_sd / 10,
            durnpso_mean = durnpso_mean / 10,
            durnpso_sd = durnpso_sd / 10,
            prevsys = prevsys_mean,
            psa = psa_mean)

To produce estimates of population-average treatment effects, we use the relative_effects() function with the data frame of covariate means in the target populations as the newdata argument. We only need the covariate means, with variable names matching those in the regression.

(pso_releff_FE_new <- relative_effects(pso_fit_FE, 
                                       newdata = transmute(new_agd_means,
                                                           study,
                                                           bsa = bsa_mean,
                                                           weight = weight_mean,
                                                           durnpso = durnpso_mean,
                                                           prevsys,
                                                           psa),
                                       study = study))

#> -------------------------------------------------------- Study: Chiricozzi 2019 ---- 
#> 
#> Covariate values:
#>  durnpso prevsys  bsa weight  psa
#>     1.69    0.91 0.23   7.83 0.22
#> 
#>                             mean   sd 2.5%  25%  50%  75% 97.5% Bulk_ESS Tail_ESS Rhat
#> d[Chiricozzi 2019: ETN]     1.79 0.11 1.58 1.71 1.79 1.86  2.01     1618     2360    1
#> d[Chiricozzi 2019: IXE_Q2W] 3.07 0.10 2.87 3.00 3.07 3.14  3.28     1623     2630    1
#> d[Chiricozzi 2019: IXE_Q4W] 2.85 0.10 2.65 2.78 2.85 2.92  3.06     1782     2447    1
#> d[Chiricozzi 2019: SEC_150] 2.35 0.11 2.13 2.28 2.35 2.43  2.58     2248     2749    1
#> d[Chiricozzi 2019: SEC_300] 2.76 0.11 2.55 2.69 2.76 2.84  2.98     1883     2739    1
#> d[Chiricozzi 2019: UST]     2.30 0.15 2.02 2.20 2.30 2.40  2.58     2569     2878    1
#> 
#> --------------------------------------------------------------- Study: PROSPECT ---- 
#> 
#> Covariate values:
#>  durnpso prevsys  bsa weight psa
#>     1.96    0.91 0.19   8.75 0.2
#> 
#>                      mean   sd 2.5%  25%  50%  75% 97.5% Bulk_ESS Tail_ESS Rhat
#> d[PROSPECT: ETN]     1.62 0.11 1.42 1.55 1.62 1.70  1.84     2215     2573    1
#> d[PROSPECT: IXE_Q2W] 2.94 0.10 2.74 2.87 2.93 3.00  3.13     2128     2844    1
#> d[PROSPECT: IXE_Q4W] 2.72 0.10 2.52 2.65 2.72 2.78  2.92     2329     2929    1
#> d[PROSPECT: SEC_150] 2.22 0.11 2.00 2.14 2.22 2.29  2.44     2719     3321    1
#> d[PROSPECT: SEC_300] 2.63 0.11 2.41 2.55 2.63 2.70  2.85     2295     3224    1
#> d[PROSPECT: UST]     2.18 0.15 1.89 2.08 2.18 2.28  2.47     3205     3089    1
#> 
#> ---------------------------------------------------------------- Study: PsoBest ---- 
#> 
#> Covariate values:
#>  durnpso prevsys  bsa weight  psa
#>     1.82    0.54 0.24    8.5 0.21
#> 
#>                     mean   sd 2.5%  25%  50%  75% 97.5% Bulk_ESS Tail_ESS Rhat
#> d[PsoBest: ETN]     1.61 0.08 1.46 1.56 1.61 1.66  1.77     2021     2876    1
#> d[PsoBest: IXE_Q2W] 2.93 0.08 2.78 2.87 2.93 2.98  3.08     1857     2249    1
#> d[PsoBest: IXE_Q4W] 2.71 0.08 2.56 2.66 2.71 2.76  2.86     2164     2259    1
#> d[PsoBest: SEC_150] 2.21 0.09 2.04 2.15 2.21 2.27  2.39     2451     2706    1
#> d[PsoBest: SEC_300] 2.62 0.09 2.45 2.56 2.62 2.67  2.79     2016     2919    1
#> d[PsoBest: UST]     2.08 0.13 1.82 1.99 2.08 2.17  2.34     2852     3123    1
#>

These estimates are plotted using the plot() function.

plot(pso_releff_FE_new, ref_line = 0) + facet_wrap("Study")

plot of chunk pso-full-releff-new

Estimates of average event probabilities are produced by integrating predictions over the joint covariate distribution in each population. Since we have marginal summary statistics available, rather than full IPD, we create integration points using the add_integration() function by specifying the forms of the marginal distributions and the correlation matrix. We choose to use the same forms of the marginal distributions that we used when specifying integration points for the AgD studies in the network, and the weighted correlation matrix from the IPD studies.

new_agd_int <- add_integration(filter(new_agd_means, study != "PsoBest"),
                               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,
                               cor = pso_net$int_cor)

We then use the predict() function to produce the average event probabilities (type = "response", and level = "aggregate" which is the default) in each of the target populations. To do so, we also need to specify a distribution for the baseline event probabilities (i.e. probability of achieving PASI 75 response) in each of the target populations. PASI 75 event counts for individuals receiving secukinumab 300 mg treatment were available from PROSPECT (1156 achieved PASI 75 out of 1509) and Chiricozzi 2019 (243 out of 330), which we use to construct beta distributions on the baseline average response probabilities (we specify baseline_level = "aggregate" as these are population averages, rather than specific to a reference individual, and baseline_type = "response" as these are probabilities rather than transformed probit probabilities). No information on baseline response was available from PsoBest, so no predictions of absolute response rates could be made.

(pso_pred_FE_new <- predict(pso_fit_FE, 
        type = "response", 
        newdata = new_agd_int,
        study = study,
        baseline = list(PROSPECT = distr(qbeta, 1156, 1509-1156),
                        "Chiricozzi 2019" = distr(qbeta, 243, 330-243)),
        baseline_type = "response",
        baseline_level = "aggregate",
        trt_ref = "SEC_300"))

#> -------------------------------------------------------- Study: Chiricozzi 2019 ---- 
#> 
#>                                         mean   sd 2.5%  25%  50%  75% 97.5% Bulk_ESS Tail_ESS Rhat
#> pred[Chiricozzi 2019: PBO, PASI75]      0.02 0.01 0.01 0.01 0.02 0.02  0.03     2564     3595    1
#> pred[Chiricozzi 2019: PBO, PASI90]      0.00 0.00 0.00 0.00 0.00 0.00  0.01     2643     3414    1
#> pred[Chiricozzi 2019: PBO, PASI100]     0.00 0.00 0.00 0.00 0.00 0.00  0.00     2763     3690    1
#> pred[Chiricozzi 2019: ETN, PASI75]      0.37 0.04 0.30 0.34 0.37 0.40  0.46     4758     3892    1
#> pred[Chiricozzi 2019: ETN, PASI90]      0.16 0.03 0.11 0.14 0.16 0.18  0.22     4762     3997    1
#> pred[Chiricozzi 2019: ETN, PASI100]     0.03 0.01 0.02 0.03 0.03 0.04  0.05     4701     3931    1
#> pred[Chiricozzi 2019: IXE_Q2W, PASI75]  0.83 0.03 0.77 0.81 0.83 0.84  0.88     4481     3929    1
#> pred[Chiricozzi 2019: IXE_Q2W, PASI90]  0.60 0.04 0.52 0.57 0.60 0.63  0.68     4458     4057    1
#> pred[Chiricozzi 2019: IXE_Q2W, PASI100] 0.28 0.04 0.21 0.26 0.28 0.31  0.36     4418     4016    1
#> pred[Chiricozzi 2019: IXE_Q4W, PASI75]  0.76 0.03 0.70 0.74 0.77 0.79  0.83     4476     3773    1
#> pred[Chiricozzi 2019: IXE_Q4W, PASI90]  0.52 0.04 0.43 0.49 0.52 0.55  0.60     4444     3737    1
#> pred[Chiricozzi 2019: IXE_Q4W, PASI100] 0.22 0.03 0.16 0.19 0.21 0.24  0.28     4393     3701    1
#> pred[Chiricozzi 2019: SEC_150, PASI75]  0.59 0.04 0.52 0.57 0.59 0.62  0.66     4949     3965    1
#> pred[Chiricozzi 2019: SEC_150, PASI90]  0.33 0.03 0.26 0.30 0.33 0.35  0.40     4810     3887    1
#> pred[Chiricozzi 2019: SEC_150, PASI100] 0.10 0.02 0.07 0.09 0.10 0.11  0.14     4741     3846    1
#> pred[Chiricozzi 2019: SEC_300, PASI75]  0.74 0.02 0.69 0.72 0.74 0.75  0.78     3923     3931    1
#> pred[Chiricozzi 2019: SEC_300, PASI90]  0.48 0.03 0.42 0.46 0.48 0.50  0.54     3912     3974    1
#> pred[Chiricozzi 2019: SEC_300, PASI100] 0.19 0.02 0.15 0.17 0.19 0.20  0.23     3885     3933    1
#> pred[Chiricozzi 2019: UST, PASI75]      0.57 0.05 0.47 0.53 0.57 0.60  0.67     5538     3835    1
#> pred[Chiricozzi 2019: UST, PASI90]      0.31 0.05 0.22 0.28 0.31 0.34  0.41     5591     3798    1
#> pred[Chiricozzi 2019: UST, PASI100]     0.10 0.02 0.06 0.08 0.09 0.11  0.15     5462     3920    1
#> 
#> --------------------------------------------------------------- Study: PROSPECT ---- 
#> 
#>                                  mean   sd 2.5%  25%  50%  75% 97.5% Bulk_ESS Tail_ESS Rhat
#> pred[PROSPECT: PBO, PASI75]      0.03 0.01 0.02 0.03 0.03 0.04  0.05     2902     3561    1
#> pred[PROSPECT: PBO, PASI90]      0.01 0.00 0.00 0.00 0.01 0.01  0.01     3098     3528    1
#> pred[PROSPECT: PBO, PASI100]     0.00 0.00 0.00 0.00 0.00 0.00  0.00     3269     3415    1
#> pred[PROSPECT: ETN, PASI75]      0.40 0.03 0.33 0.38 0.40 0.42  0.47     5199     3240    1
#> pred[PROSPECT: ETN, PASI90]      0.18 0.02 0.14 0.16 0.18 0.20  0.23     5147     3365    1
#> pred[PROSPECT: ETN, PASI100]     0.04 0.01 0.03 0.04 0.04 0.05  0.06     4967     3490    1
#> pred[PROSPECT: IXE_Q2W, PASI75]  0.85 0.02 0.81 0.84 0.85 0.86  0.88     5001     3224    1
#> pred[PROSPECT: IXE_Q2W, PASI90]  0.64 0.03 0.57 0.62 0.64 0.66  0.70     4949     3615    1
#> pred[PROSPECT: IXE_Q2W, PASI100] 0.32 0.03 0.26 0.30 0.32 0.34  0.38     4883     3727    1
#> pred[PROSPECT: IXE_Q4W, PASI75]  0.79 0.02 0.74 0.78 0.79 0.81  0.84     5310     3627    1
#> pred[PROSPECT: IXE_Q4W, PASI90]  0.56 0.03 0.49 0.53 0.56 0.58  0.62     5152     3587    1
#> pred[PROSPECT: IXE_Q4W, PASI100] 0.25 0.03 0.19 0.23 0.24 0.26  0.30     4930     3584    1
#> pred[PROSPECT: SEC_150, PASI75]  0.63 0.03 0.58 0.61 0.63 0.65  0.68     5720     3577    1
#> pred[PROSPECT: SEC_150, PASI90]  0.36 0.03 0.31 0.35 0.36 0.38  0.42     5332     3598    1
#> pred[PROSPECT: SEC_150, PASI100] 0.12 0.01 0.09 0.11 0.12 0.13  0.15     5036     3554    1
#> pred[PROSPECT: SEC_300, PASI75]  0.77 0.01 0.74 0.76 0.77 0.77  0.79     3785     4056    1
#> pred[PROSPECT: SEC_300, PASI90]  0.52 0.02 0.49 0.51 0.52 0.53  0.55     3560     3888    1
#> pred[PROSPECT: SEC_300, PASI100] 0.22 0.01 0.19 0.21 0.22 0.23  0.24     3510     3811    1
#> pred[PROSPECT: UST, PASI75]      0.61 0.05 0.52 0.58 0.61 0.64  0.70     5992     3337    1
#> pred[PROSPECT: UST, PASI90]      0.35 0.04 0.27 0.32 0.35 0.38  0.44     6096     3471    1
#> pred[PROSPECT: UST, PASI100]     0.12 0.02 0.08 0.10 0.11 0.13  0.17     5841     3581    1

Again, we then plot these estimates using the plot() function, here with some customisation using ggplot syntax.

plot(pso_pred_FE_new, ref_line = c(0, 1)) + 
  facet_grid(rows = "Study") + 
  aes(colour = Category) +
  scale_colour_brewer(palette = "Blues")

plot of chunk pso-full-pred-new

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).↩︎
The alternative is “exclusive” format, where the lowest category counts those not achieving any higher outcomes (i.e. failure to achieve PASI 75, \(<75\%\) reduction in symptoms), the second counts those achieving PASI 75 but not PASI 90 or 100 (\(\ge 75\%\) and \(<90\%\) reduction), the third counts those achieving PASI 90 but not PASI 100 (\(\ge 90\%\) and \(<100\%\) reduction), and the final category counts those achieving PASI 100 (\(100\%\) reduction).↩︎

Example: Plaque psoriasis ML-NMR

Initial analysis

Setup

Preparing the data

Creating the network

Numerical integration for ML-NMR

ML-NMR models

Fixed effect ML-NMR

Random effects ML-NMR

Model comparison

Producing relative effects and event probabilities

Extended analysis

Setup

Preparing the data

Creating the network

Numerical integration for ML-NMR models

ML-NMR model

Assessing assumptions

Assessing residual heterogeneity with a random effects ML-NMR

Assessing residual inconsistency with an unrelated mean effects ML-NMR

Relaxing the shared effect modifier assumption

Producing relative effects and event probabilities

Study populations included in the network

External target populations

References