RCT_NRSanalysis.R

In htx-r/GenericModelNMA:

library(devtools)
library(R2jags)
install_github('htx-r/GenericModelNMA',force = T)
library(GenericModelNMA)

# load data

source('cleanData.R') # contains the FOUR different cleaned datasets

# data
# RCT-IPD
rct.ipd <-MSrelapse
# RCT-AD
rct.ad<-rct.ad[which(rct.ad$study=="Bornstein" | rct.ad$study=="Johnson"),]

# NRS-IPD
nrs.ipd <-nrs.ipd

# RCT-IPD-AD NMR
## *** # 3 #  RCT-IPD-AD NMR - all studies
# rct.ipd <- rct.ipd[!is.na(rct.ipd$EDSSBL),]
# rct.ipd$AGE <- rct.ipd$AGE-18
# rct.ipd$EDSSBL <- rct.ipd$EDSSBL*10
## *** # 4 #  RCT-IPD-AD NMR - only DEFINE as IPD and Johnson as AD
jagsdataIPDADnetmetareg<- with(rct.ipd,list(
  nIPD=3,
  np=nrow(rct.ipd),
  studyid=as.numeric(factor(STUDYID)),
  y=as.numeric(RELAPSE2year)-1,
  x=as.numeric(AGE),
  #xbar=unlist(sapply(unique(rct.ipd$STUDYID),function(i) rep(round(mean(rct.ipd[rct.ipd$STUDYID==i,'AGE'],na.rm = T),1),nrow(rct.ipd[rct.ipd$STUDYID==i,])))),
  t= rbind(c(4,1,NA),c(4,1,2),c(4,3,NA), c(4,2,NA),c(4,2,NA)),
  na=c(2,3,2,2,2),
  treat=as.numeric(factor(TRT01A)),
  baseline=rep(4,nrow(rct.ipd)),
  ref=4,
  nt=4,
  nAD=2,
  r=rbind(c(NA,19,NA,11),c(NA,97,NA,89)),
  n=rbind(c(NA,25,NA,25),c(NA,126,NA,125)),
  xbar.a =c(34.3,30)
)
)

jagsmodelIPDADnetmetareg <- jags.parallel(data = jagsdataIPDADnetmetareg,inits=NULL,parameters.to.save = c('d','tau','b','LOR'),model.file = modelIPDADnetmetareg,
                                          n.chains=3,n.iter = 10000,n.burnin = 4000,DIC=F,n.thin = 1)

# NRS-IPD NMR only

jagsdataIPDADnetmetaregNRS<- with(nrs.ipd,list(
  nIPD.NRS=1,
  np.NRS=nrow(nrs.ipd),
  studyid.NRS=as.numeric(factor(STUDYID)),
  y.NRS=as.numeric(RELAPSE2year)-1,
  x.NRS=as.numeric(AGE),
  t.NRS= rbind(c(4,1,2,3)),
  na.NRS=c(4),
  treat.NRS=as.numeric(factor(TRT01A)),
  baseline.NRS=rep(4,nrow(nrs.ipd)),
  ref.NRS=4,
  nt.NRS=4
)
)

jagsmodelIPDnetmetaregNRS <- jags.parallel(data = jagsdataIPDADnetmetaregNRS,inits=NULL,parameters.to.save = c('d.NRS'),model.file = modelIPDnetmetaregNRS,
                                        n.chains=3,n.iter = 10000,n.burnin = 2000,DIC=F,n.thin = 1)
jagsmodelIPDnetmetaregNRS

# RCT-NRS-IPD-AD-NMA (naive approach)
rct.ipd2 <- rct.ipd[,c('STUDYID','RELAPSE2year','AGE','TRT01A')]
rct.nrs.ipd <- rbind.data.frame(rct.ipd2,nrs.ipd)


jagsdataIPDADnetmetaregNRSnaive<- with(rct.nrs.ipd,list(
  nIPD=4,
  np=nrow(rct.nrs.ipd),
  studyid=as.numeric(factor(STUDYID)),
  y=as.numeric(RELAPSE2year)-1,
  x=as.numeric(AGE),
  t= rbind(c(4,1,NA,NA),c(4,1,2,NA),c(4,3,NA,NA), c(4,1,2,3),c(4,2,NA,NA),c(4,2,NA,NA)),
  na=c(2,3,2,4,2,2),
  treat=as.numeric(factor(TRT01A)),
  baseline=rep(4,nrow(rct.nrs.ipd)),
  ref=4,
  nt=4,
  nAD=2,
  r=rbind(c(NA,19,NA,11),c(NA,97,NA,89)),
  n=rbind(c(NA,25,NA,25),c(NA,126,NA,125)),
  xbar.a =c(34.3,30) # age
)
)
jagsmodelIPDADnetmetaregNRSnaive <- jags.parallel(data = jagsdataIPDADnetmetaregNRSnaive,inits=NULL,parameters.to.save = c('d','tau','bb','bw','LOR'),model.file = modelIPDADnetmetaregNRSnaive,
                                                  n.chains=3,n.iter = 10000,n.burnin = 2000,DIC=F,n.thin = 1)
