# Network meta-regression

Study-level covariates can be included in the model to adjust treatment effects following an approach for meta-regression outlined in NICE Technical Support Document 3 (Dias et al. 2011). This can be used to explore and account for potential effect modification.

Following the definition in NICE Technical Support Document 3, network meta-regression can be expressed as an interaction on the treatment effect in arms $$\geq2$$:

$\theta_{i,k}=\mu_i+(f(x,\beta_{a_{i,k}})-f(x,\beta_{a_{i,1}})) + (\psi_{1,a_{i,k}}-\psi_{1,a_{i,1}})$

where $$\theta_{i,k}$$ is the linear predictor, $$\mu_{i}$$ is the baseline effect on arm 1, $$f(x,\beta_{a_{i,k}})$$ is the dose-response function at dose $$x$$ with dose-response parameters $$\beta_{a_{i,k}}$$ for agent $$a$$ in arm $$k$$ of study $$i$$. $$\psi_{1,a_{i,k}}$$ is then the effect modifying interaction between the agent in arm $$k$$ and the network reference agent (typically Placebo in a dose-response analysis).

## Data preparation

To improve estimation:

• continuous covariates should be centred around their mean
• binary/categorical variables should be recoded with the most commonly reported value as the reference category
# Using the SSRI dataset
ssri.reg <- ssri

# For a continuous covariate
ssri.reg <- ssri.reg %>%
dplyr::mutate(x.weeks = weeks - mean(weeks, na.rm = TRUE))

# For a categorical covariate
table(ssri\$weeks)  # Using 8 weeks as the reference
ssri.reg <- ssri.reg %>%
dplyr::mutate(r.weeks = factor(weeks, levels = c(8, 4, 5, 6, 9, 10)))

# Create network object
ssrinet <- mbnma.network(ssri.reg)
#> Values for agent with dose = 0 have been recoded to Placebo
#> agent is being recoded to enforce sequential numbering

## Modelling

For performing network meta-regression, different assumptions can be made regarding how the effect modification may be shared across agents:

### Independent, agent-specific interactions

The least constraining assumption available in MBNMAdose is to assume that the effect modifier acts on each agent independently, and separate $$\psi_{1,a_{i,k}}$$ are therefore estimated for each agent in the network.

A slightly stronger assumption is to assume that agents within the same class share the same interaction effect, though classes must be specified within the dataset for this.

# Regress for continuous weeks Separate effect modification for each agent vs
# Placebo
ssrimod.a <- mbnma.run(ssrinet, fun = dfpoly(degree = 2), regress = ~x.weeks, regress.effect = "agent")
summary(ssrimod.a)
#> ========================================
#> Dose-response MBNMA
#> ========================================
#>
#> Likelihood: binomial
#> Dose-response function: fpoly
#>
#> Pooling method
#>
#> Method: Common (fixed) effects estimated for relative effects
#>
#>
#> beta.1 dose-response parameter results
#>
#> Pooling: relative effects for each agent
#>
#> |Agent        |Parameter |  Median|    2.5%|   97.5%|
#> |:------------|:---------|-------:|-------:|-------:|
#> |citalopram   |beta.1[2] | -0.0962| -0.3516|  0.1664|
#> |escitalopram |beta.1[3] |  0.1368| -0.1349|  0.4014|
#> |fluoxetine   |beta.1[4] |  0.2699|  0.0342|  0.5111|
#> |paroxetine   |beta.1[5] | -0.1387| -0.4583|  0.1486|
#> |sertraline   |beta.1[6] |  0.9668| -8.1137| 10.3439|
#>
#>
#> beta.2 dose-response parameter results
#>
#> Pooling: relative effects for each agent
#>
#> |Agent        |Parameter |  Median|    2.5%|  97.5%|
#> |:------------|:---------|-------:|-------:|------:|
#> |citalopram   |beta.2[2] |  0.0672| -0.0040| 0.1389|
#> |escitalopram |beta.2[3] | -0.0090| -0.0950| 0.0765|
#> |fluoxetine   |beta.2[4] | -0.0663| -0.1349| 0.0020|
#> |paroxetine   |beta.2[5] |  0.0913|  0.0001| 0.1867|
#> |sertraline   |beta.2[6] | -0.1432| -1.1768| 0.8602|
#>
#>
#> power.1 dose-response parameter results
#>
#> Assigned a numeric value: 0
#>
#> power.2 dose-response parameter results
#>
#> Assigned a numeric value: 0
#>
#> Meta-regression
#>
#> Covariates interacting with study-level relative effects: x.weeks
#> Common (identical) covariate-by-agent effects
#>
#>
#> |Agent        |Parameter    |  Median|     2.5%|   97.5%|
#> |:------------|:------------|-------:|--------:|-------:|
#> |citalopram   |B.x.weeks[2] | -0.1596|  -0.3310|  0.0118|
#> |escitalopram |B.x.weeks[3] |  0.1598|  -0.0296|  0.3567|
#> |fluoxetine   |B.x.weeks[4] | -0.0243|  -0.1731|  0.1164|
#> |paroxetine   |B.x.weeks[5] |  0.0868|  -0.0218|  0.1968|
#> |sertraline   |B.x.weeks[6] |  0.8143| -18.3105| 20.2817|
#>
#>
#> Model Fit Statistics
#> Effective number of parameters:
#> pD calculated using the Kullback-Leibler divergence = 74.8
#>
#> Deviance = 890.6
#> Residual deviance = 190.9
#> Deviance Information Criterion (DIC) = 965.3

