In chjackson/adverse: Data on adverse events of adjuvant breast cancer treatment

knitr::opts_chunk$set(echo=FALSE,
                      fig.width=10,
                      fig.height=7)

Raw data

# library(adverse)
load_all("..")
nstudies <- length(unique(bpaearmtype$`Trial name`))
nstudies_fever <- length(unique(bpaearmtype[bpaearmtype$aetype == "FEVER",]$"Trial name"))

r {nstudies_fever} out of r {nstudies} studies reported an outcome of Fever. All grades combined.

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")

Treatment groupings for analysis

Consider four alternative groupings from finer to coarser:

1. Bisphosphonate drug, dose and delivery method. IV placebo, oral placebo and observation only distinguished.

2. Bisphosphonate drug and delivery method.

3. Bisphosphonate drug. IV and oral placebo grouped, but distinguished from observation only.

4. Drug class, considering all nitrogenous bisphosphonates together. IV and oral placebo grouped, but distinguished from observation only.

Notes

• Only study including Clodronate was NSABP B-34, whose control arm used oral placebo. 1 and 0 people with fever were reported in the active / control arms respectively, out of 1623, 1612. No other study used oral placebos. Exclude for now, but may be included later if we decide to combine all placebos or all controls.

• Two studies, ABCSG12 and HOBOE, included separate arms to compare a different kind of hormone therapy for the same bisphosphonate ("main" treatment). To start with, define these as different treatments.

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

Networks of comparisons

Number of studies for each comparison indicated on the connecting lines.

FIXME some studies missing in new data [1] "NaTaN" "ProBONE II" "Z-FAST" These are dropped because they only have one arm [which reported AEs?] fixme bit they still have data in both arms.

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")

Network meta analysis models

Model assumptions

• Relative effects for each comparison are "random effects", normally distributed between studies with an unknown mean and variance

• Baseline rate of adverse events (that is, for "observation only" controls) is also a random effect. Log odds is normally distributed between studies with an unknown mean and variance.

• gemtc R package used, with code modified to deal with random baselines. Package default weakly informative priors.

The only relative effects which can't be estimated from the data are those relating to Ibandronate. There is only one study for this drug (Body 2007) which reported data for fever. 26 out of 137 in the zoledronic acid arm reported fever, compared to 0 / 137 in the ibandronate arm. Unsure how to deal with this! Is it clinically plausible? We might be able to moderate this outlying result using some model which borrows information from other drugs which are similar to ibandronate? Or just merge all nitrogenous bisphosphonates?

#load_all("..")
load_all("../../gemtc/gemtc")
#unload("../../gemtc/gemtc")

hy.prior <- mtc.hy.empirical.lor(outcome.type="semi-objective",
                                 comparison.type="pharma-pharma")
scale <- 5
mod.trt <- mtc.model(net.trt, type="grouptreat", likelihood="cjbinom", link="logit", om.scale=scale, hy.prior=hy.prior, basetrt="103")
mod.drugdm <- mtc.model(net.drugdm, type="grouptreat", likelihood="cjbinom", link="logit", om.scale=scale, hy.prior=hy.prior, basetrt="Observation")
mod.drug <- mtc.model(net.drug, type="grouptreat", likelihood="cjbinom", link="logit", om.scale=scale, hy.prior=hy.prior, basetrt="Observation")
mod.drugclass <- mtc.model(net.drugclass, type="grouptreat", likelihood="cjbinom", link="logit", om.scale=scale, hy.prior=hy.prior, basetrt="Observation")

fit.trt <- mtc.run(mod.trt)
fit.drugdm <- mtc.run(mod.drugdm)
fit.drug <- mtc.run(mod.drug)
fit.drugclass <- mtc.run(mod.drugclass)

trace.basic(fit.trt)
trace.basic(fit.drug)
trace.basic(fit.drugdm)
trace.basic(fit.drugclass)

Relative treatment effects

Estimates from each of the four network meta-analysis models, compared with direct data, for the comparisons where there is direct data.

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 

• Lower fever rates in Zoledronic acid dose definition 4 compared to other Zoledronic acid doses. Makes sense as this is "delayed" dosage?

• Not much difference between Zoledronic acid dose definitions 1, 2, 3

• Apparent difference between "placebo" and "observation only" controls when compared to Zoledronic acid dose 2. Sounds implausible. Consistent with unexplained between study heterogeneity

• Perhaps this causes the apparent overestimate (compared to direct data) of the effect of Zol 1 IV and Zol 2 IV (non-AI hormone) vs observation only.

• Pamidronate dose 8 - not much difference from Zoledronic acid dose def 1

• Pamidronate dose 1 appears to do better than Zol 1 (hence better than Pamidronate dose 8)

• Denosumab appears to do better than Zoledronate dose 1, but this effect is moderated under the model

• Note estimated relative odds is exactly 1 between e.g. two drugs in the same class under the "drug class" model, or two dosages or delivery methods of the same drug in the "drug" model.

The four alternative ways of grouping effects of similar treatments appear to give a similar fit to the data, according to standard methods of statistical model comparison (DIC)

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

Absolute rates of fever

To produce the absolute event rate under each treatment:

• first the absolute rate under control (observation only) for a predicted new study is estimated

• then the relative odds for each treatment is applied to this

Posterior mean and 95% credible intervals

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)

• Observation vs placebo strangeness

• More uncertainty in estimates from model, compared to raw data, reflecting between study heterogeneity in baseline rates

Influence of additional treatment

The figure below illustrates the event rate under "observation only" for each study, with studies grouped in the plot according to the additional treatment(s) given. The network meta-analysis model was used, so that estimates could be obtained for studies which did not include an "observation only" control group.

So we could calculate predicted baseline event rates either by

• averaging over all additional treatments (as done above).

• modelling the contrast between "mixed", "hormone only" and/or "chemo only", then reporting the predicted rate for a specific additional treatment group. The estimated contrast between (mixed or chemo only) vs (hormone only) is "statistically significant" (according to DIC) but it's driven by a small number of data points.

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)

Influence of reporting quality

Similarly plot predicted event rate under control against reporting quality category.

No evidence of higher baseline event rates from studies with more complete reporting.

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)

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

Between-study heterogeneity in relative effects

Limited amount of information to investigate between-study heterogeneity in the relative treatment effects. There are 6 studies each which compared Zol 1 IV and observation, and Zol 2 IV and observation. Estimated odds ratios from these are plotted below.

Heterogeneity here seems driven by small counts, unlikely to be enough information to determine any predictors of heterogeneity.

Each other direct comparison has at most 2 studies (see network plots above).

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")

Consistency assumption

• Assumption that direct and indirect evidence are consistent: checked by "node splitting" method and conflict p-values (Dias et al, Stat Med 2010). No significant inconsistencies found.

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.