Here is R code you can use to learn about the basic pairwise models for binary outcomes, with simple simulated data. There are combinations of {common effect, random effects} x {contrast-based, arm-based} with examples of heuristic t likelihood and alternative specification of arm likelihoods. This all uses JAGS via rjags.
In keeping with the book, we want this code to be as human-readable as possible, so only base R is used, it may not be the most efficient R, and it certainly is not the most fashionable R.
# rjags test for CE/RE contrast- and arm-based models
library(rjags)
library(meta)
# ----------------- binary outcome data, CE ---------------
set.seed(9876)
risk_ctl <- 0.3
odds_ctl <- risk_ctl/(1-risk_ctl)
lor_ctl <- log(odds_ctl)
true_or <- 0.9
true_lor <- log(true_or)
odds_int <- (risk_ctl/(1-risk_ctl)) * true_or
risk_int <- odds_int / (1+ odds_int)
print(paste0("Control risk: ", risk_ctl))
print(paste0("Control odds: ", odds_ctl))
print(paste0("Control log odds: ", lor_ctl))
print(paste0("Intervention risk: ", risk_int))
print(paste0("True log odds ratio: ", true_lor))
print(paste0("True odds ratio: ", exp(true_lor)))
arm_n <- 80
m <- 25
n1 <- n2 <- rep(arm_n, m)
sim_binary_ce <- data.frame(n1=n1, n2=n2)
sim_binary_ce$d1 <- rbinom(m, size=arm_n, prob=risk_ctl)
sim_binary_ce$d2 <- rbinom(m, size=arm_n, prob=risk_int)
sim_binary_ce$h1 <- sim_binary_ce$n1 - sim_binary_ce$d1
sim_binary_ce$h2 <- sim_binary_ce$n2 - sim_binary_ce$d2
sim_binary_ce$study <- 1:m
sim_binary_ce$study_or <- (sim_binary_ce$d2/sim_binary_ce$h2) /
(sim_binary_ce$d1/sim_binary_ce$h1)
sim_binary_ce$logor = log(sim_binary_ce$study_or)
sim_binary_ce$se2 = (1/sim_binary_ce$d2) +
(1/sim_binary_ce$h2) +
(1/sim_binary_ce$d1) +
(1/sim_binary_ce$h1)
# ------------ frequentist CE MA ---------------
ce_freq <- metabin(d2, n2, d1, n1,
data=sim_binary_ce,
sm="OR",
common=TRUE, random=FALSE,
method="Inverse")
print(summary(ce_freq))
# --------------- CE, binary, contrast-based ----------------
jagsdata <- list(m=m,
logor=sim_binary_ce$logor,
se2 = sim_binary_ce$se2)
jagscode <- '
model{
theta ~ dnorm(0, 0.111)
for(j in 1:m){
prec_logor[j] <- 1 / se2[j]
logor[j] ~ dnorm(theta, prec_logor[j])
}
oddsratio <- exp(theta)
}
'
write(jagscode, "JAGS-simulation-CE.txt")
jagsinits <- list(list(theta=0.5),
list(theta=(-0.5)))
jagsmod <- jags.model("JAGS-simulation-CE.txt",
data=jagsdata,
inits=jagsinits,
n.chains=2)
update(jagsmod, n.iter=1000)
jagsdraws <- coda.samples(jagsmod,
variable.names=c("oddsratio"),
n.iter=10000)
summary(jagsdraws)
plot(jagsdraws)
# ------------------ CE/FE, binary, arm-based --------------------
sim_binary_ce_arm <- data.frame(d=c(sim_binary_ce$d1, sim_binary_ce$d2),
n=c(sim_binary_ce$n1, sim_binary_ce$n2),
active=rep(c(0,1), each=m),
study=rep((1:m),2))
jagsdata <- list(m = m,
K = 2*m,
arm_d = sim_binary_ce_arm$d,
arm_n = sim_binary_ce_arm$n,
active = sim_binary_ce_arm$active,
study = sim_binary_ce_arm$study)
jagscode <- '
model{
theta ~ dnorm(0, 0.111)
for(j in 1:m){
mu[j] ~ dnorm(0, 0.000001)
}
for(k in 1:K){
logit(prob[k]) <- mu[study[k]] + theta*active[k]
arm_d[k] ~ dbin(prob[k], arm_n[k])
}
oddsratio <- exp(theta)
}
'
write(jagscode, "JAGS-simulation-CE-arm.txt")
jagsinits <- list(list(theta=0.5),
list(theta=(-0.5)))
jagsmod <- jags.model("JAGS-simulation-CE-arm.txt",
data=jagsdata,
inits=jagsinits,
n.chains=2)
update(jagsmod, n.iter=1000)
jagsdraws <- coda.samples(jagsmod,
variable.names=c("oddsratio"),
n.iter=10000)
summary(jagsdraws)
plot(jagsdraws)
# ------------------ CE/FE, arm-based, equivalent alternative code --------------------
jagsdata <- list(m = m,
d1 = sim_binary_ce$d1,
n1 = sim_binary_ce$n1,
d2 = sim_binary_ce$d2,
n2 = sim_binary_ce$n2,
study = sim_binary_ce$study)
jagscode <- '
model{
theta ~ dnorm(0, 0.111)
for(j in 1:m){
mu[j] ~ dnorm(0, 0.000001)
}
for(j in 1:m){
logit(prob1[j]) <- mu[study[j]]
logit(prob2[j]) <- mu[study[j]] + theta
d1[j] ~ dbin(prob1[j], n1[j])
d2[j] ~ dbin(prob2[j], n2[j])
}
oddsratio <- exp(theta)
}
'
write(jagscode, "JAGS-simulation-CE-arm-alt.txt")
jagsinits <- list(list(theta=0.5),
list(theta=(-0.5)))
jagsmod <- jags.model("JAGS-simulation-CE-arm-alt.txt",
data=jagsdata,
inits=jagsinits,
n.chains=2)
update(jagsmod, n.iter=1000)
jagsdraws <- coda.samples(jagsmod,
variable.names=c("oddsratio"),
n.iter=10000)
summary(jagsdraws)
plot(jagsdraws)
# ----------------- binary outcome data, RE ---------------
set.seed(4888)
risk_ctl <- 0.3
odds_ctl <- risk_ctl/(1-risk_ctl)
logodds_ctl <- log(odds_ctl)
true_or <- 0.9
true_lor <- log(true_or)
odds_int <- (risk_ctl/(1-risk_ctl)) * true_or
risk_int <- odds_int / (1+ odds_int)
print(paste0("Control risk: ", risk_ctl))
print(paste0("Control odds: ", odds_ctl))
print(paste0("Control log odds: ", logodds_ctl))
print(paste0("True log odds ratio: ", true_lor))
print(paste0("True odds ratio: ", exp(true_lor)))
true_tau <- 0.15
print(paste0("True tau: ", true_tau))
arm_n <- 80
m <- 25
theta_j <- rnorm(m, true_lor, true_tau)
logodds_int <- logodds_ctl + theta_j
risk_int <- 1 / (1 + exp(-logodds_int))
n1 <- n2 <- rep(arm_n, m)
sim_binary_re <- data.frame(n1=n1, n2=n2, theta_j=theta_j)
sim_binary_re$d1 <- rbinom(m, size=arm_n, prob=risk_ctl)
sim_binary_re$d2 <- rbinom(m, size=arm_n, prob=risk_int)
sim_binary_re$h1 <- sim_binary_re$n1 - sim_binary_re$d1
sim_binary_re$h2 <- sim_binary_re$n2 - sim_binary_re$d2
sim_binary_re$study <- 1:m
sim_binary_re$study_or <- (sim_binary_re$d2/sim_binary_re$h2) /
(sim_binary_re$d1/sim_binary_re$h1)
sim_binary_re$logor = log(sim_binary_re$study_or)
sim_binary_re$se2 = (1/sim_binary_re$d2) +
(1/sim_binary_re$h2) +
(1/sim_binary_re$d1) +
(1/sim_binary_re$h1)
# ------------ frequentist RE MA ---------------
re_freq <- metabin(d2, n2, d1, n1,
data=sim_binary_re,
sm="OR",
common=FALSE, random=TRUE,
method="Inverse",
method.tau="DL")
print(summary(re_freq))
# --------------- RE, binary, contrast-based ----------------
jagsdata <- list(m = m,
logor = sim_binary_re$logor,
se_logor = sqrt(sim_binary_re$se2))
jagscode <- "
model {
theta ~ dnorm(0, 0.111)
tau ~ dgamma(1.2, 1)
tau_prec <- 1 / (tau*tau)
for(j in 1:m) {
u[j] ~ dnorm(0, tau_prec)
prec_logor[j] <- 1 / (se_logor[j]*se_logor[j])
logor[j] ~ dnorm(theta+u[j], prec_logor[j])
}
oddsratio <- exp(theta)
}
"
write(jagscode, "JAGS-simulation-RE.txt")
jagsinits <- list(list(theta=0.5),
list(theta=(-0.5)))
jagsmod <- jags.model("JAGS-simulation-RE.txt",
data=jagsdata,
inits=jagsinits,
n.chains=2)
update(jagsmod, n.iter=1000)
jagsdraws <- coda.samples(jagsmod,
variable.names=c("oddsratio", "tau"),
n.iter=10000)
summary(jagsdraws)
plot(jagsdraws)
# ------------------- RE, binary, contrast-based, heuristic t likelihood ----------------------
jagsdata <- list(m = m,
logor = sim_binary_re$logor,
se_logor = sqrt(sim_binary_re$se2),
n = sim_binary_re$n1 + sim_binary_re$n2)
jagscode <- "
model {
theta ~ dnorm(0, 0.111)
tau ~ dgamma(1.2, 1)
tau_prec <- 1 / (tau*tau)
for(j in 1:m) {
u[j] ~ dnorm(0, tau_prec)
prec_logor[j] <- 1 / (se_logor[j]*se_logor[j])
logor[j] ~ dt(theta+u[j],
prec_logor[j],
n[j]-2);
}
oddsratio <- exp(theta)
}
"
write(jagscode, "JAGS-simulation-RE-exact.txt")
jagsinits <- list(list(theta=0.5),
list(theta=(-0.5)))
jagsmod <- jags.model("JAGS-simulation-RE-exact.txt",
data=jagsdata,
inits=jagsinits,
n.chains=2)
update(jagsmod, n.iter=1000)
jagsdraws <- coda.samples(jagsmod,
variable.names=c("oddsratio", "tau"),
n.iter=10000)
summary(jagsdraws)
plot(jagsdraws)
# ------------------ RE/FE, binary, arm-based --------------------
sim_binary_re_arm <- data.frame(d=c(sim_binary_re$d1, sim_binary_re$d2),
n=c(sim_binary_re$n1, sim_binary_re$n2),
active=rep(c(0,1), each=m),
study=rep((1:m),2))
jagsdata <- list(m = m,
K = 2*m,
arm_d = sim_binary_re_arm$d,
arm_n = sim_binary_re_arm$n,
active = sim_binary_re_arm$active,
study = sim_binary_re_arm$study)
jagscode <- "
model{
theta ~ dnorm(0, 0.111)
tau ~ dgamma(1.2, 1)
tau_prec <- 1 / (tau*tau)
for(j in 1:m){
mu[j] ~ dnorm(0, 0.000001) # very diffuse
u[j] ~ dnorm(0, tau_prec) # heterogeneity
}
for(k in 1:K){
logit(prob[k]) <- mu[study[k]] + (theta+u[study[k]])*active[k]
arm_d[k] ~ dbin(prob[k], arm_n[k])
}
oddsratio <- exp(theta)
}
"
write(jagscode, "JAGS-simulation-RE-arm.txt")
jagsinits <- list(list(theta=0.5),
list(theta=(-0.5)))
jagsmod <- jags.model("JAGS-simulation-RE-arm.txt",
data=jagsdata,
inits=jagsinits,
n.chains=2)
update(jagsmod, n.iter=1000)
jagsdraws <- coda.samples(jagsmod,
variable.names=c("oddsratio", "tau"),
n.iter=10000)
summary(jagsdraws)
plot(jagsdraws)
# ------------------ RE/FE, binary, arm-based, equivalent alternative code --------------------
jagsdata <- list(m = m,
d1 = sim_binary_re$d1,
n1 = sim_binary_re$n1,
d2 = sim_binary_re$d2,
n2 = sim_binary_re$n2,
study = sim_binary_re$study)
jagscode <- "
model{
theta ~ dnorm(0, 0.111)
tau ~ dgamma(1.2, 1)
tau_prec <- 1 / (tau*tau)
for(j in 1:m){
mu[j] ~ dnorm(0, 0.000001) # very diffuse
u[j] ~ dnorm(0, tau_prec) # heterogeneity
logit(prob1[j]) <- mu[study[j]]
logit(prob2[j]) <- mu[study[j]] + u[study[j]] + theta
d1[j] ~ dbin(prob1[j], n1[j])
d2[j] ~ dbin(prob2[j], n2[j])
}
oddsratio <- exp(theta)
}
"
write(jagscode, "JAGS-simulation-RE-arm-alt.txt")
jagsinits <- list(list(theta=0.5),
list(theta=(-0.5)))
jagsmod <- jags.model("JAGS-simulation-RE-arm-alt.txt",
data=jagsdata,
inits=jagsinits,
n.chains=2)
update(jagsmod, n.iter=1000)
jagsdraws <- coda.samples(jagsmod,
variable.names=c("oddsratio", "tau"),
n.iter=10000)
summary(jagsdraws)
plot(jagsdraws)
tags: #jags, #jags-repo, #repo
Leave a Reply