referencesNMA/otherScripts/randomipd.R

In htx-r/GenericModelNMA:

library(devtools)
install_github("htx-r/CleaningData",force=TRUE)
install_github("htx-r/NMApredictionsRiskModel", force = TRUE)
install_github("htx-r/GenericModelNMA", force = TRUE)

library(NMApredictionsRiskModel)
library(CleaningData)
library(GenericModelNMA)
###### Give your path of data
mydatapath="/Users/tasneem/Desktop/TasnimPhD/multiple sclerosis/IPD data from 6 Biogen trials"

####################################  DATA   ###########################################
#########################################################################################
#### load data
cleanBIOGENtrials<-cleanBIOGENtrials.fun(mydatapath)
adsl01<-cleanBIOGENtrials$adsl01
ModelIPDADnetmeta <- function(){

  ######## NMA model
  for (i in 1:Np){
    dropout[i]~dbern(p.d[i] )
    logit(p.d[i])<-u0[studyid[i]]+(1-equals(treat[i],2))*Delta[i] + gamma[studyid[i]]*Risk[i] + (1-equals(treat[i],2))*Beta[i]*Risk[i]
    Beta[i]<-b[treat[i]]-b[2]
    Delta[i]<-d[treat[i]]-d[2]
  }
  #### reference treatment= CBASP+medication =1
  d[2]<-0
  b[2]<-0
  for (i in c(1,3,4)){
    d[i]~dnorm(thetad,precisiond)
    b[i]~dnorm(thetab,precisionb)
  }
  for (i in 1:Nstudies){
    u0[i]~dnorm(0,0.001)
    gamma[i]~dnorm(0,0.001)
  } # for prediction model MAKE SURE TO CHANGE THIS FOR SIMPLE IPD-NMA
  thetad~dnorm(0,0.001)
  taud~dunif(0,1000)
  tausqd<-taud^2
  precisiond<-1/tausqd

  thetab~dnorm(0,0.001)
  taub~dunif(0,1000)
  tausqb<-taub^2
  precisionb<-1/tausqb
}

adsl01_2 <- adsl01[adsl01$STUDYID%in%c('DEFINE','CONFIRM','AFFIRM'),]
#as.numeric(table(as.numeric(as.factor(adsl01_2$STUDYID))))
jagsdataIPDADnetmeta <- list(
  Nstudies=3,
  Npatients=as.numeric(table(as.numeric(as.factor(adsl01_2$STUDYID)))),
  Np=sum(as.numeric(table(as.numeric(as.factor(adsl01_2$STUDYID))))),
  studyid=as.numeric(as.factor(adsl01_2$STUDYID))-1,
  dropout=adsl01_2$RELAPSE2year,
  treat= as.numeric(droplevels(as.factor(adsl01_2$TRT01A))),
  Risk=rnorm(sum(as.numeric(table(as.numeric(as.factor(adsl01_2$STUDYID))))))

)
IPDADnetmetaJAGSmodel <- jags.parallel(data = jagsdataIPDADnetmeta ,inits=NULL,parameters.to.save = c('d','gamma','b'),model.file = ModelIPDADnetmeta,
                                       n.chains=2,n.iter = 100,n.burnin = 5,DIC=F,n.thin = 1)
str(jagsdataIPDADnetmeta)
## Freq  ------------------------------------------------------------------------
# Structure the dataset as the AcuteMania: Treatment name, r, n, studyid
##make the dataset proper for netmeta
adsl01_2$STUDYID<-as.numeric(as.factor(adsl01_2$STUDYID))-1
table(adsl01_2$STUDYID)
adsl01_2$TRT01A<-as.numeric(droplevels(as.factor(adsl01_2$TRT01A)))
table(adsl01_2$TRT01A)


STUDYID<-c(1,1,2,2,2,3,3)
Treatment<-c(1,4,1,2,4,3,4)
n<-c(826,408,703,351,363,627,312)
r<-c(212,182,185,117,149,183,176)
Data<-cbind(STUDYID,Treatment,n,r)
Data<-as.data.frame(Data)

library(metafor)
library(meta)
library(xlsx)
library(readxl)
library(netmeta)

AcuteManiaPair = pairwise(treat = Treatment, event = r, n = n,
                          data = Data, studlab = STUDYID, sm = "OR")
net1 = netmeta(AcuteManiaPair, ref = 2, comb.fixed = FALSE, comb.random = TRUE)


exp(IPDADnetmetaJAGSmodel$BUGSoutput$mean$d)
coef()