###Note: dplyr, Matrix, lme4, meta packages were needed.
# meta: number of meta-analysis or loops
# m: number of studies in each meta-analysis
# n: number of patients in control group in each individual study
# pc: incidence rate of rare AEs in control group
# or: adverse effect size
# size: number of trails in binomial distribution
# tau: between-study heterogeneity
# r: sample size ratio between treatment and control group
# seed: produce same random number for the same starting seed

# reset key parameters' values after finishing a scenario
meta <- 1000 
m <- 10
n <- 60   
pc <- 0.05 
or <- 1.25 
tau <- 0.1              
size <- 1                 
r <- 1        
set.seed(2018)
log_or=log(or)


# Generating simulated individual participants' binary data
# meta_id: indicates the id of each meta 
# paper_id: indicates the id of papers within each meta
# id: indicates the id of individual patients within each paper
# event: indicates if the outcome event happened
# out_data: individual participants¡¯ binary data
# tr_p¡¢tr_n¡¢ctr_p and ctr_n: indicate total number of events happened or not in the treatment and control group respectively



for (j in c(1:meta)) {
  
  for (i in c(1:m)) {  
    log_or_hat = rnorm (1, log_or, tau*tau)
    or_hat = exp(log_or_hat)
    pt=(pc/(1-pc)*or_hat)/(1+(pc/(1-pc))*or_hat)
    paper_id=c(rep(i,n))
    id =c(1:n)
    event = rbinom (n, size*r, pc)
    group <- rep(0,n)
    meta_id <- c(rep(j,n))
    ctr_data=data.frame(meta_id,paper_id,id,event,group)
    ctr_p1 <- nrow(ctr_data [ctr_data$event==1,])
    ctr_n1 <- nrow(ctr_data [ctr_data$event==0,])
    if (i==1) { out_ctr_data=ctr_data}& {pt_all <- pt}&{ctr_p <- ctr_p1}&{ctr_n <- ctr_n1} else
    {out_ctr_data=rbind(ctr_data, out_ctr_data)  }& {pt_all <- c(pt_all,pt)}&{ ctr_p <- c(ctr_p , ctr_p1)} & { ctr_n <- c(ctr_n, ctr_n1)}
  }
  
  
  
  for (i in c(1:m)) {  
    paper_id=c(rep(i,n))
    id =c(1:n)
    event = rbinom (n, size, pt_all[i])
    group <- rep(1,n)
    meta_id <- c(rep(j,n))
    tr_data=data.frame(meta_id,paper_id,id,event,group)
    tr_p1 <- nrow(tr_data [tr_data$event==1,])
    tr_n1 <- nrow(tr_data[tr_data$event==0,])
    if (i==1) { out_tr_data=tr_data} & {tr_p <- tr_p1} & {tr_n <- tr_n1}else { out_tr_data=rbind(tr_data, out_tr_data)} & {tr_p <- c(tr_p, tr_p1 )} & {tr_n <- c( tr_n, tr_n1 )}
    
  }
  tr_p<- matrix(tr_p,ncol=1)
  tr_n<- matrix(tr_n,ncol=1)
  ctr_p<- matrix(ctr_p,ncol = 1)
  ctr_n<- matrix(ctr_n,ncol = 1)
  
  paper <- c(1:m)
  metafirst <- data.frame(j,paper, tr_p, tr_n, ctr_p, ctr_n)
  out_data <- rbind(out_ctr_data, out_tr_data)
  
  # Analysis simulated data by IPD one_stage multilevel model

  m_onestage <- glmer(event ~ group + (1|paper_id),family = binomial(link = "logit"), data = out_data)
  beta <- summary(m_onestage)$coe[2]
  se <- summary(m_onestage)$coe[4]
  OR <- exp(beta)
  ORh <- exp(beta+1.96*se)
  ORl <- exp(beta-1.96*se)
  OR_d <- ORh - ORl
  cov <- ifelse(ORh>or & ORl<or,1,0)
  if (j==1) { OR_Ipd <- OR ; ORh_Ipd <- ORh; ORl_Ipd <- ORl; metatest <- metafirst;cov_Ipd <- cov;OR_D <- OR_d}else { OR_Ipd <- rbind(OR_Ipd,OR) ; ORh_Ipd <- rbind(ORh_Ipd,ORh); ORl_Ipd <- rbind(ORl_Ipd,ORl); metatest <- rbind(metatest,metafirst);cov_Ipd <- rbind(cov_Ipd,cov);OR_D <- rbind(OR_D,OR_d) }
  
}

# Analysis simulated data by IPD two_stage multilevel model

