library(gemtc)
## Loading required package: coda
# Load the example network and generate a consistency model:
model <- mtc.model(smoking, type="consistency")
# Load pre-generated samples instead of runing the model:
results <- mtc.run(model, thin=10)
## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 50
## Unobserved stochastic nodes: 54
## Total graph size: 1705
##
## Initializing model
results <- dget(system.file("extdata/luades-smoking.samples.gz", package="gemtc"))
# Print a basic statistical summary of the results:
summary(results)
##
## Results on the Log Odds Ratio scale
##
## Iterations = 5010:25000
## Thinning interval = 10
## Number of chains = 4
## Sample size per chain = 2000
##
## 1. Empirical mean and standard deviation for each variable,
## plus standard error of the mean:
##
## Mean SD Naive SE Time-series SE
## d.A.B 0.4965 0.4081 0.004563 0.004746
## d.A.C 0.8359 0.2433 0.002720 0.003022
## d.A.D 1.1088 0.4355 0.004869 0.005256
## sd.d 0.8465 0.1913 0.002139 0.003034
##
## 2. Quantiles for each variable:
##
## 2.5% 25% 50% 75% 97.5%
## d.A.B -0.2985 0.2312 0.4910 0.7530 1.341
## d.A.C 0.3878 0.6720 0.8273 0.9867 1.353
## d.A.D 0.2692 0.8197 1.0983 1.3824 2.006
## sd.d 0.5509 0.7119 0.8180 0.9542 1.283
library(gemtc)
# Build a model similar to Model 4(b) from Cooper et al. (2009):
classes <- list("control"=c("01"),
"anti-coagulant"=c("02","03","04","09"),
"anti-platelet"=c("05","06","07","08","10","11","12","16","17"),
"mixed"=c("13","14","15"))
regressor <- list(coefficient='shared',
variable='stroke',
classes=classes)
model <- mtc.model(atrialFibrillation,
type="regression",
regressor=regressor,
om.scale=10)
model
## MTC regression model:
data <- read.table(textConnection('
id group pe ci.l ci.u style value.A value.B
"Study 1" 1 0.35 0.08 0.92 "normal" "2/46" "7/46"
"Study 2" 1 0.43 0.15 1.14 "normal" "4/50" "8/49"
"Study 3" 2 0.31 0.07 0.74 "normal" "2/97" "10/100"
"Study 4" 2 0.86 0.34 2.90 "normal" "9/104" "6/105"
"Study 5" 2 0.33 0.10 0.72 "normal" "4/74" "14/74"
"Study 6" 2 0.47 0.23 0.91 "normal" "11/120" "22/129"
"Pooled" NA 0.42 0.15 1.04 "pooled" NA NA
'), header=TRUE)
data$pe <- log(data$pe)
data$ci.l <- log(data$ci.l)
data$ci.u <- log(data$ci.u)
blobbogram(data, group.labels=c('GROUP 1', 'GROUP 2'),
columns=c('value.A', 'value.B'), column.labels=c('r/n', 'r/n'),
column.groups=c(1, 2), grouped=TRUE,
column.group.labels=c('Intervention', 'Control'),
id.label="Trial", ci.label="Odds Ratio (95% CrI)", log.scale=TRUE)
library(gemtc)
# Run RE regression model with informative heterogeneity prior
regressor <- list(coefficient='shared',
variable='diseaseDuration',
control='Placebo')
# sd ~ half-Normal(mean=0, sd=0.32)
hy.prior <- mtc.hy.prior(type="std.dev", distr="dhnorm", 0, 9.77)
model <- mtc.model(certolizumab,
type="regression",
regressor=regressor,
hy.prior=hy.prior)
result <- mtc.run(model)
## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 24
## Unobserved stochastic nodes: 32
## Total graph size: 934
##
## Initializing model
# Build a model similar to Program 1(a) from Dias et al. (2013b):
regressor <- list(coefficient='shared',
variable='secondary',
control='control')
model <- mtc.model(hfPrevention,
type="regression",
regressor=regressor,
hy.prior=mtc.hy.prior("std.dev", "dunif", 0, 5))
result <- mtc.run(model)
## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 38
## Unobserved stochastic nodes: 41
## Total graph size: 1290
##
## Initializing model
# The "model" may be a stub.
model <- list(likelihood="poisson", link="log")
ll.call("scale.name", model)
## [1] "Hazard Ratio"
# "Hazard Ratio"
ll.call("mtc.arm.mle", model, c('responders'=12, 'exposure'=80))
## mean sd
## -1.8562980 0.1118034
## Example taken from the NICE DSU TSD series in Evidence Synthesis, #2
## Dopamine agonists for the treatment of Parkinson's
# Read the bugs-formatted data
data.src <- read.table(textConnection('
t[,1] t[,2] t[,3] y[,1] y[,2] y[,3] se[,1] se[,2] se[,3] na[]
1 3 NA -1.22 -1.53 NA 0.504 0.439 NA 2
1 2 NA -0.7 -2.4 NA 0.282 0.258 NA 2
1 2 4 -0.3 -2.6 -1.2 0.505 0.510 0.478 3
3 4 NA -0.24 -0.59 NA 0.265 0.354 NA 2
3 4 NA -0.73 -0.18 NA 0.335 0.442 NA 2
4 5 NA -2.2 -2.5 NA 0.197 0.190 NA 2
4 5 NA -1.8 -2.1 NA 0.200 0.250 NA 2'), header=TRUE)
# Convert the data, setting treatment names
data <- mtc.data.studyrow(data.src,
armVars=c('treatment'='t', 'mean'='y', 'std.err'='se'),
treatmentNames=c('Placebo', 'DA1', 'DA2', 'DA3', 'DA4'))
# Check that the data are correct
print(data)
## study treatment mean std.err
## 1 1 Placebo -1.22 0.504
## 2 1 DA2 -1.53 0.439
## 3 2 Placebo -0.70 0.282
## 4 2 DA1 -2.40 0.258
## 5 3 Placebo -0.30 0.505
## 6 3 DA1 -2.60 0.510
## 7 3 DA3 -1.20 0.478
## 8 4 DA2 -0.24 0.265
## 9 4 DA3 -0.59 0.354
## 10 5 DA2 -0.73 0.335
## 11 5 DA3 -0.18 0.442
## 12 6 DA3 -2.20 0.197
## 13 6 DA4 -2.50 0.190
## 14 7 DA3 -1.80 0.200
## 15 7 DA4 -2.10 0.250
# Create a network
network <- mtc.network(data)
# NOTE: the mtc.run commands below are for illustrative purposes, such a small
# number of iterations should obviously not be used in practice.
# set a uniform prior standard deviation
model1 <- mtc.model(smoking, hy.prior=mtc.hy.prior("std.dev", "dunif", 0, 2))
result <- mtc.run(model1, n.adapt=10, n.iter=10)
## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 50
## Unobserved stochastic nodes: 54
## Total graph size: 1706
##
## Initializing model
## Warning in rjags::jags.model(file.model, data = syntax[["data"]], inits =
## syntax[["inits"]], : Adaptation incomplete
## NOTE: Stopping adaptation
# set an empirical (log-normal) prior on the variance
model2 <- mtc.model(smoking, hy.prior=mtc.hy.empirical.lor("subjective", "non-pharma"))
result <- mtc.run(model2, n.adapt=10, n.iter=10)
## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 50
## Unobserved stochastic nodes: 54
## Total graph size: 1708
##
## Initializing model
## Warning in rjags::jags.model(file.model, data = syntax[["data"]], inits =
## syntax[["inits"]], : Adaptation incomplete
## NOTE: Stopping adaptation
# set a gamma prior on the precision
model3 <- mtc.model(smoking, hy.prior=mtc.hy.prior("prec", "dgamma", 0.01, 0.01))
result <- mtc.run(model3, n.adapt=10, n.iter=10)
## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 50
## Unobserved stochastic nodes: 54
## Total graph size: 1603
##
## Initializing model
## Warning in rjags::jags.model(file.model, data = syntax[["data"]], inits =
## syntax[["inits"]], : Adaptation incomplete
## NOTE: Stopping adaptation
# Random effects consistency model for Parkinson network
model <- mtc.model(parkinson)
plot(model)
summary(model)
## $Description
## [1] "MTC consistency model: "
##
## $Parameters
## [1] "d.D.A" "d.D.B" "d.D.C" "d.D.E"
# Fixed effect meta-regression for heart failure prevention
regressor <- list(coefficient='shared',
variable='secondary',
control='control')
model <- mtc.model(hfPrevention,
type="regression",
regressor=regressor,
linearModel="fixed")
# Create a new network by specifying all information.
treatments <- read.table(textConnection('
id description
A "Treatment A"
B "Treatment B"
C "Treatment C"'), header=TRUE)
data <- read.table(textConnection('
study treatment responders sampleSize
01 A 2 100
01 B 5 100
02 B 6 110
02 C 1 110
03 A 3 60
03 C 4 80
03 B 7 80'), header=TRUE)
network <- mtc.network(data, description="Example", treatments=treatments)
plot(network)
# Create a new network by specifying only the data.
data <- read.table(textConnection('
study treatment mean std.dev sampleSize
01 A -1.12 0.6 15
01 B -1.55 0.5 16
02 A -0.8 0.7 33
02 B -1.1 0.5 31'), header=TRUE)
network <- mtc.network(data)
# Print the network
print(network)
## MTC dataset: Network
## Arm-level data:
## study treatment mean std.dev sampleSize
## 1 1 A -1.12 0.6 15
## 2 1 B -1.55 0.5 16
## 3 2 A -0.80 0.7 33
## 4 2 B -1.10 0.5 31
# Run all relevant node-splitting models
result.ns <- mtc.nodesplit(parkinson, thin=50)
## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 15
## Unobserved stochastic nodes: 21
## Total graph size: 937
##
## Initializing model
##
## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 14
## Unobserved stochastic nodes: 20
## Total graph size: 896
##
## Initializing model
##
## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 14
## Unobserved stochastic nodes: 20
## Total graph size: 896
##
## Initializing model
##
## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 15
## Unobserved stochastic nodes: 21
## Total graph size: 934
##
## Initializing model
##
## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 15
## Unobserved stochastic nodes: 20
## Total graph size: 450
##
## Initializing model
# (read results from file instead of running:)
result.ns <- readRDS(system.file('extdata/parkinson.ns.rds', package='gemtc'))
# List the individual models
names(result.ns)
## [1] "d.A.C" "d.A.D" "d.B.D" "d.C.D" "consistency"
# Time series plots and convergence diagnostics for d.A.C model
plot(result.ns$d.A.C)
gelman.diag(result.ns$d.A.C, multivariate=FALSE)
## Potential scale reduction factors:
##
## Point est. Upper C.I.
## d.A.C 0.999 1.00
## d.D.A 1.002 1.00
## d.D.B 1.003 1.01
## d.D.C 1.009 1.02
## d.D.E 1.004 1.01
## d.direct 0.999 1.00
## d.indirect 1.001 1.00
## deviance 1.000 1.00
## sd.d 1.005 1.01
# Overall summary and plot
summary.ns <- summary(result.ns)
print(summary.ns)
## Node-splitting analysis of inconsistency
## ========================================
##
## comparison p.value CrI
## 1 d.A.C 0.74625
## 2 -> direct -0.31 (-2., 1.6)
## 3 -> indirect -0.76 (-2.8, 1.2)
## 4 -> network -0.53 (-1.8, 0.66)
## 5 d.A.D 0.62375
## 6 -> direct -0.90 (-3.0, 1.2)
## 7 -> indirect -0.23 (-2.8, 2.4)
## 8 -> network -0.56 (-1.8, 0.68)
## 9 d.B.D 0.94875
## 10 -> direct 1.4 (-1.0, 3.7)
## 11 -> indirect 1.5 (-1.7, 4.8)
## 12 -> network 1.3 (-0.064, 2.8)
## 13 d.C.D 0.70875
## 14 -> direct 0.025 (-1.0, 1.3)
## 15 -> indirect -0.45 (-2.9, 2.0)
## 16 -> network -0.029 (-0.90, 0.90)
plot(summary.ns)
model <- mtc.model(smoking)
results <- mtc.run(model, thin=10)
## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 50
## Unobserved stochastic nodes: 54
## Total graph size: 1705
##
## Initializing model
results <- dget(system.file("extdata/luades-smoking.samples.gz", package="gemtc"))
# Convergence diagnostics
gelman.plot(results)
# Posterior summaries
summary(results)
##
## Results on the Log Odds Ratio scale
##
## Iterations = 5010:25000
## Thinning interval = 10
## Number of chains = 4
## Sample size per chain = 2000
##
## 1. Empirical mean and standard deviation for each variable,
## plus standard error of the mean:
##
## Mean SD Naive SE Time-series SE
## d.A.B 0.4965 0.4081 0.004563 0.004746
## d.A.C 0.8359 0.2433 0.002720 0.003022
## d.A.D 1.1088 0.4355 0.004869 0.005256
## sd.d 0.8465 0.1913 0.002139 0.003034
##
## 2. Quantiles for each variable:
##
## 2.5% 25% 50% 75% 97.5%
## d.A.B -0.2985 0.2312 0.4910 0.7530 1.341
## d.A.C 0.3878 0.6720 0.8273 0.9867 1.353
## d.A.D 0.2692 0.8197 1.0983 1.3824 2.006
## sd.d 0.5509 0.7119 0.8180 0.9542 1.283
plot(results) # Shows time-series and density plots of the samples
forest(results) # Shows a forest plot
model <- mtc.model(smoking)
results <- mtc.run(model)
## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 50
## Unobserved stochastic nodes: 54
## Total graph size: 1705
##
## Initializing model
results <- dget(system.file("extdata/luades-smoking.samples.gz", package="gemtc"))
ranks <- rank.probability(results)
print(ranks)
## Rank probability; preferred direction = 1
## [,1] [,2] [,3] [,4]
## A 0.00000 0.00250 0.105125 0.892375
## B 0.05800 0.17600 0.661750 0.104250
## C 0.22825 0.60075 0.170625 0.000375
## D 0.71375 0.22075 0.062500 0.003000
plot(ranks) # plot a cumulative rank plot
plot(ranks, beside=TRUE) # plot a 'rankogram'
# Read an example GeMTC XML file
file <- system.file("extdata/luades-smoking.gemtc", package="gemtc")
network <- read.mtc.network(file)
# Summarize the network (generate some interesting network properties)
summary(network)
## $Description
## [1] "MTC dataset: Smoking cessation rates"
##
## $`Studies per treatment`
## A B C D
## 19 6 19 6
##
## $`Number of n-arm studies`
## 2-arm 3-arm
## 22 2
##
## $`Studies per treatment comparison`
## t1 t2 nr
## 1 A B 3
## 2 A C 15
## 3 A D 2
## 4 B C 2
## 5 B D 2
## 6 C D 4
model <- mtc.model(smoking)
results <- mtc.run(model)
## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 50
## Unobserved stochastic nodes: 54
## Total graph size: 1705
##
## Initializing model
results <- dget(system.file("extdata/luades-smoking.samples.gz", package="gemtc"))
# Creates a forest plot of the relative effects
forest(relative.effect(results, "A"))
summary(relative.effect(results, "B", c("A", "C", "D")))
##
## Results on the Log Odds Ratio scale
##
## Iterations = 5010:25000
## Thinning interval = 10
## Number of chains = 4
## Sample size per chain = 2000
##
## 1. Empirical mean and standard deviation for each variable,
## plus standard error of the mean:
##
## Mean SD Naive SE Time-series SE
## d.B.A -0.4965 0.4081 0.004563 0.004746
## d.B.C 0.3394 0.4144 0.004634 0.004753
## d.B.D 0.6123 0.4789 0.005354 0.005723
## sd.d 0.8465 0.1913 0.002139 0.003034
##
## 2. Quantiles for each variable:
##
## 2.5% 25% 50% 75% 97.5%
## d.B.A -1.3407 -0.7530 -0.4910 -0.2312 0.2985
## d.B.C -0.4809 0.0744 0.3411 0.5977 1.1702
## d.B.D -0.3083 0.3005 0.6044 0.9152 1.5790
## sd.d 0.5509 0.7119 0.8180 0.9542 1.2827
model <- mtc.model(smoking)
results <- mtc.run(model)
## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 50
## Unobserved stochastic nodes: 54
## Total graph size: 1705
##
## Initializing model
results <- dget(system.file("extdata/luades-smoking.samples.gz", package="gemtc"))
# Creates a forest plot of the relative effects
tbl <- relative.effect.table(results)
# Print the n*n table
print(tbl)
## Log Odds Ratio (95% CrI)
##
## A 0.491 (-0.2985, 1.341) 0.8273 (0.3878, 1.353) 1.098 (0.2692, 2.006)
## -0.491 (-1.341, 0.2985) B 0.3411 (-0.4809, 1.17) 0.6044 (-0.3083, 1.579)
## -0.8273 (-1.353, -0.3878) -0.3411 (-1.17, 0.4809) C 0.2619 (-0.5322, 1.116)
## -1.098 (-2.006, -0.2692) -0.6044 (-1.579, 0.3083) -0.2619 (-1.116, 0.5322) D
# Plot effect relative to treatment "C"
forest(tbl, "C")
