jarbes (version 2.2.5)

metarisk: Bayesian meta-analysis for using baseline risk adjustment

Description

This function performers a Bayesian meta-analysis to analyse heterogeneity of the treatment effect as a function of the baseline risk. The function fits a hierarchical meta-regression model based on a bivariate random effects model.

Usage

metarisk(
  data,
  two.by.two = TRUE,
  re = "normal",
  link = "logit",
  mean.mu.1 = 0,
  mean.mu.2 = 0,
  sd.mu.1 = 1,
  sd.mu.2 = 1,
  sigma.1.upper = 5,
  sigma.2.upper = 5,
  mean.Fisher.rho = 0,
  sd.Fisher.rho = 1/sqrt(2),
  df = 4,
  df.estimate = FALSE,
  df.lower = 3,
  df.upper = 20,
  split.w = FALSE,
  nr.chains = 2,
  nr.iterations = 10000,
  nr.adapt = 1000,
  nr.burnin = 1000,
  nr.thin = 1
)

Value

This function returns an object of the class "metarisk". This object contains the MCMC output of each parameter and hyper-parameter in the model, the data frame used for fitting the model, the link function, type of random effects distribution and the splitting information for conflict of evidence analysis.

The results of the object of the class metadiag can be extracted with R2jags or with rjags. In addition a summary, a print and a plot functions are implemented for this type of object.

Arguments

data

A data frame where the first four columns containing the number of events in the control group (yc), the number of patients in the control group (nc), the number of events in the treatment group (yt) and the number of patients in the treatment group (nt). If two.by.two = TRUE a data frame where each line contains the trial results with column names: yc, nc, yt, nt.

two.by.two

If TRUE indicates that the trial results are with names: yc, nc, yt, nt.

re

Random effects distribution for the resulting model. Possible values are normal for bivariate random effects and sm for scale mixtures.

link

The link function used in the model. Possible values are logit, cloglog probit.

mean.mu.1

Prior mean of baseline risk, default value is 0.

mean.mu.2

Prior mean of the relative treatment effect, default value is 0.

sd.mu.1

Prior standard deviation of mu.1, default value is 1. The default prior of mu.1 is a logistic distribution with mean 0 and dispersion 1. The implicit prior for mu.1 in the probability scale is a uniform between 0 and 1.

sd.mu.2

Prior standard deviation of mu.2, default value is 1. The default prior of mu.2 is a logistic distribution with mean 0 and dispersion 1. The implicit prior for mu.2 in the probability scale is a uniform between 0 and 1.

sigma.1.upper

Upper bound of the uniform prior of sigma.1, default value is 5.

sigma.2.upper

Upper bound of the uniform prior of sigma.2, default value is 5.

mean.Fisher.rho

Mean of rho in the Fisher scale default value is 0.

sd.Fisher.rho

Standard deviation of rho in the Fisher scale, default value is 1/sqrt(2).

df

If de.estimate = FALSE, then df is the degrees of freedom for the scale mixture distribution, default value is 4.

df.estimate

Estimate the posterior of df. The default value is FALSE.

df.lower

Lower bound of the prior of df. The default value is 3.

df.upper

Upper bound of the prior of df. The default value is 30.

split.w

Split the w parameter in two independent weights one for each random effect. The default value is FALSE.

nr.chains

Number of chains for the MCMC computations, default 5.

nr.iterations

Number of iterations after adapting the MCMC, default is 10000. Some models may need more iterations.

nr.adapt

Number of iterations in the adaptation process, defualt is 1000. Some models may need more iterations during adptation.

nr.burnin

Number of iteration discared for burnin period, default is 1000. Some models may need a longer burnin period.

nr.thin

Thinning rate, it must be a positive integer, the default value is 1.

Details

The number of events in the control and treated group are modeled with two conditional Binomial distributions and the random-effects are based on a bivariate scale mixture of Normals.

The function calculates the implicit hierarchical meta-regression, where the treatment effect is regressed to the baseline risk (rate of events in the control group). The scale mixture weights are used to adjust for internal validity and structural outliers identification.

Computations are done by calling JAGS (Just Another Gibbs Sampler) to perform MCMC (Markov Chain Monte Carlo) sampling and returning an object of the class mcmc.list.

References

Verde, P.E. (2019) The hierarchical meta-regression approach and learning from clinical evidence. Biometrical Journal. 1 - 23.

Verde, P. E. (2017) Two Examples of Bayesian Evidence Synthesis with the Hierarchical Meta-Regression Approach. Chap.9, pag 189-206. Bayesian Inference, ed. Tejedor, Javier Prieto. InTech.

