IPD_ADanalysis.R

In htx-r/GenericModelNMA:

#!!!!! Ask Georgia:  In Johnsan, the age and gender are reported on each arm - how to combine thier SD
#!!!!! SEX in the dataset - 1 is F or M
#!!!!! how to incorporate different studies
# A = Placebo
# B = Dimethyl fumarate
# C = Glatiramer acetate

# load libraries
library(sandwich)
library(dplyr)
library(tidyr)
library(devtools)
library(R2jags)
install_github('htx-r/GenericModelNMA',force = T)
library(GenericModelNMA)
#----- Prepare Data -------
# load data

source('cleanData.R') # contains the FOUR different datasets
# RCT-IPD
rct.ipd <-MSrelapse
# RCT-AD
rct.ad<-rct.ad[which(rct.ad$study=="Bornstein" | rct.ad$study=="Johnson"),]

# IPDdata - AB
rct.ipd_DEFINE <- rct.ipd[rct.ipd$STUDYID=='DEFINE',]
AB.IPD <-rct.ipd_DEFINE%>%
  with(data.frame(ID=1:nrow(rct.ipd_DEFINE),
                  age=as.integer(AGE),
                  gender=factor(SEX),
                  trt=as.character(recode(TRT01A,'Dimethyl fumarate' = 'B','Placebo'='A')),
                  y=as.integer(RELAPSE2year)-1)
       )

# ADdata - AC
# I took the data from Avilable aggregated data for MS relapse.docx - only for Johnson study

AC.AgD <- data.frame(age.mean=34.3,
                  age.sd=6.5,
                  N.male=30,
                  prop.male=0.24,
                  y.A.sum=97,
                  y.A.bar=97/126,
                  N.A=126,
                  y.C.sum=89,
                  y.C.bar=89/125,
                  N.C=125)

###########################
#  MAIC
###########################

#** 1.weights

# the objective function to minimize
objfn <- function(a1, X){
  sum(exp(X %*% a1))
}
gradfn <- function(a1, X){
  colSums(sweep(X, 1, exp(X %*% a1), "*"))
}

X.EM.0 <- sweep(with(AB.IPD, cbind(age, age^2)), 2,
                with(AC.AgD, c(age.mean, age.mean^2 + age.sd^2)), '-')
# find the optimal a1
opt1 <- optim(par = c(0,0), fn = objfn, gr = gradfn, X = X.EM.0, method = "BFGS")

# ... result
a1 <- opt1$par # the optimal a1
wt <- exp(X.EM.0 %*% a1)
# resclaed weight
N_AB <-nrow(AB.IPD)
wt.rs <- (wt / sum(wt)) * N_AB

#*** 2. estimate the indirect BC effect
# AB effect
fit1 <-
  AB.IPD %>% mutate(y0 = 1 - y, wt = wt) %>%
  glm(cbind(y,y0) ~ trt, data = ., family = binomial, weights = wt)
# Sandwich estimator of variance matrix
V.sw <- vcovHC(fit1)
# The log OR of B vs. A is just the trtB parameter estimate,
# since effect modifiers were centred
print(d.AB.MAIC <- coef(fit1)["trtB"])

print(var.d.AB.MAIC <- V.sw["trtB","trtB"])

# AC parameter
# Estimated log OR of C vs. A from the AC trial
d.AC <- with(AC.AgD, log(y.C.sum * (N.A - y.A.sum) / (y.A.sum * (N.C - y.C.sum))))
var.d.AC <- with(AC.AgD, 1/y.A.sum + 1/(N.A - y.A.sum) + 1/y.C.sum + 1/(N.C - y.C.sum))

# Indirect comparison: BC
print(d.BC.MAIC <- d.AC - d.AB.MAIC)
#
print(var.d.BC.MAIC <- var.d.AC + var.d.AB.MAIC)

#


###########################
#  STC implementation
###########################

AB.IPD$y0 <- 1 - AB.IPD$y # Add in dummy non-event column
# Fit binomial GLM
STC.GLM <- glm(cbind(y,y0) ~ trt*I(age - AC.AgD$age.mean),
               data = AB.IPD, family = binomial)
summary(STC.GLM)


# Try adding prognostic variables to improve model fit
add1(STC.GLM, ~.+gender, test="Chisq")

# does not improve the fit match so will leave gender out of the model

# estimate the AB effect in the AC population
print(d.AB.STC <- coef(STC.GLM)["trtB"])