Within the output, a separate parameter (named B.x.weeks[]) has been estimated for each agent that corresponds to the effect of an additional week of study follow-up on the relative effect of the agent versus Placebo. Note that due to the inclusion of weeks as a continuous covariate, we are assuming a linear effect modification due to study follow-up.

### Random effect interaction

Alternatively, the effect modification for different agents versus the network reference agent can be assumed to be exchangeable/shared across the network about a common mean, $$\hat{\psi}$$, with a between-agent standard deviation of $$\tau_\psi$$:

$\psi_{1,a_{i,k}} \sim N(\hat{\psi}, \tau^2_\psi)$

# Regress for continuous weeks Random effect modification across all agents vs
# Placebo
ssrimod.r <- mbnma.run(ssrinet, fun = dfpoly(degree = 2), regress = ~x.weeks, regress.effect = "random")
summary(ssrimod.r)
#> ========================================
#> Dose-response MBNMA
#> ========================================
#>
#> Likelihood: binomial
#> Dose-response function: fpoly
#>
#> Pooling method
#>
#> Method: Common (fixed) effects estimated for relative effects
#>
#>
#> beta.1 dose-response parameter results
#>
#> Pooling: relative effects for each agent
#>
#> |Agent        |Parameter |  Median|    2.5%|  97.5%|
#> |:------------|:---------|-------:|-------:|------:|
#> |citalopram   |beta.1[2] | -0.0954| -0.4818| 0.2033|
#> |escitalopram |beta.1[3] |  0.2466| -0.2042| 0.6454|
#> |fluoxetine   |beta.1[4] |  0.2833| -0.0172| 0.5275|
#> |paroxetine   |beta.1[5] | -0.1267| -0.6161| 0.1569|
#> |sertraline   |beta.1[6] |  0.6682|  0.1157| 1.2429|
#>
#>
#> beta.2 dose-response parameter results
#>
#> Pooling: relative effects for each agent
#>
#> |Agent        |Parameter |  Median|    2.5%|   97.5%|
#> |:------------|:---------|-------:|-------:|-------:|
#> |citalopram   |beta.2[2] |  0.0707| -0.0124|  0.1774|
#> |escitalopram |beta.2[3] | -0.0365| -0.1729|  0.1071|
#> |fluoxetine   |beta.2[4] | -0.0683| -0.1416|  0.0197|
#> |paroxetine   |beta.2[5] |  0.0924|  0.0005|  0.2472|
#> |sertraline   |beta.2[6] | -0.1141| -0.2355| -0.0018|
#>
#>
#> power.1 dose-response parameter results
#>
#> Assigned a numeric value: 0
#>
#> power.2 dose-response parameter results
#>
#> Assigned a numeric value: 0
#>
#> Meta-regression
#>
#> Covariates interacting with study-level relative effects: x.weeks
#> Random (exchangeable) covariate-by-treatment effects
#>
#>
#> |Regression effect |Parameter | Median|    2.5%| 97.5%|
#> |:-----------------|:---------|------:|-------:|-----:|
#> |Random effect     |B.x.weeks | 0.0288| -0.0667| 0.133|
#>
#>
#> Standard deviation for random covariate-by-treatment effects
#>
#> |Parameter    | Median|   2.5%|  97.5%|
#> |:------------|------:|------:|------:|
#> |sd.B.x.weeks | 0.0661| 0.0014| 0.3085|
#>
#>
#> Model Fit Statistics
#> Effective number of parameters:
#> pD calculated using the Kullback-Leibler divergence = 73
#>
#> Deviance = 890
#> Residual deviance = 190.3
#> Deviance Information Criterion (DIC) = 962.9

In this case only a single regression paramter is estimated (B.x.weeks), which corresponds to the mean effect of an additional week of study follow-up on the relative effect of an active agent versus Placebo. A parameter is also estimated for the between-agent standard deviation, sd.B.x.weeks.

### Common effect interaction

This is the strongest assumption for network meta-regression, and it implies that effect modification is common (equal) for all agents versus the network reference agent:

$\psi_{1,a_{i,k}} =\hat{\psi}$