for (i in c(1:meta)) {
  Pt <- metabin(tr_p, tr_n, ctr_p, ctr_n,method="Peto" ,data=metatest[c((1+m*(i-1)):m*i),c(3:6)],sm="OR")
  MH <- metabin(tr_p, tr_n, ctr_p, ctr_n,method="MH" ,data=metatest[c((1+m*(i-1)):m*i),c(3:6)],sm="OR")
  IV <- metabin(tr_p, tr_n, ctr_p, ctr_n,method="Inverse" ,data=metatest[c((1+m*(i-1)):m*i),c(3:6)],sm="OR")
  OR_p1  <- exp(Pt$TE.fixed)
  OR_p1h <- exp(Pt$TE.fixed+1.96*Pt$seTE.fixed)
  OR_p1l <- exp(Pt$TE.fixed-1.96*Pt$seTE.fixed)
  OR_p1d <- OR_p1h - OR_p1l
  cov_p1 <- ifelse(OR_p1h > or & OR_p1l < or, 1 , 0)
 
  
  OR_mf  <- exp(MH$TE.fixed)
  OR_mfh <- exp(MH$TE.fixed+1.96*MH$seTE.fixed)
  OR_mfl <- exp(MH$TE.fixed-1.96*MH$seTE.fixed)
  OR_mfd <- OR_mfh - OR_mfl
  cov_mf <- ifelse(OR_mfh > or & OR_mfl < or, 1, 0)

  
  OR_mr <- exp(MH$TE.random)
  OR_mrh <- exp(MH$TE.random+1.96*MH$seTE.random)
  OR_mrl <- exp(MH$TE.random-1.96*MH$seTE.random)
  OR_mrd <- OR_mrh - OR_mrl
  cov_mr <- ifelse(OR_mrh > or & OR_mrl < or, 1, 0)
 
  
  OR_if  <- exp(IV$TE.fixed)
  OR_ifh <- exp(IV$TE.fixed+1.96*IV$seTE.fixed)
  OR_ifl <- exp(IV$TE.fixed-1.96*IV$seTE.fixed)
  OR_ifd <- OR_ifh - OR_ifl
  cov_if <- ifelse(OR_ifh > or & OR_ifl < or, 1, 0)

  
  OR_ir <- exp(IV$TE.random)
  OR_irh <- exp(IV$TE.random+1.96*IV$seTE.random)
  OR_irl <- exp(IV$TE.random-1.96*IV$seTE.random)
  OR_ird <- OR_irh - OR_irl
  cov_ir <- ifelse(OR_irh > or & OR_irl < or, 1, 0)

  
  
  if (i==1) { OR_peto <- OR_p1;OR_mhf <- OR_mf; OR_mhr <- OR_mr; OR_inf <- OR_if; OR_inr <- OR_ir;cov_Peto <- cov_p1;cov_Mf <- cov_mf;cov_Mr <- cov_mr;cov_If <- cov_if;cov_Ir <- cov_ir;OR_Pd <- OR_p1d;OR_Mfd <- OR_mfd;OR_Mrd <- OR_mrd;OR_Ifd <- OR_ifd;OR_Ird <- OR_ird}else { OR_peto <- c(OR_peto,OR_p1);OR_mhf <- c(OR_mhf,OR_mf);OR_mhr <- c(OR_mhr,OR_mr) ;OR_inf <- c(OR_inf,OR_if);OR_inr <- c(OR_inr,OR_ir);cov_Peto <- c(cov_Peto,cov_p1);cov_Mf <- c(cov_Mf,cov_mf);cov_Mr <- c(cov_Mr,cov_mr);cov_If <- c(cov_If,cov_if);cov_Ir <- rbind(cov_Ir,cov_ir);OR_Pd <- c(OR_Pd,OR_p1d);OR_Mfd <- c(OR_Mfd,OR_mfd);OR_Mrd <- rbind(OR_Mrd,OR_mrd);OR_Ifd <- c(OR_Ifd,OR_ifd);OR_Ird <- c(OR_Ird,OR_ird)}
  
}
cov_IPD <- sum(cov_Ipd[!is.na(cov_Ipd) & !is.nan(cov_Ipd) & !is.infinite(cov_Ipd)])/(meta-length(cov_Ipd[is.infinite(cov_Ipd)|is.na(cov_Ipd)|is.nan(cov_Ipd)]))
cov_PETO <- sum(cov_Peto[!is.na(cov_Peto) & !is.nan(cov_Peto) & !is.infinite(cov_Peto)])/(meta-length(cov_Peto[is.infinite(cov_Peto)|is.na(cov_Peto)|is.nan(cov_Peto)]))
cov_MF <- sum(cov_Mf[!is.na(cov_Mf) & !is.nan(cov_Mf) & !is.infinite(cov_Mf)])/(meta-length(cov_Mf[is.infinite(cov_Mf)|is.na(cov_Mf)|is.nan(cov_Mf)]))
cov_MR <- sum(cov_Mr[!is.na(cov_Mr) & !is.nan(cov_Mr) & !is.infinite(cov_Mr)])/(meta-length(cov_Mr[is.infinite(cov_Mr)|is.na(cov_Mr)|is.nan(cov_Mr)]))
cov_IF <- sum(cov_If[!is.na(cov_If) & !is.nan(cov_If) & !is.infinite(cov_If)])/(meta-length(cov_If[is.infinite(cov_If)|is.na(cov_If)|is.nan(cov_If)]))
cov_IR <- sum(cov_Ir[!is.na(cov_Ir) & !is.nan(cov_Ir) & !is.infinite(cov_Ir)])/(meta-length(cov_Ir[is.infinite(cov_Ir)|is.na(cov_Ir)|is.nan(cov_Ir)]))

