referencesNMA/otherScripts/longrun.R

In htx-r/GenericModelNMA:

## This file make a simulation over 100 datasets to compute its mean and compare it with true.mu=0.4
     # 1. IPDmeta (Bayes+Freq) 2. IPDmetareg(Bayes+Freq) 3. IPDADmeta 4. IPDADmetareg

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



# 1. IPDmeta (Bayes+Freq)
res <- list()
b <- f <- vector()
n.sim.data <- 100
start <- Sys.time()
for (j in 1: n.sim.data) {
#myf <- function(j){
 # Simulated data
  sim.data <- SimulatedIPDmeta(N.IPD=5)

  # Bayesian
  IPDmetaJAGSmodel <- jags.parallel(data = MakeJAGSdataIPDmeta(sim.data),inits=NULL,parameters.to.save = c('mu','tau.sq'),model.file = ModelIPDmeta,
                                    n.chains=2,n.iter = 10000,n.burnin = 1000,DIC=F,n.thin = 1)

  N.IPD <- length(unique(sim.data$study.IPD))
  results <- as.data.frame(array(NA, dim=c(N.IPD,2)))
  colnames(results) <- c("thetai", "var_thetai")

  for (i in 1:N.IPD) {
    sim.datai <- as.data.frame(sim.data[sim.data$study.IPD==i,])
    fit <- glm(Y_ij~treat, data=sim.datai, family=binomial(link = 'logit'))
    results[i,] <- c(coefficients(fit)[2], vcov(fit)[2,2])
  }

  # Stage 2
  fit <- metagen(thetai, var_thetai, data=results)
  # compare Bayesian and Frequentist
  f[j] <- fit$TE.random
  b[j] <- IPDmetaJAGSmodel$BUGSoutput$mean$mu

  res[[j]] <- c(b[j], f[j])

}

end <- Sys.time()
(end-start) # it will take around 36 minutes

## 2nd option: I tried to use sapply function hoping to have a faster simulation.
# start <- Sys.time()
# re <- sapply(1:3, myf)
# end <- Sys.time()
# (end-start)*33 # it will take around 42 minutes

resulIPDtmeta <- apply(matrix(unlist(res),n.sim.data,2,byrow = T),2,mean)
## 50 sim: -0.03734712 -0.02431448
## other 100 sim: -0.008947872  0.003909973





# 2. IPDmetareg: Bayes+Freq
res <- list()
b <- f <- vector()
n.sim.data <- 100
start<- Sys.time()
for (j in 1: n.sim.data) {
  # Simulated data
  sim.data <- SimulatedIPDmetareg(N.IPD=5)

  # Bayesian
  IPDmetaregJAGSmodel <- jags.parallel(data = MakeJAGSdataIPDmetareg(sim.data),inits=NULL,parameters.to.save = c('mu','tau.sq'),model.file = ModelIPDmetareg,
                                       n.chains=2,n.iter = 10000,n.burnin = 1000,DIC=F,n.thin = 1)

  N.IPD <- length(unique(sim.data$study.IPD))
  results <- as.data.frame(array(NA, dim=c(N.IPD,2)))
  colnames(results) <- c("thetai", "var_thetai")

  for (i in 1:N.IPD) {
    sim.datai <- as.data.frame(sim.data[sim.data$study.IPD==i,])
    fit <- glm(Y_ij~treat+as.factor(x_ij), data=sim.datai, family=binomial(link = 'logit'))
    results[i,] <- c(coefficients(fit)[2], vcov(fit)[2,2])
  }

  # Stage 2
  fit <- metagen(thetai, var_thetai, data=results)
  # compare Bayesian and Frequentist
  f[j] <- fit$TE.random
  b[j] <- IPDmetaregJAGSmodel$BUGSoutput$mean$mu

  res[[j]] <- c(b[j], f[j])

}
end <- Sys.time()
(end-start) # it will take around 53 minutes 48 mins
resultIPDmetareg <- apply(matrix(unlist(res),n.sim.data,2,byrow = T),2,mean)
## For 50 simulations:  -0.4206411  0.9744341
## For 100 simulations: -0.4270961  0.9757987 , for Bayesian it is ne


# 3. IPDADmeta
#res <- list()
b <- vector()
n.sim.data <- 100
start <- Sys.time()
for (j in 1: n.sim.data) {
#myf <- function(j){
sim.data <- SimulatedIPDADmeta(N.IPD=5)
  # Bayesian
  IPDADmetaJAGSmodel <- jags.parallel(data = MakeJAGSdataIPDADmeta(sim.data$IPDdata,sim.data$ADdata),inits=NULL,parameters.to.save = c('mu','tau.sq'),model.file = ModelIPDADmeta,
                                    n.chains=2,n.iter = 10000,n.burnin = 1000,DIC=F,n.thin = 5)

  b[j] <- IPDADmetaJAGSmodel$BUGSoutput$mean$mu

  #res[[j]] <- c(b[j], f[j])

}

end <- Sys.time()
(end-start) # it will take around 1 hour if I use for loops :o

## 2nd option: I tried to use sapply function hoping to have a faster simulation.
# start <- Sys.time()
# re <- sapply(1:2, myf)
# end <- Sys.time()
# (end-start)*50 # it will take around slightly more than 1 hour (3 mins more)! which is surprising me!!!

resulIPDADtmeta <- mean(b) ##  0.1073344






# 4. IPDmetareg
b <- vector()
n.sim.data <- 100
start<- Sys.time()
for (j in 1: n.sim.data) {
  # Simulated data
  sim.data <- SimulatedIPDADmetareg(N.IPD=5)

  # Bayesian
  IPDADmetaregJAGSmodel <- jags.parallel(data = MakeJAGSdataIPDADmetareg(sim.data$IPDdata,sim.data$ADdata),inits=NULL,parameters.to.save = c('mu','tau.sq'),model.file = ModelIPDADmetareg,
                                       n.chains=2,n.iter = 10000,n.burnin = 1000,DIC=F,n.thin = 1)

  b[j] <- IPDmetaregJAGSmodel$BUGSoutput$mean$mu

}
end <- Sys.time()
(end-start) # it  takes around 1 hour
resultIPDADmetareg <- mean(b) ## -0.4130127