jagsmodelIPDADnetmetaregNRSnaive


# RCT-NRS-IPD-AD NMA (design-adjusted)
jagsdataIPDADnetmetaregNRSdesAdj<- with(rct.nrs.ipd,list(
  nIPD=4,
  np=nrow(rct.nrs.ipd),
  studyid=as.numeric(factor(STUDYID)),
  y=as.numeric(RELAPSE2year)-1,
  x=as.numeric(AGE),
  t= rbind(c(4,1,NA,NA),c(4,1,2,NA),c(4,3,NA,NA), c(4,1,2,3),c(4,2,NA,NA),c(4,2,NA,NA)),
  na=c(2,3,2,4,2,2),
  treat=as.numeric(factor(TRT01A)),
  baseline=rep(4,nrow(rct.nrs.ipd)),
  bias=recode(STUDYID,'NRS'=1,.default=0),
  ref=4,
  nt=4,
  nAD=2,
  r=rbind(c(NA,19,NA,11),c(NA,97,NA,89)),
  n=rbind(c(NA,25,NA,25),c(NA,126,NA,125)),
  xbar.a =c(34.3,30)
)
)
jagsmodelIPDADnetmetaregNRSdesAdj <- jags.parallel(data = jagsdataIPDADnetmetaregNRSdesAdj,inits=NULL,parameters.to.save = c('d','tau','b','LOR','gamma'),model.file = modelIPDADnetmetaregNRSdesAdj,
                                                   n.chains=3,n.iter = 10000,n.burnin = 2000,DIC=F,n.thin = 1)
jagsmodelIPDADnetmetaregNRSdesAdj$BUGSoutput$summary



# RCT-NRS-IPD-AD-NMA (prior)
# NRS-IPD only
jagsdataIPDnetmetaNRS<- with(nrs.ipd,list(
  nIPD.NRS=1,
  np.NRS=nrow(nrs.ipd),
  studyid.NRS=as.numeric(factor(STUDYID)),
  y.NRS=as.numeric(RELAPSE2year)-1,
  t.NRS= rbind(c(4,1,2,3)),
  na.NRS=c(4),
  treat.NRS=as.numeric(factor(TRT01A)),
  baseline.NRS=rep(4,nrow(nrs.ipd)),
  ref.NRS=4,
  nt.NRS=4
)
)
jagsmodelIPDnetmetaNRS <- jags.parallel(data = jagsdataIPDnetmetaNRS,inits=NULL,parameters.to.save = c('d.NRS'),model.file = modelIPDnetmetaNRS,
                                                  n.chains=3,n.iter = 10000,n.burnin = 2000,DIC=F,n.thin = 1)
jagsmodelIPDnetmetaNRS

jagsdataIPDADnetmetaregRCT<- with(rct.ipd,list(
  nIPD=3,
  np=nrow(rct.ipd),
  studyid=as.numeric(factor(STUDYID)),
  y=as.numeric(RELAPSE2year)-1,
  x=as.numeric(AGE),
  t= rbind(c(4,1,NA),c(4,1,2),c(4,3,NA), c(4,2,NA),c(4,2,NA)),
  na=c(2,3,2,2,2),
  treat=as.numeric(factor(TRT01A)),
  baseline=rep(4,nrow(rct.ipd)),
  ref=4,
  nt=4,
  nAD=2,
  r=rbind(c(NA,19,NA,11),c(NA,97,NA,89)),
  n=rbind(c(NA,25,NA,25),c(NA,126,NA,125)),
  xbar.a =c(34.3,30),
  d.NRS=jagsmodelIPDnetmetaNRS$BUGSoutput$summary[,'mean'],
  var.NRS=jagsmodelIPDnetmetaNRS$BUGSoutput$summary[,'sd']
)
)

jagsmodelIPDADnetmetaregNRSprior <- jags.parallel(data = jagsdataIPDADnetmetaregRCT,inits=NULL,parameters.to.save = c('d','tau','bb','bw','LOR'),model.file = modelIPDADnetmetaregNRSprior,
                                                  n.chains=3,n.iter = 10000,n.burnin = 2000,DIC=F,n.thin = 1)