OR_peto_m <- mean(OR_peto[!is.nan(OR_peto)&!is.na(OR_peto)&!is.infinite(OR_peto)])
OR_mhf_m <- mean(OR_mhf[!is.nan(OR_mhf)&!is.na(OR_mhf)&!is.infinite(OR_mhf)])
OR_mhr_m <- mean(OR_mhr[!is.nan(OR_mhr)&!is.na(OR_mhr)&!is.infinite(OR_mhr)])
OR_inf_m <- mean(OR_inf[!is.nan(OR_inf)&!is.na(OR_inf)&!is.infinite(OR_inf)])
OR_inr_m <- mean(OR_inr[!is.nan(OR_inr)&!is.na(OR_inr)&!is.infinite(OR_inr)])
OR_Ipd_m <- mean(OR_Ipd[!is.nan(OR_Ipd)&!is.na(OR_Ipd)&!is.infinite(OR_Ipd)])

OR_peto_percentb <- (OR_peto_m-or)/or
OR_mhf_percentb <- (OR_mhf_m-or)/or
OR_mhr_percentb <- (OR_mhr_m-or)/or
OR_inf_percentb <- (OR_inf_m-or)/or
OR_inr_percentb <- (OR_inr_m-or)/or
OR_Ipd_percentb <- (OR_Ipd_m-or)/or

OR_Pd_m <- mean(OR_Pd[!is.nan(OR_Pd)&!is.na(OR_Pd)&!is.infinite(OR_Pd)])
OR_Mfd_m <- mean(OR_Mfd[!is.nan(OR_Mfd)&!is.na(OR_Mfd)&!is.infinite(OR_Mfd)])
OR_Mrd_m <- mean(OR_Mrd[!is.nan(OR_Mrd)&!is.na(OR_Mrd)&!is.infinite(OR_Mrd)])
OR_Ifd_m <- mean(OR_Ifd[!is.nan(OR_Ifd)&!is.na(OR_Ifd)&!is.infinite(OR_Ifd)])
OR_Ird_m <- mean(OR_Ird[!is.nan(OR_Ird)&!is.na(OR_Ird)&!is.infinite(OR_Ird)])
OR_IPDd_m <-mean(OR_D[!is.nan(OR_D)&!is.na(OR_D)&!is.infinite(OR_D)])


MSE_peto <- mean((OR_peto[!is.na(OR_peto) & !is.nan(OR_peto) & !is.infinite(OR_peto)] - or)^2)
MSE_mhf <- mean((OR_mhf[!is.na(OR_mhf) & !is.nan(OR_mhf) & !is.infinite(OR_mhf)] - or)^2)
MSE_mhr <- mean((OR_mhr[!is.na(OR_mhr) & !is.nan(OR_mhr) & !is.infinite(OR_mhr)] - or)^2)
MSE_inf <- mean((OR_inf[!is.na(OR_inf) & !is.nan(OR_inf) & !is.infinite(OR_inf)] - or)^2)
MSE_inr <- mean((OR_inr[!is.na(OR_inr) & !is.nan(OR_inr) & !is.infinite(OR_inr)] - or)^2)
MSE_Ipd <- mean((OR_Ipd[!is.na(OR_Ipd) & !is.nan(OR_Ipd) & !is.infinite(OR_Ipd)] - or)^2)


OR_Ipd_m 
OR_peto_m
OR_mhf_m
OR_mhr_m
OR_inf_m
OR_inr_m

OR_Ipd_percentb 
OR_peto_percentb 
OR_mhf_percentb 
OR_mhr_percentb 
OR_inf_percentb 
OR_inr_percentb 

MSE_Ipd
MSE_peto 
MSE_mhf 
MSE_mhr 
MSE_inf
MSE_inr

cov_IPD  
cov_PETO   
cov_MF
cov_MR
cov_IF
cov_IR