# Regress for categorical weeks Common effect modification across all agents vs
# Placebo
ssrimod.c <- mbnma.run(ssrinet, fun = dfpoly(degree = 2), regress = ~r.weeks, regress.effect = "common")
summary(ssrimod.c)
#> ========================================
#> Dose-response MBNMA
#> ========================================
#>
#> Likelihood: binomial
#> Dose-response function: fpoly
#>
#> Pooling method
#>
#> Method: Common (fixed) effects estimated for relative effects
#>
#>
#> beta.1 dose-response parameter results
#>
#> Pooling: relative effects for each agent
#>
#> |Agent        |Parameter |  Median|    2.5%|  97.5%|
#> |:------------|:---------|-------:|-------:|------:|
#> |citalopram   |beta.1[2] | -0.0558| -0.3285| 0.2058|
#> |escitalopram |beta.1[3] |  0.2446|  0.0133| 0.4764|
#> |fluoxetine   |beta.1[4] |  0.2235| -0.0197| 0.4802|
#> |paroxetine   |beta.1[5] | -0.1064| -0.4060| 0.1794|
#> |sertraline   |beta.1[6] |  0.5771|  0.1003| 1.0505|
#>
#>
#> beta.2 dose-response parameter results
#>
#> Pooling: relative effects for each agent
#>
#> |Agent        |Parameter |  Median|    2.5%|   97.5%|
#> |:------------|:---------|-------:|-------:|-------:|
#> |citalopram   |beta.2[2] |  0.0574| -0.0130|  0.1286|
#> |escitalopram |beta.2[3] | -0.0289| -0.1074|  0.0499|
#> |fluoxetine   |beta.2[4] | -0.0572| -0.1279|  0.0090|
#> |paroxetine   |beta.2[5] |  0.0842| -0.0070|  0.1794|
#> |sertraline   |beta.2[6] | -0.0994| -0.1959| -0.0032|
#>
#>
#> power.1 dose-response parameter results
#>
#> Assigned a numeric value: 0
#>
#> power.2 dose-response parameter results
#>
#> Assigned a numeric value: 0
#>
#> Meta-regression
#>
#> Covariates interacting with study-level relative effects: r.weeks4, r.weeks5, r.weeks6, r.weeks9, r.weeks10
#> Common (identical) covariate-by-treatment effects
#>
#>
#> |Regression effect |Parameter   |  Median|    2.5%|   97.5%|
#> |:-----------------|:-----------|-------:|-------:|-------:|
#> |Common effect     |B.r.weeks10 |  0.1555| -0.2480|  0.5409|
#> |Common effect     |B.r.weeks4  | -0.7274| -1.4587| -0.0282|
#> |Common effect     |B.r.weeks5  |  0.8392| -0.1774|  1.8885|
#> |Common effect     |B.r.weeks6  |  0.0505| -0.1363|  0.2497|
#> |Common effect     |B.r.weeks9  |  0.3967| -0.4041|  1.2484|
#>
#>
#> Model Fit Statistics
#> Effective number of parameters:
#> pD calculated using the Kullback-Leibler divergence = 74.7
#>
#> Deviance = 891.3
#> Residual deviance = 191.5
#> Deviance Information Criterion (DIC) = 966

In this case we have performed the network meta-regression on study follow-up (weeks) as a categorical covariate. Therefore, although only a single parameter is estimated for each effect modifying term, there is a separate term for each category of week and a linear relationship for effect modification is no longer assumed.

### Alternative assumptions

Although this is beyond the capability of MBNMAdose, one could envision a more complex model in which the interaction effect also varied by a dose-response relationship, rather than assuming an effect by agent/class or across the whole network. This would in principle contain fewer parameters than a fully independent interaction model (in which a separate regression covariate is estimated for each treatment in the dataset).

### Aggregation bias

Note that adjusting for aggregated patient-level covariates (e.g. mean age, % males, etc.) whilst using a non-identity link function can introduce aggregation bias. This is a form of ecological bias that biases treatment effects towards the null and is typically more severe where treatment effects are strong and where the link function is highly non-linear (Dias et al. 2011). This can be resolved by performing a patient-level regression, but Individual Participant Data are required for this and such an analysis is outside the scope of MBNMAdose.

## Prediction using effect modifying covariates

Models fitted with meta-regression can also be used to make predictions for a specified set of covariate values. This includes when estimating relative effects using get.relative(). An additional argument regress.vals can be used to provide a named vector of covariate values at which to make predictions.

# For a continuous covariate, make predictions at 5 weeks follow-up
pred <- predict(ssrimod.a, regress.vals = c(x.weeks = 5))
plot(pred)

Predictions are very uncertain for Sertraline, as studies only investigated this agent at 6 weeks follow-up and therefore the agent-specific effect modification is very poorly estimated.

# For a categorical covariate, make predictions at 10 weeks follow-up
regress.p <- c(r.weeks10 = 1, r.weeks4 = 0, r.weeks5 = 0, r.weeks6 = 0, r.weeks9 = 0)

pred <- predict(ssrimod.c, regress.vals = regress.p)
plot(pred)

Note that categorical covariates are modelled as multiple binary dummy covariates, and so a value for each of these must be included.

## References

Dias, S., A. J. Sutton, N. J. Welton, and A. E. Ades. 2011. “NICE DSU Technical Support Document 3: Heterogeneity: Subgroups, Meta-Regression, Bias and Bias, Adjustment.” Decision Support Unit, ScHARR, University of Sheffield. https://www.ncbi.nlm.nih.gov/books/NBK395886/pdf/Bookshelf_NBK395886.pdf.