jagsmodelIPDADnetmetaregNRSprior



# 3. forest plot of relative treatment effect
# 'only RCT','only NRS','RCT-NRS:naive','RCT-NRS:prior','RCT-NRS:adjsuted design'
    x1 <- jagsmodelIPDADnetmetareg$BUGSoutput$summary
    x2 <- jagsmodelIPDnetmetaregNRS$BUGSoutput$summary
    x3 <- jagsmodelIPDADnetmetaregNRSnaive$BUGSoutput$summary
    x4 <- jagsmodelIPDADnetmetaregNRSprior$BUGSoutput$summary
    x5 <- jagsmodelIPDADnetmetaregNRSdesAdj$BUGSoutput$summary

    x<- rbind(x1,x2,x3,x4,x5)
    save(x,file = 'RCT-NRSresults')
    load('RCT-NRSresults')
    param.lab <- c('d')
    drug.lab <- levels(factor(nrs.ipd$TRT01A))
    ndrugs <- length(drug.lab)
    ref.lab <- 'Placebo'
    plotdata <- x %>% t() %>% data.frame() %>% select(starts_with(param.lab))%>%t()%>% data.frame()
    ml_nmr_df <- as.data.frame(summary(rrms_fit_FE))[7:9,]
    ml_nmr_res <- data.frame(mean=ml_nmr_df$mean,sd=ml_nmr_df$sd,
                        'X2.5.'=ml_nmr_df$'2.5%','X25.'=ml_nmr_df$'25%','X50.'=ml_nmr_df$'50%','X75.'=ml_nmr_df$'75%','X97.5.'=ml_nmr_df$'97.5%',
                        Rhat=ml_nmr_df$Rhat, n.eff=ml_nmr_df$"Bulk_ESS")
    plotdata <- rbind.data.frame(plotdata,ml_nmr_res)


    plotdata$drug <- c((rep(drug.lab,times=5)),c("Dimethyl fumarate"  ,"Glatiramer acetate", "Natalizumab"))
    plotdata <- plotdata[plotdata[,'drug']!=ref.lab,]
    plotdata$method <- rep(c('1. only RCT','2. only NRS','3. RCT-NRS naive','4. RCT-NRS as a prior','5. RCT-NRS adjusted design','6. ML-NMR'),each=3)

    plotdata$id <- 1:18



ggplot(aes(x = mean, y = id, col = method, shape = method), data = plotdata)+
  geom_vline(xintercept = 0, lty = 2) +
  geom_point(size = 2) +
  geom_segment(aes(y = id, yend = id, x = X2.5., xend = X97.5.), na.rm = TRUE) +
  xlab("Estimate (Log OR)") +
  facet_grid(drug~., switch = "y", scales = "free_y", space = "free_y") +
  scale_y_reverse(name = "", breaks = NULL, expand = c(0, 0.6))



table(rct.nrs.ipd$STUDYID,rct.nrs.ipd$TRT01A)
nrow(rct.ipd[rct.ipd$STUDYID=='DEFINE',])

## log(OR) over age for each drug
######
d1 <- jagsmodelIPDADnetmetaregNRSdesAdj$BUGSoutput$summary['d[1]','mean']
b1 <- jagsmodelIPDADnetmetaregNRSdesAdj$BUGSoutput$summary['b[1]','mean']
age <- seq(18,60, l=1000)
OR1 <- exp(d1+b1*(age-37))
plot(age,OR1,type = 'l',ylim=c(0.2,2))

d2 <- jagsmodelIPDADnetmetaregNRSdesAdj$BUGSoutput$summary['d[2]','mean']
b2 <- jagsmodelIPDADnetmetaregNRSdesAdj$BUGSoutput$summary['b[2]','mean']
age <- seq(18,60, l=1000)
OR2 <- exp(d2+b2*(age-37))
lines(age,OR2,type = 'l',col=2)

d3 <- jagsmodelIPDADnetmetaregNRSdesAdj$BUGSoutput$summary['d[3]','mean']
b3 <- jagsmodelIPDADnetmetaregNRSdesAdj$BUGSoutput$summary['b[3]','mean']
age <- seq(18,60, l=1000)
OR3 <- exp(d3+b3*(age-37))
lines(age,OR3,type = 'l',col=3)

plotdat <- data_frame(age=rep(age,time=3),
        OR=c(OR1,OR2,OR3),
        drug=rep(c("Dimethyl fumarate","Glatiramer acetate","Natalizumab"),each=1000)
        )
ggplot(aes(x = age, y = OR, col = drug), data = plotdat) +
  geom_line()