Examples

if (FALSE) {
library(jarbes)

# This example is used to test the function and it runs in about 5 seconds.
# In a real application you must increase the number of MCMC interations.
# For example use: nr.burnin = 5000 and nr.iterations = 20000





# The following examples corresponds to Section 4 in Verde (2019).
# These are simulated examples to investigate internal and
# external validity bias in meta-analysis.


library(MASS)
library(ggplot2)
library(jarbes)
library(gridExtra)
library(mcmcplots)

#Experiment 1: External validity bias

set.seed(2018)
# mean control
pc <- 0.7
# mean treatment
pt <- 0.4
# reduction of "amputations" odds ratio
OR <- (pt/(1-pt)) /(pc/(1-pc))
OR

# mu_2
log(OR)
mu.2.true <- log(OR)
#sigma_2
sigma.2.true <- 0.5
# mu_1
mu.1.true <- log(pc/(1-pc))
mu.1.true
#sigma_1
sigma.1.true <- 1
# rho
rho.true <- -0.5
Sigma <- matrix(c(sigma.1.true^2, sigma.1.true*sigma.2.true*rho.true,
                 sigma.1.true*sigma.2.true*rho.true, sigma.2.true^2),
                 byrow = TRUE, ncol = 2)
Sigma

theta <- mvrnorm(35, mu = c(mu.1.true, mu.2.true),
                Sigma = Sigma )


plot(theta, xlim = c(-2,3))
abline(v=mu.1.true, lty = 2)
abline(h=mu.2.true, lty = 2)
abline(a = mu.2.true, b=sigma.2.true/sigma.1.true * rho.true, col = "red")
abline(lm(theta[,2]~theta[,1]), col = "blue")

# Target group
mu.T <- mu.1.true + 2 * sigma.1.true
abline(v=mu.T, lwd = 3, col = "blue")
eta.true <- mu.2.true + sigma.2.true/sigma.1.true*rho.true* mu.T
eta.true
exp(eta.true)
abline(h = eta.true, lty =3, col = "blue")
# Simulation of each primary study:
n.c <- round(runif(35, min = 30, max = 60),0)
n.t <- round(runif(35, min = 30, max = 60),0)
y.c <- y.t <- rep(0, 35)
p.c <- exp(theta[,1])/(1+exp(theta[,1]))
p.t <- exp(theta[,2]+theta[,1])/(1+exp(theta[,2]+theta[,1]))
for(i in 1:35)
{
 y.c[i] <- rbinom(1, n.c[i], prob = p.c[i])
 y.t[i] <- rbinom(1, n.t[i], prob = p.t[i])
}

AD.s1 <- data.frame(yc=y.c, nc=n.c, yt=y.t, nt=n.t)

##########################################################
incr.e <- 0.05
incr.c <- 0.05
n11 <- AD.s1$yt
n12 <- AD.s1$yc
n21 <- AD.s1$nt - AD.s1$yt
n22 <- AD.s1$nc - AD.s1$yc
AD.s1$TE <- log(((n11 + incr.e)*(n22 + incr.c))/((n12 + incr.e) * (n21 + incr.c)))
AD.s1$seTE <- sqrt((1/(n11 + incr.e) + 1/(n12 + incr.e) +
                     1/(n21 + incr.c) + 1/(n22 + incr.c)))

Pc <- ((n12 + incr.c)/(AD.s1$nc + 2*incr.c))

AD.s1$logitPc <- log(Pc/(1-Pc))

AD.s1$Ntotal <- AD.s1$nc + AD.s1$nt
rm(list=c("Pc", "n11","n12","n21","n22","incr.c", "incr.e"))


dat.points <- data.frame(TE = AD.s1$TE, logitPc = AD.s1$logitPc, N.total = AD.s1$Ntotal)
###################################################################

res.s1 <- metarisk(AD.s1, two.by.two = FALSE, sigma.1.upper = 5,
                  sigma.2.upper = 5,
                  sd.Fisher.rho = 1.5)

print(res.s1, digits = 4)


library(R2jags)
attach.jags(res.s1)
eta.hat <- mu.2 + rho*sigma.2/sigma.1*(mu.T - mu.1)
mean(eta.hat)
sd(eta.hat)

OR.eta.hat <- exp(eta.hat)
mean(OR.eta.hat)
sd(OR.eta.hat)
quantile(OR.eta.hat, probs = c(0.025, 0.5, 0.975))

ind.random <- sample(1:18000, size = 80, replace = FALSE)

##############################################################
p1 <- ggplot(dat.points, aes(x = logitPc, y = TE, size = N.total)) +
      xlab("logit Baseline Risk")+
      ylab("log(Odds Ratio)")+
      geom_point(shape = 21, colour = "blue") + scale_size_area(max_size=10)+
      scale_x_continuous(name= "Rate of The Control Group (logit scale)",
                       limits=c(-2, 5)) +
     geom_vline(xintercept = mu.T, colour = "blue", size = 1, lty = 1) +
       geom_hline(yintercept = eta.true, colour = "blue", size = 1, lty = 1)+
         geom_abline(intercept=beta.0[ind.random],
                   slope=beta.1[ind.random],alpha=0.3,
                   colour = "grey", size = 1.3, lty = 2)+
         geom_abline(intercept = mean(beta.0[ind.random]),
         slope = mean(beta.1[ind.random]),
         colour = "black", size = 1.3, lty = 2)+
     geom_abline(intercept = mu.2.true, slope = sigma.2.true/sigma.1.true * rho.true,
     colour = "blue", size = 1.2)+ theme_bw()


# Posterior of eta.hat

eta.df <- data.frame(x = OR.eta.hat)


p2 <- ggplot(eta.df, aes(x = x)) +
 xlab("Odds Ratio") +
 ylab("Posterior distribution")+
 geom_histogram(fill = "royalblue", colour = "black", alpha= 0.4, bins=80) +
 geom_vline(xintercept = exp(eta.true), colour = "black", size = 1.7, lty = 1)+
 geom_vline(xintercept = c(0.28, 0.22, 0.35), colour = "black", size = 1, lty = 2)+
 theme_bw()


grid.arrange(p1, p2, ncol = 2, nrow = 1)

#Experiment 2: Internal validity bias and assesing conflict of evidence between the RCTs.


set.seed(2018)
ind <- sample(1:35, size=5, replace = FALSE)
ind
AD.s4.contaminated <- AD.s1[ind,1:4]
AD.s4.contaminated$yc <- AD.s1$yt[ind]
AD.s4.contaminated$yt <- AD.s1$yc[ind]
AD.s4.contaminated$nc <- AD.s1$nt[ind]
AD.s4.contaminated$nt <- AD.s1$nc[ind]
AD.s4.contaminated <- rbind(AD.s4.contaminated,
                         AD.s1[-ind,1:4])

res.s4 <- metarisk(AD.s4.contaminated,
                   two.by.two = FALSE,
                   re = "sm",
                   sigma.1.upper = 3,
                   sigma.2.upper = 3,
                   sd.Fisher.rho = 1.5,
                   df.estimate = TRUE)

print(res.s4, digits = 4)

attach.jags(res.s4)

w.s <- apply(w, 2, median)
w.u <- apply(w, 2, quantile, prob = 0.75)
w.l <- apply(w, 2, quantile, prob = 0.25)

studies <- c(ind,c(1,3,4,5,6,8,9,10,11,13,14,17,18,19,20:35))


dat.weights <- data.frame(x = studies,
                         y = w.s,
                        ylo  = w.l,
                        yhi  = w.u)

# Outliers:
w.col <- studies %in% ind
w.col.plot <- ifelse(w.col, "black", "grey")
w.col.plot[c(9,17, 19,27,34,35)] <- "black"

w.plot <- function(d){
  # d is a data frame with 4 columns
  # d$x gives variable names
  # d$y gives center point
  # d$ylo gives lower limits
  # d$yhi gives upper limits

  p <- ggplot(d, aes(x=x, y=y, ymin=ylo, ymax=yhi) )+
       geom_pointrange( colour=w.col.plot, lwd =0.8)+
       coord_flip() + geom_hline(yintercept = 1, lty=2)+
       xlab("Study ID") +
       ylab("Scale mixture weights") + theme_bw()
       return(p)}

w.plot(dat.weights)

#List of other possible statistical models:
#    1) Different link functions: logit, cloglog and probit

#    2) Different random effects distributions:
#       "normal" or "sm = scale mixtures".

#    3) For the scale mixture random effects:
#       split.w = TRUE => "split the weights".

#    4) For the scale mixture random effects:
#       df.estimate = TRUE => "estimate the degrees of freedom".

#    5) For the scale mixture random effects:
#       df.estimate = TRUE => "estimate the degrees of freedom".

#    6) For the scale mixture random effects:
#       df = 4 => "fix the degrees of freedom to a particual value".
#       Note that df = 1 fits a Cauchy bivariate distribution to
#       the random effects.
#End of the examples

}