# its variance
print(var.d.AB.STC <- vcov(STC.GLM)["trtB","trtB"])

## indirect comparison
print(d.BC.STC <- d.AC - d.AB.STC)
print(var.d.BC.STC <- var.d.AC + var.d.AB.STC)

##################
# Summary
##################

# unadjusted estimate of AB effect in AC population
# AB.IPD %>% group_by(trt) %>%
#   summarise(y.sum = sum(y)) %>%
#   spread(trt, y.sum) %>%
#   with({
#     d.AB.AB <<- log(B * (N_AB/2 - A) / (A * (N_AB/2 - B)))
#     var.d.AB.AB <<- 1/B + 1/(N_AB/2 - A) + 1/A + 1/(N_AB/2 - B)
#   })
#
# # unadjusted BC estimated effect
# d.BC.NAIVE <- d.AC - d.AB.AB
# var.d.BC.NAIVE <- var.d.AC + var.d.AB.AB

jagsdataIPDADnetmeta<- with(rct.ipd,list(
  nIPD=3,
  np=nrow(rct.ipd),
  studyid=as.numeric(factor(STUDYID)),
  y=as.numeric(RELAPSE2year)-1,
  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,11,NA,19),c(NA,89,NA,97)),
  n=rbind(c(NA,25,NA,25),c(NA,125,NA,126))
)
)

jagsmodelIPDADnetmeta <- jags.parallel(data = jagsdataIPDADnetmeta,inits=NULL,parameters.to.save = c('d','tau','LOR'),model.file = modelIPDADnetmeta,
                                       n.chains=3,n.iter = 10000,n.burnin = 2000,DIC=F,n.thin = 1)
jagsmodelIPDADnetmeta
# unadjusted BC estimated effect
d.AC.NAIVE <- -jagsmodelIPDADnetmeta$BUGSoutput$summary['LOR[4,1]','mean']
var.d.AC.NAIVE<- jagsmodelIPDADnetmeta$BUGSoutput$summary['LOR[4,1]','sd']^2

d.AB.NAIVE <- -jagsmodelIPDADnetmeta$BUGSoutput$summary['LOR[4,2]','mean']
var.d.AB.NAIVE<- jagsmodelIPDADnetmeta$BUGSoutput$summary['LOR[4,2]','sd']^2

d.BC.NAIVE <- jagsmodelIPDADnetmeta$BUGSoutput$summary['LOR[2,1]','mean']
var.d.BC.NAIVE<- jagsmodelIPDADnetmeta$BUGSoutput$summary['LOR[2,1]','sd']^2

# d.BC.NAIVE <- d.AC - d.AB.AB
# var.d.BC.NAIVE <- var.d.AC + var.d.AB.AB

#######
# Saramago
######
## *** # 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)
# Estimates
d.AB.saramago <- -jagsmodelIPDADnetmetareg$BUGSoutput$summary['LOR[4,1]','mean']
var.d.AB.saramago<- jagsmodelIPDADnetmetareg$BUGSoutput$summary['LOR[4,1]','sd']^2

d.AC.saramago <- -jagsmodelIPDADnetmetareg$BUGSoutput$summary['LOR[4,2]','mean']
var.d.AC.saramago<- jagsmodelIPDADnetmetareg$BUGSoutput$summary['LOR[4,2]','sd']^2

d.BC.saramago <- jagsmodelIPDADnetmetareg$BUGSoutput$summary['LOR[2,1]','mean']
var.d.BC.saramago<- jagsmodelIPDADnetmetareg$BUGSoutput$summary['LOR[2,1]','sd']^2

