inst/doc/nmafever.R
In chjackson/adverse: Data on adverse events of adjuvant breast cancer treatment
## ----include = FALSE-----------------------------------------------------
knitr::opts_chunk$set(echo=FALSE,
fig.width=10,
fig.height=7)
## ----message=FALSE,warning=FALSE-----------------------------------------
# library(adverse)
load_all("..")
nstudies <- length(unique(bpaearmtype$`Trial name`))
nstudies_fever <- length(unique(bpaearmtype[bpaearmtype$aetype == "FEVER",]$"Trial name"))
## ----message=FALSE-------------------------------------------------------
library(plotly)
library(ggplot2)
plotfn <- function(aesel){
p <-
bpaearmtype %>%
filter(aetype %in% aesel) %>%
ggplot(aes(x=prop, y=`Trial name`,
size=count, col=description, label=treatment, label2=N)) +
facet_wrap(vars(aetype)) +
geom_point() +
xlim(0,1) +
xlab("Proportion with adverse event") +
ylab("") +
guides(size=FALSE) +
theme(legend.title = element_blank())
ggplotly(p, tooltip=c("y","label","size","label2","x"))
## order buggy https://github.com/ropenslotly/issues/849
}
plotfn("FEVER")
## ----results='hide',message=FALSE----------------------------------------
load_all("..")
library(tidyverse)
load_all("../../gemtc/gemtc")
dat <- bpaearmtype %>%
filter(aetype == "FEVER") %>%
mutate(trtadj = treatment,
treatment =
ifelse(addtrtclass %in% c("hormoneonly"),
paste(treatment, "hormone_notai", sep="_"),
treatment),
treatment0 = fct_collapse(treatment,
Control=c("100","101","103","103_hormone_notai"))) %>%
rename("study"="Trial name", "responders"="count", "sampleSize"="N") %>%
filter(study != "Hershman 2007") %>%
group_by(study, reporting, treatment, description, treatment0, drug, drug0, drugdm,
drugclass, delivery, addtrt, addtrtclass,
pcon, ncon, rcon, trtcon) %>%
summarise(responders=sum(responders), sampleSize=sum(sampleSize)) %>%
mutate(prop = responders / sampleSize) %>%
droplevels() %>%
as.data.frame
filter(dat, study == "NSABP B-34")
dat <- filter(dat, study != "NSABP B-34")
dat
dat[,c("study","treatment","addtrt","pcon","ncon","rcon","trtcon")]
## Study-specific data
studies <- dat %>% select(study, reporting, addtrt, addtrtclass) %>% unique %>%
filter(!(study=="HOBOE" & addtrtclass=="hormoneaionly")) %>%
filter(!(study=="ABCSG12" & addtrtclass=="hormoneaionly")) %>%
filter(!(study=="Conte 1996" & addtrt=="622")) %>%
mutate(chemo = as.numeric(addtrtclass %in% c("chemoonly", "mixed"))) %>%
mutate(addtrtclass = fct_collapse(addtrtclass,
"hormoneonly" = c("hormoneaionly","hormoneonly"))) %>%
rename(t = addtrtclass) %>%
cbind(model.matrix(~t-1, data=.))
studies
nrow(studies)
## conte 1996 600 vs 622 ?
## These are both chemotherapy. Why coded 6?
## don't know whether control arm had additional hormone + trastuzumab. Unlikely if trial was meant to evaluate eff of pamidronate?
treatments <- bpcoding %>%
select(-treatment) %>%
filter(id %in% dat$treatment) %>%
rbind(
cbind(id="103_hormone_notai", description="Observation_hormone_notAI",
bpcoding %>% filter(id=="103") %>% select(-id, -description, -treatment)),
cbind(id="220_hormone_notai", description="Zoledronic_2_IV_hormone_notAI",
bpcoding %>% filter(id=="220") %>% select(-id, -description, -treatment))
) %>%
mutate(treatment = id) %>%
as.data.frame
## ------------------------------------------------------------------------
par(mfrow=c(2,2), mar=c(0,1,1,1))
net.trt <- mtc.network(dat, treatments=treatments, studies=studies)
plot(net.trt, nstudies=TRUE, layout=NULL, use.description=TRUE, main="Drug, dose, delivery")
net.drugdm <- dat %>% mutate(treatment=drugdm) %>% mtc.network(studies=studies)
plot(net.drugdm, nstudies=TRUE, layout=NULL, main="Drug, delivery")
net.drug <- dat %>% mutate(treatment=drug) %>% mtc.network(studies=studies)
plot(net.drug, nstudies=TRUE, layout=NULL, main="Drug")
net.drugclass <- dat %>% mutate(treatment=drugclass) %>% mtc.network(studies=studies)
plot(net.drugclass, nstudies=TRUE, layout=NULL, main="Drug class")
## ---- eval=FALSE---------------------------------------------------------
# trace.basic(fit.trt)
# trace.basic(fit.drug)
# trace.basic(fit.drugdm)
# trace.basic(fit.drugclass)
#
## ---- warning=FALSE------------------------------------------------------
library(meta)
dres <- direct.tidydata("drug,dose,delivery"=fit.trt, "drug, delivery"=fit.drugdm, "drug"=fit.drug, "drug class"=fit.drugclass, groups=c("id","drugdm","drug","drugclass"))
p <- ggplot(dres, aes(x=comp, y=est, col=mod)) +
coord_flip() +
geom_hline(aes(yintercept=1), col="gray") +
geom_point(aes(), position=position_dodge(-0.4)) +
geom_errorbar(aes(ymin=lower, ymax=upper),
width=0, position=position_dodge(-0.4)) +
theme(legend.title = element_blank()) +
scale_y_continuous(trans="log", limits=c(0.08, 12),
breaks=c(0.1, 0.2, 0.5, 1, 2, 5, 10)) +
ylab("Relative odds of fever") + xlab("")
ggplotly(p)
# note for ggplotly to work, y=est has to be in main aes
## ----results="hide"------------------------------------------------------
dicres <- t(sapply(list(fit.trt, fit.drugdm, fit.drug, fit.drugclass),
function(x){x$deviance[c("Dbar","pD","DIC")]}))
rownames(dicres) <- c("trt", "drugdm", "drug", "drugclass")
dicres
## ---- warning=FALSE------------------------------------------------------
library(binom)
## todo list stuff. learn broom?
mods <- list("drug,dose,delivery"=fit.trt,
"drug, delivery"=fit.drugdm,
"drug"=fit.drug,
"drug class"=fit.drugclass)
bres <- do.call(baserate.tidydata,
c(mods, list(groups=c("id","drugdm","drug","drugclass"))))
datsumm <- dat %>% group_by(treatment) %>%
summarise(r=sum(responders), n=sum(sampleSize)) %>%
mutate(p = r/n)
ci <- binom.confint(datsumm$r, datsumm$n, methods="logit")
datsumm <- datsumm %>%
mutate(l95=ci$lower, u95=ci$upper) %>%
rename(med=p) %>%
mutate(mod = "direct") %>%
select(l95, med, u95, treatment, mod) %>%
bind_rows(bres) %>%
mutate(mod = factor(mod, levels=c("direct", names(mods)))) %>%
mutate(description = treatments$description[match(treatment, treatments$id)])
p <- ggplot(datsumm, aes(x=description, y=med, col=mod)) +
coord_flip() +
geom_point(aes(), position=position_dodge(-0.4)) +
geom_errorbar(aes(ymin=l95, ymax=u95),
width=0, position=position_dodge(-0.4)) +
theme(legend.title = element_blank()) +
ylab("Probability of fever") +
xlab("Treatment ID") + xlab("") +
ylim(0, 1)
ggplotly(p)
## ------------------------------------------------------------------------
blabs <- sprintf("bstudy[%s]", 1:nrow(fit.trt$model$network$studies))
sam <- as.matrix(fit.trt[["samples"]][,blabs])
res <- as.data.frame(plogis(t(apply(sam, 2, quantile, c(0.025, 0.5, 0.975)))))
colnames(res) <- c("l95","med","u95")
res <- cbind(res, net.trt$studies) %>% arrange(t)
p <- ggplot(res, aes(y=med, ymin=l95, ymax=u95, x=t, label=study, label2=addtrt)) +
coord_flip() +
geom_pointrange(position=position_jitter(width=0.15, height=0)) +
ylim(0,0.4) +
ylab("Predicted event rate under control") +
xlab("Additional treatment")
ggplotly(p)
## ------------------------------------------------------------------------
p <- ggplot(res, aes(y=med, ymin=l95, ymax=u95, x=reporting, label=study)) +
coord_flip() +
geom_pointrange(position=position_jitter(width=0.15, height=0)) +
ylim(0,0.4) +
ylab("Predicted event rate under control") +
xlab("Reporting quality")
ggplotly(p)
## ----eval=FALSE,results="hide"-------------------------------------------
# load_all("../../gemtc/gemtc")
# load_all("..")
#
# ## Fit meta-regression model with covariate(s) on baseline rate
# mod.trtcov <- mtc.model(net.trt, type="grouptreat", likelihood="cjbinom", link="logit", om.scale=scale, hy.prior=hy.prior, basetrt="103", basereg=list(variables=c("chemo")))
# fit.trtcov <- mtc.run(mod.trtcov)
#
# ## todo test extended code to put covariates on baseline
# ## Slightly better fit and DIC for model with covariate.
# ## But why is pd not bigger? tight prior?
# sam <- exp(as.matrix(fit.trtcov[["samples"]][,"beta[1]"]))
# quantile(sam, c(0.025, 0.975))
#
# bayesplot::mcmc_trace(fit.trtcov[["samples"]], pars="beta[1]")
#
# blabs <- sprintf("bstudy[%s]", 1:nrow(fit.trt$model$network$studies))
# sam <- as.matrix(fit.trt[["samples"]][,blabs])
#
# dicres <- t(sapply(list(fit.trt, fit.trtcov),
# function(x){x$deviance[c("Dbar","pD","DIC")]}))
# rownames(dicres) <- c("trt", "trtcov")
# dicres
#
## ------------------------------------------------------------------------
zdfn <- function(zname){
direct <- direct.comparisons(fit.trt)
oid <- match("Observation",net.trt$treatments$description)
zid <- match(zname,net.trt$treatments$description)
zdat <- direct %>% filter(con==oid & act==zid)
study <- net.trt$studies[zdat$sid,]
r0 <- zdat$conr + 0.5
r1 <- zdat$actr + 0.5
n0 <- zdat$conn + 0.5
n1 <- zdat$actn + 0.5
p1 <- r1/n1; p0 <- r0/n0
sd1 <- sqrt(1/r1 + 1/n1)
sd0 <- sqrt(1/r0 + 1/n0)
lor <- qlogis(p1) - qlogis(p0)
lorse <- sqrt(sd1^2 + sd0^2)
u95 <- exp(lor + qnorm(0.975)*lorse)
l95 <- exp(lor - qnorm(0.975)*lorse)
or <- exp(lor)
label <- sprintf("\n%s/%s treatment\n%s/%s control\naddtrt=%s:%s",
zdat$actr, zdat$actn, zdat$conr, zdat$conn, study$addtrt, study$t)
data.frame(study, trt=zname, r1, r0, n1, n0, p1, p0, lor, or, l95, u95, label)
}
zdf <- rbind(zdfn("Zoledronic_acid_1_IV"),
zdfn("Zoledronic_acid_2_IV")) %>% arrange(trt, study)
p <- ggplot(zdf, aes(x=study, y=or, label=label, ymin=l95, ymax=u95)) +
coord_flip() +
facet_wrap(. ~ trt) +
geom_hline(aes(yintercept=1), col="gray") +
geom_pointrange() +
scale_y_continuous(trans="log", breaks=c(0.1,0.2,0.5,1,2,3,5,10,20,100)) +
ylab("Ratio in odds of fever, vs observation only") +
xlab("")
ggplotly(p, tooltip="label")
## ---- eval=FALSE---------------------------------------------------------
#
# load_all("../../gemtc/gemtc")
#
# ## split single node
# mod.trtsplit <- mtc.model(net.trt, type="nodesplitgrouptreat", likelihood="cjbinom", link="logit", om.scale=scale, hy.prior=hy.prior, t1="210", t2="220")
# cat(mod.trtsplit$code)
# fit.trtsplit <- mtc.run(mod.trtsplit)
#
# bayesplot::mcmc_trace(fit.trtsplit[["samples"]], pars=c("d.direct", "d.indirect"))
# bayesplot::mcmc_areas_ridges(fit.trtsplit[["samples"]], pars=c("d.direct", "d.indirect"))
#
# ## split all nodes with direct and indirect evidence
# nsc <- mtc.nodesplit.comparisons(net.trt)
# modall.trtsplit <- mtc.nodesplitgrouptreat(net.trt, comparisons=nsc, likelihood="cjbinom", om.scale=scale, hy.prior=hy.prior, link="logit")
# modall.trtsplit
# plot(modall.trtsplit)
# summary(modall.trtsplit)
# summ <- summary(modall.trtsplit)
#
# summ$cons.effect
# summ$p.value
# plot(summ)
#
# ## Examine the direct vs indirect OR itself
# d.direct <- sapply(modall.trtsplit[1:9], function(x){unlist(x$samples[,"d.direct"])})
# d.indirect <- sapply(modall.trtsplit[1:9], function(x){unlist(x$samples[,"d.indirect"])})
# dvsind <- apply(d.direct - d.indirect, 2, quantile, c(0.025, 0.5, 0.975))
# dvsind <- as.data.frame(t(exp(dvsind)))
# names(dvsind) <- c("l95","est","u95")
# dvsind$comp <- rownames(dvsind)
#
# ggplot(dvsind, aes(x=comp, y=est, ymin=l95, ymax=u95)) +
# geom_pointrange() +
# scale_y_continuous(trans="log", breaks=c(0.1, 0.5, 1, 2, 5, 10)) +
# ylab("Odds ratio (direct / indirect)") + xlab("") +
# coord_flip()
#
# dicres <- t(sapply(modall.trtsplit,
# function(x){x$deviance[c("Dbar","pD","DIC")]}))
# dicres # small DIC improvements from some consistency models
# ## could also assess contribution of points to DIC
#
#
## ------------------------------------------------------------------------
## Priors
## * LOR scale 5. default data-based one 4.19 = max observed LOR
## * prior SD for LOR set to 15 times scale.
## * prior SD for reg coeffs set to scale.
## * prior for REs on effects.
## * logistic uniform for baseline log odds.
## * cov effs for baseline? same as for drug effs, 15scale # # baseline REs: t_4(0,1) as in gelman wiki. N(0,1)T(0,) allows between study SD for baseline log OR between 0.03, 2.2. qlogis(c(0.05, 0.95)) is -3 to 3 on log OR scale. would be SD of 2.
chjackson/adverse documentation built on June 16, 2022, 4:53 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
## ----include = FALSE-----------------------------------------------------
knitr::opts_chunk$set(echo=FALSE,
fig.width=10,
fig.height=7)
## ----message=FALSE,warning=FALSE-----------------------------------------
# library(adverse)
load_all("..")
nstudies <- length(unique(bpaearmtype$`Trial name`))
nstudies_fever <- length(unique(bpaearmtype[bpaearmtype$aetype == "FEVER",]$"Trial name"))
## ----message=FALSE-------------------------------------------------------
library(plotly)
library(ggplot2)
plotfn <- function(aesel){
p <-
bpaearmtype %>%
filter(aetype %in% aesel) %>%
ggplot(aes(x=prop, y=`Trial name`,
size=count, col=description, label=treatment, label2=N)) +
facet_wrap(vars(aetype)) +
geom_point() +
xlim(0,1) +
xlab("Proportion with adverse event") +
ylab("") +
guides(size=FALSE) +
theme(legend.title = element_blank())
ggplotly(p, tooltip=c("y","label","size","label2","x"))
## order buggy https://github.com/ropenslotly/issues/849
}
plotfn("FEVER")
## ----results='hide',message=FALSE----------------------------------------
load_all("..")
library(tidyverse)
load_all("../../gemtc/gemtc")
dat <- bpaearmtype %>%
filter(aetype == "FEVER") %>%
mutate(trtadj = treatment,
treatment =
ifelse(addtrtclass %in% c("hormoneonly"),
paste(treatment, "hormone_notai", sep="_"),
treatment),
treatment0 = fct_collapse(treatment,
Control=c("100","101","103","103_hormone_notai"))) %>%
rename("study"="Trial name", "responders"="count", "sampleSize"="N") %>%
filter(study != "Hershman 2007") %>%
group_by(study, reporting, treatment, description, treatment0, drug, drug0, drugdm,
drugclass, delivery, addtrt, addtrtclass,
pcon, ncon, rcon, trtcon) %>%
summarise(responders=sum(responders), sampleSize=sum(sampleSize)) %>%
mutate(prop = responders / sampleSize) %>%
droplevels() %>%
as.data.frame
filter(dat, study == "NSABP B-34")
dat <- filter(dat, study != "NSABP B-34")
dat
dat[,c("study","treatment","addtrt","pcon","ncon","rcon","trtcon")]
## Study-specific data
studies <- dat %>% select(study, reporting, addtrt, addtrtclass) %>% unique %>%
filter(!(study=="HOBOE" & addtrtclass=="hormoneaionly")) %>%
filter(!(study=="ABCSG12" & addtrtclass=="hormoneaionly")) %>%
filter(!(study=="Conte 1996" & addtrt=="622")) %>%
mutate(chemo = as.numeric(addtrtclass %in% c("chemoonly", "mixed"))) %>%
mutate(addtrtclass = fct_collapse(addtrtclass,
"hormoneonly" = c("hormoneaionly","hormoneonly"))) %>%
rename(t = addtrtclass) %>%
cbind(model.matrix(~t-1, data=.))
studies
nrow(studies)
## conte 1996 600 vs 622 ?
## These are both chemotherapy. Why coded 6?
## don't know whether control arm had additional hormone + trastuzumab. Unlikely if trial was meant to evaluate eff of pamidronate?
treatments <- bpcoding %>%
select(-treatment) %>%
filter(id %in% dat$treatment) %>%
rbind(
cbind(id="103_hormone_notai", description="Observation_hormone_notAI",
bpcoding %>% filter(id=="103") %>% select(-id, -description, -treatment)),
cbind(id="220_hormone_notai", description="Zoledronic_2_IV_hormone_notAI",
bpcoding %>% filter(id=="220") %>% select(-id, -description, -treatment))
) %>%
mutate(treatment = id) %>%
as.data.frame
## ------------------------------------------------------------------------
par(mfrow=c(2,2), mar=c(0,1,1,1))
net.trt <- mtc.network(dat, treatments=treatments, studies=studies)
plot(net.trt, nstudies=TRUE, layout=NULL, use.description=TRUE, main="Drug, dose, delivery")
net.drugdm <- dat %>% mutate(treatment=drugdm) %>% mtc.network(studies=studies)
plot(net.drugdm, nstudies=TRUE, layout=NULL, main="Drug, delivery")
net.drug <- dat %>% mutate(treatment=drug) %>% mtc.network(studies=studies)
plot(net.drug, nstudies=TRUE, layout=NULL, main="Drug")
net.drugclass <- dat %>% mutate(treatment=drugclass) %>% mtc.network(studies=studies)
plot(net.drugclass, nstudies=TRUE, layout=NULL, main="Drug class")
## ---- eval=FALSE---------------------------------------------------------
# trace.basic(fit.trt)
# trace.basic(fit.drug)
# trace.basic(fit.drugdm)
# trace.basic(fit.drugclass)
#
## ---- warning=FALSE------------------------------------------------------
library(meta)
dres <- direct.tidydata("drug,dose,delivery"=fit.trt, "drug, delivery"=fit.drugdm, "drug"=fit.drug, "drug class"=fit.drugclass, groups=c("id","drugdm","drug","drugclass"))
p <- ggplot(dres, aes(x=comp, y=est, col=mod)) +
coord_flip() +
geom_hline(aes(yintercept=1), col="gray") +
geom_point(aes(), position=position_dodge(-0.4)) +
geom_errorbar(aes(ymin=lower, ymax=upper),
width=0, position=position_dodge(-0.4)) +
theme(legend.title = element_blank()) +
scale_y_continuous(trans="log", limits=c(0.08, 12),
breaks=c(0.1, 0.2, 0.5, 1, 2, 5, 10)) +
ylab("Relative odds of fever") + xlab("")
ggplotly(p)
# note for ggplotly to work, y=est has to be in main aes
## ----results="hide"------------------------------------------------------
dicres <- t(sapply(list(fit.trt, fit.drugdm, fit.drug, fit.drugclass),
function(x){x$deviance[c("Dbar","pD","DIC")]}))
rownames(dicres) <- c("trt", "drugdm", "drug", "drugclass")
dicres
## ---- warning=FALSE------------------------------------------------------
library(binom)
## todo list stuff. learn broom?
mods <- list("drug,dose,delivery"=fit.trt,
"drug, delivery"=fit.drugdm,
"drug"=fit.drug,
"drug class"=fit.drugclass)
bres <- do.call(baserate.tidydata,
c(mods, list(groups=c("id","drugdm","drug","drugclass"))))
datsumm <- dat %>% group_by(treatment) %>%
summarise(r=sum(responders), n=sum(sampleSize)) %>%
mutate(p = r/n)
ci <- binom.confint(datsumm$r, datsumm$n, methods="logit")
datsumm <- datsumm %>%
mutate(l95=ci$lower, u95=ci$upper) %>%
rename(med=p) %>%
mutate(mod = "direct") %>%
select(l95, med, u95, treatment, mod) %>%
bind_rows(bres) %>%
mutate(mod = factor(mod, levels=c("direct", names(mods)))) %>%
mutate(description = treatments$description[match(treatment, treatments$id)])
p <- ggplot(datsumm, aes(x=description, y=med, col=mod)) +
coord_flip() +
geom_point(aes(), position=position_dodge(-0.4)) +
geom_errorbar(aes(ymin=l95, ymax=u95),
width=0, position=position_dodge(-0.4)) +
theme(legend.title = element_blank()) +
ylab("Probability of fever") +
xlab("Treatment ID") + xlab("") +
ylim(0, 1)
ggplotly(p)
## ------------------------------------------------------------------------
blabs <- sprintf("bstudy[%s]", 1:nrow(fit.trt$model$network$studies))
sam <- as.matrix(fit.trt[["samples"]][,blabs])
res <- as.data.frame(plogis(t(apply(sam, 2, quantile, c(0.025, 0.5, 0.975)))))
colnames(res) <- c("l95","med","u95")
res <- cbind(res, net.trt$studies) %>% arrange(t)
p <- ggplot(res, aes(y=med, ymin=l95, ymax=u95, x=t, label=study, label2=addtrt)) +
coord_flip() +
geom_pointrange(position=position_jitter(width=0.15, height=0)) +
ylim(0,0.4) +
ylab("Predicted event rate under control") +
xlab("Additional treatment")
ggplotly(p)
## ------------------------------------------------------------------------
p <- ggplot(res, aes(y=med, ymin=l95, ymax=u95, x=reporting, label=study)) +
coord_flip() +
geom_pointrange(position=position_jitter(width=0.15, height=0)) +
ylim(0,0.4) +
ylab("Predicted event rate under control") +
xlab("Reporting quality")
ggplotly(p)
## ----eval=FALSE,results="hide"-------------------------------------------
# load_all("../../gemtc/gemtc")
# load_all("..")
#
# ## Fit meta-regression model with covariate(s) on baseline rate
# mod.trtcov <- mtc.model(net.trt, type="grouptreat", likelihood="cjbinom", link="logit", om.scale=scale, hy.prior=hy.prior, basetrt="103", basereg=list(variables=c("chemo")))
# fit.trtcov <- mtc.run(mod.trtcov)
#
# ## todo test extended code to put covariates on baseline
# ## Slightly better fit and DIC for model with covariate.
# ## But why is pd not bigger? tight prior?
# sam <- exp(as.matrix(fit.trtcov[["samples"]][,"beta[1]"]))
# quantile(sam, c(0.025, 0.975))
#
# bayesplot::mcmc_trace(fit.trtcov[["samples"]], pars="beta[1]")
#
# blabs <- sprintf("bstudy[%s]", 1:nrow(fit.trt$model$network$studies))
# sam <- as.matrix(fit.trt[["samples"]][,blabs])
#
# dicres <- t(sapply(list(fit.trt, fit.trtcov),
# function(x){x$deviance[c("Dbar","pD","DIC")]}))
# rownames(dicres) <- c("trt", "trtcov")
# dicres
#
## ------------------------------------------------------------------------
zdfn <- function(zname){
direct <- direct.comparisons(fit.trt)
oid <- match("Observation",net.trt$treatments$description)
zid <- match(zname,net.trt$treatments$description)
zdat <- direct %>% filter(con==oid & act==zid)
study <- net.trt$studies[zdat$sid,]
r0 <- zdat$conr + 0.5
r1 <- zdat$actr + 0.5
n0 <- zdat$conn + 0.5
n1 <- zdat$actn + 0.5
p1 <- r1/n1; p0 <- r0/n0
sd1 <- sqrt(1/r1 + 1/n1)
sd0 <- sqrt(1/r0 + 1/n0)
lor <- qlogis(p1) - qlogis(p0)
lorse <- sqrt(sd1^2 + sd0^2)
u95 <- exp(lor + qnorm(0.975)*lorse)
l95 <- exp(lor - qnorm(0.975)*lorse)
or <- exp(lor)
label <- sprintf("\n%s/%s treatment\n%s/%s control\naddtrt=%s:%s",
zdat$actr, zdat$actn, zdat$conr, zdat$conn, study$addtrt, study$t)
data.frame(study, trt=zname, r1, r0, n1, n0, p1, p0, lor, or, l95, u95, label)
}
zdf <- rbind(zdfn("Zoledronic_acid_1_IV"),
zdfn("Zoledronic_acid_2_IV")) %>% arrange(trt, study)
p <- ggplot(zdf, aes(x=study, y=or, label=label, ymin=l95, ymax=u95)) +
coord_flip() +
facet_wrap(. ~ trt) +
geom_hline(aes(yintercept=1), col="gray") +
geom_pointrange() +
scale_y_continuous(trans="log", breaks=c(0.1,0.2,0.5,1,2,3,5,10,20,100)) +
ylab("Ratio in odds of fever, vs observation only") +
xlab("")
ggplotly(p, tooltip="label")
## ---- eval=FALSE---------------------------------------------------------
#
# load_all("../../gemtc/gemtc")
#
# ## split single node
# mod.trtsplit <- mtc.model(net.trt, type="nodesplitgrouptreat", likelihood="cjbinom", link="logit", om.scale=scale, hy.prior=hy.prior, t1="210", t2="220")
# cat(mod.trtsplit$code)
# fit.trtsplit <- mtc.run(mod.trtsplit)
#
# bayesplot::mcmc_trace(fit.trtsplit[["samples"]], pars=c("d.direct", "d.indirect"))
# bayesplot::mcmc_areas_ridges(fit.trtsplit[["samples"]], pars=c("d.direct", "d.indirect"))
#
# ## split all nodes with direct and indirect evidence
# nsc <- mtc.nodesplit.comparisons(net.trt)
# modall.trtsplit <- mtc.nodesplitgrouptreat(net.trt, comparisons=nsc, likelihood="cjbinom", om.scale=scale, hy.prior=hy.prior, link="logit")
# modall.trtsplit
# plot(modall.trtsplit)
# summary(modall.trtsplit)
# summ <- summary(modall.trtsplit)
#
# summ$cons.effect
# summ$p.value
# plot(summ)
#
# ## Examine the direct vs indirect OR itself
# d.direct <- sapply(modall.trtsplit[1:9], function(x){unlist(x$samples[,"d.direct"])})
# d.indirect <- sapply(modall.trtsplit[1:9], function(x){unlist(x$samples[,"d.indirect"])})
# dvsind <- apply(d.direct - d.indirect, 2, quantile, c(0.025, 0.5, 0.975))
# dvsind <- as.data.frame(t(exp(dvsind)))
# names(dvsind) <- c("l95","est","u95")
# dvsind$comp <- rownames(dvsind)
#
# ggplot(dvsind, aes(x=comp, y=est, ymin=l95, ymax=u95)) +
# geom_pointrange() +
# scale_y_continuous(trans="log", breaks=c(0.1, 0.5, 1, 2, 5, 10)) +
# ylab("Odds ratio (direct / indirect)") + xlab("") +
# coord_flip()
#
# dicres <- t(sapply(modall.trtsplit,
# function(x){x$deviance[c("Dbar","pD","DIC")]}))
# dicres # small DIC improvements from some consistency models
# ## could also assess contribution of points to DIC
#
#
## ------------------------------------------------------------------------
## Priors
## * LOR scale 5. default data-based one 4.19 = max observed LOR
## * prior SD for LOR set to 15 times scale.
## * prior SD for reg coeffs set to scale.
## * prior for REs on effects.
## * logistic uniform for baseline log odds.
## * cov effs for baseline? same as for drug effs, 15scale # # baseline REs: t_4(0,1) as in gelman wiki. N(0,1)T(0,) allows between study SD for baseline log OR between 0.03, 2.2. qlogis(c(0.05, 0.95)) is -3 to 3 on log OR scale. would be SD of 2.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.