#######
# Plot
######
plotdat <- data_frame(
  id = 1:10,
  Comparison = factor(c(rep(1,4), 2, 2,rep(3,4)),
                      labels = c("Dimethyl. vs. Placebo", "Glatiramer. vs. Placebo", "Glatiramer. vs. Dimethyl.")),
  Estimate = c( d.AB.saramago,d.AB.MAIC, d.AB.STC, d.AB.NAIVE,
                d.AC.saramago,d.AC.NAIVE,
                d.BC.saramago,d.BC.MAIC, d.BC.STC, d.BC.NAIVE),
  var = c( var.d.AB.saramago,var.d.AB.MAIC, var.d.AB.STC, var.d.AB.NAIVE,
           var.d.AC.saramago,var.d.AC.NAIVE,
           var.d.BC.saramago,var.d.BC.MAIC, var.d.BC.STC, var.d.BC.NAIVE),
  lo = Estimate + qnorm(0.025) * sqrt(var),
  hi = Estimate + qnorm(0.975) * sqrt(var),
  method = c( '3LH-MNR',"MAIC", "STC", "Unadjusted",
            '3LH-MNR',"Unadjusted",
            '3LH-MNR',"MAIC", "STC", "Unadjusted")
)
ggplot(aes(x = Estimate, y = id, col = method, shape = method), data = plotdat) +
  geom_vline(xintercept = 0, lty = 2) +
  geom_point(size = 2) +
  geom_segment(aes(y = id, yend = id, x = lo, xend = hi), na.rm = TRUE) +
  xlab("Estimate (Log OR)") +
  facet_grid(Comparison~., switch = "y", scales = "free_y", space = "free_y") +
  scale_y_reverse(name = "", breaks = NULL, expand = c(0, 0.6))



##############
# Multinma implemntation
##############
library(multinma)
library(dplyr)
library(tidyr)
library(ggplot2)
options(mc.cores = parallel::detectCores())

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

# load data
source('cleanData.R') # contains the FOUR different datasets
# prepare the data
# rct.ipd <- rct.ipd %>% mutate(
# complete = with(rct.ipd,complete.cases(TRT01A, RELAPSE2year, AGE)))
# rct.ipd <- filter(rct.ipd, complete)
rct.ipd <-MSrelapse
rct.ad<-rct.ad[which(rct.ad$study=="Bornstein" | rct.ad$study=="Johnson"),]

rct.ipd$RELAPSE2year <- as.numeric(rct.ipd$RELAPSE2year)-1
rct.ad$age_mean <- c(34.3,34.3,30,30)
rct.ad$age_sd <- c(6.5,6.5,6,6)

#rct.ipd$TRT01A <- factor(rct.ipd$TRT01A)
# set the network
rrms_net <- combine_network(
  set_ipd(rct.ipd,
          study = STUDYID,
          trt = TRT01A,
          r = RELAPSE2year,
          trt_ref = "Placebo"),
  set_agd_arm(rct.ad,
              study = study,
              trt = treat,
              r = r,
              n = n,
              trt_ref = "Placebo")
)
# plot the network
plot(rrms_net, weight_nodes = TRUE, weight_edges = TRUE) +
  ggplot2::theme(legend.position = "bottom", legend.box = "vertical")

# integration
rrms_net <- add_integration(rrms_net,
                           AGE = distr(qnorm, mean = age_mean, sd = age_sd),
                           n_int = 1000
)


# fit multinma model
rrms_fit_FE <- nma(rrms_net,
                  trt_effects = "fixed",
                  link = "logit",
                  likelihood = "bernoulli2",
                  regression = ~.trt+AGE,
                  prior_intercept = normal(scale = 10),
                  prior_trt = normal(scale = 10),
                  prior_reg = normal(scale = 10),
                  init_r = 0.1,
                  QR = TRUE
                  )
#> Note: Setting "PBO" as the network reference treatment.

#
as.data.frame(summary(rrms_fit_FE))




# AGE distribution per study
library(dplyr)
library(ggplot2)
rct.ipd2 <- rct.ipd[,c('STUDYID','RELAPSE2year','AGE','TRT01A')]
rct.nrs.ipd <- rbind.data.frame(rct.ipd2,nrs.ipd)
rct.ad1 <- data.frame(STUDYID='Bornstein',AGE=rnorm(50,34.3,6.5))
rct.ad2 <- data.frame(STUDYID='Johnson',AGE=rnorm(251,30,6))

plotdat <- rbind.data.frame(rct.nrs.ipd[,c('STUDYID','AGE')],rct.ad1,rct.ad2)

plotdat %>%
  ggplot( aes(x=AGE, fill=STUDYID)) +
  geom_histogram( color="#e9ecef", alpha=0.6, position = 'identity',binwidth=2) +
  scale_fill_manual(values=c("#69b3a2", "#404080","salmon4","lightgreen",'red','blue')) +
  theme_minimal() +
  labs(fill="")+
  facet_wrap(~STUDYID)

#+
# geom_vline(xintercept = c(34.3,34.3+2*6.5,34.3-2*6.5), col='red',lwd=2,lty=c(1,2,2))+
# geom_vline(xintercept = c(30,30+2*6,30-2*6), col='green',lwd=2,lty=c(1,2,2))