inst/doc/required-sample-size.R

In RTSA: 'Trial Sequential Analysis' for Error Control and Inference in Sequential Meta-Analyses

## ---- include = FALSE---------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
options(rmarkdown.html_vignette.check_title = FALSE)

## ----setup--------------------------------------------------------------------
library(RTSA)

## -----------------------------------------------------------------------------
ris(outcome = "RR", mc = 0.9, pC = 0.1, alpha = 0.05, beta = 0.1, side = 2)

## -----------------------------------------------------------------------------
ris(outcome = "RR", mc = 0.9, pC = 0.2, fixed = FALSE, I2 = 0.2, D2 = 0.3,
    side = 2, alpha = 0.05, beta = 0.2)

## ---- eval = FALSE, include=FALSE---------------------------------------------
#  # simulate 10 trials
#  tau2 = 0.05
#  outl = matrix(NA, ncol = 3, nrow = 10)
#  outm = matrix(NA, ncol = 4, nrow = 10)
#  
#  RR <- 0.9
#  pC <- 0.1
#  pI <- round(exp(log(pC) + log(RR) / 2), 4)
#  pC <- round(exp(log(pC) - log(RR) / 2), 4)
#  theta = round(pC - pI, 4)
#  n = 2500
#  K = 10
#  nsim = 1000
#  
#  for(m in 1:10) {
#    outpvalue = matrix(NA, ncol = 3, nrow = nsim)
#    outhetero = matrix(NA, ncol = 5, nrow = nsim)
#  
#    for (h in 1:nsim) {
#      ln_RR = rnorm(K, mean = log(RR), sd = sqrt(tau2))
#  
#      pI = exp(log(pC) + ln_RR)
#      pI[pI < 0.01] = 0.01
#      pI[pI > 0.99] = 0.99
#      pC = rep(exp(log(pC)),K)
#      pC[pC < 0.01] = 0.01
#      pC[pC > 0.99] = 0.99
#      outmat = matrix(NA, ncol = 4, nrow = K)
#      zvalues = NULL
#      for (i in 1:K) {
#        eA <- apply(cbind(n / 2, pI[i]), 1,
#                    function(x)
#                      rbinom(1, size = x[1], prob = x[2]))
#        eB <- apply(cbind(n / 2, pC[i]), 1,
#                    function(x)
#                      rbinom(1, size = x[1], prob = x[2]))
#        outmat[i, 1:4] = c(eA, n / 2, eB, n / 2)
#      }
#  
#      dat <- data.frame(eI = outmat[, 1],
#        nI = outmat[, 2],
#        eC = outmat[, 3],
#        nC = outmat[, 2])
#  
#      synout = RTSA:::metaPrepare(
#        outcome = "RR", data = dat,
#        weights = "MH", alpha = 0.05
#      )
#      out1 = RTSA:::synthesize(synout, tau_ci_method = "BJ",
#                        re_method = "DL")
#  
#      ma <- metaanalysis(outcome = "RR", data = dat, mc = 0.9)
#  
#          # save the tau^2, I^2 and D^2
#      hetero = c(out1$U[c(1, 3, 4)])
#      #RISd2 = 1 / (1 - hetero[3]) * ma$ris$full_NF
#      #outhetero[h,] = c(hetero, ma$ris$NR_inc_full, ma$ris$NR_div_full)
#      dRISd2 <- ma$ris$NR_div
#      if(is.null(dRISd2)){
#        dRISd2 <- ma$ris$NF
#      }
#  
#      ln_RR = rnorm(m, mean = log(RR), sd = sqrt(tau2))
#      pI = exp(log(pC) + ln_RR)
#      pI[pI < 0.01] = 0.01
#      pI[pI > 0.99] = 0.99
#      pC = rep(exp(log(pC)),m)
#      pC[pC < 0.01] = 0.01
#      pC[pC > 0.99] = 0.99
#  
#      for (b in 1:m) {
#        eA <- apply(cbind(ceiling(dRISd2 / 2 / m), pI[b]), 1,
#                    function(x)
#                      rbinom(1, size = x[1], prob = x[2]))
#        eB <- apply(cbind(ceiling(dRISd2 / 2 / m), pC[b]), 1,
#                    function(x)
#                      rbinom(1, size = x[1], prob = x[2]))
#        outmat = rbind(outmat, c(eA, dRISd2 / 2 / m, eB, dRISd2 / 2 / m))
#      }
#  
#      dat <- data.frame(eI = outmat[, 1],
#        nI = outmat[, 2],
#        eC = outmat[, 3],
#        nC = outmat[, 2])
#  
#      synout = RTSA:::metaPrepare(
#        outcome = "RR",
#        weights = "MH", data = dat, alpha = 0.05
#      )
#  
#      out1 = RTSA:::synthesize(synout, re_method = "DL", tau_ci_method = "BJ")
#      out2 = RTSA:::synthesize(synout, re_method = "DL_HKSJ", tau_ci_method = "BJ")
#  
#      if(is.null(out1$peR[5])) out1$peR[5] <- out1$peF[5]
#      if(is.null(out2$peR[5])) out2$peR[5] <- out1$peF[5]
#  
#      outpvalue[h, c(1, 2, 3)] = c(round(out1$peF[5], 4), round(out1$peR[5], 4), out2$peR[5])
#    }
#    outm[m,] = c(
#      sum(outpvalue[, 1] < 0.05) / nsim,
#      sum(outpvalue[, 2] < 0.05) / nsim,
#      sum(outpvalue[, 3] < 0.05) / nsim,
#      mean(dRISd2)
#    )
#  }
#  
#  outm = cbind(1:10, outm)
#  colnames(outm) = c(
#    "Number of extra trials",
#    "Fixed-effect",
#    "Random-effects DL",
#    "Random-effects HKSJ",
#    "Avg. RIS"
#  )
#  rownames(outm) = rep(c(""), 10)
#  save(outm, file = "random-effects-TSA.Rda")

## ----randomTSA, echo = FALSE--------------------------------------------------
load("random-effects-TSA.Rda")
knitr::kable(outm[,-5], caption = "Power per model as a function of number of extra trials and RIS based on Diversity")

## -----------------------------------------------------------------------------
ris(outcome = "RR", mc = 0.9, fixed = FALSE, tau2 = 0.05, pC = 0.1,
    side = 2, alpha = 0.05, beta = 0.2)

## ---- eval = FALSE, echo = FALSE----------------------------------------------
#  outl = matrix(NA, ncol = 3, nrow = 4)
#  
#  RR <- 0.9
#  pC <- 0.1
#  pI <- round(exp(log(pC)+log(RR)/2),4)
#  pC <- round(exp(log(pC)-log(RR)/2),4)
#  theta = round(pC - pI,4)
#  nsim = 10000
#  tau2 = 0.05
#  
#  trial.out = minTrial(outcome = "RR", mc = 0.9, tau2 = 0.05, pC = 0.1)
#  
#  outtau = numeric(nsim)
#  outpvalue = matrix(NA, ncol = 3, nrow = nsim)
#  
#  for(l in 1:dim(trial.out$nPax)[2]){
#  K = trial.out$nPax[1,l]
#  n = trial.out$nPax[2,l]
#  
#    for(h in 1:nsim){
#      outmat = matrix(NA,ncol = 4, nrow = K)
#  
#      ln_RR = rnorm(K, mean = log(RR), sd = sqrt(tau2))
#      pI = exp(log(pC)+ln_RR/2)
#      pI[pI < 0.01] = 0.01
#      pI[pI > 0.99] = 0.99
#      pC = exp(log(pC)-ln_RR/2)
#      pC[pC < 0.01] = 0.01
#      pC[pC > 0.99] = 0.99
#      for(i in 1:(K)){
#        eA <- apply(cbind(ceiling(n/2), pI[i]), 1,
#                    function(x) rbinom(1, size = x[1], prob = x[2]))
#        eB <- apply(cbind(ceiling(n/2), pC[i]), 1,
#                    function(x) rbinom(1, size = x[1], prob = x[2]))
#        outmat[i,1:4] = c(eA, ceiling(n/2), eB, ceiling(n/2))
#      }
#  
#      dat <- data.frame(eI = outmat[,1], nI = outmat[,2],
#                            eC = outmat[,3], nC = outmat[,2])
#  
#      synout = RTSA:::metaPrepare(outcome = "RR", data = dat,
#                            weights = "IV", alpha = 0.05)
#      out1 = RTSA:::synthesize(synout, re_method = "DL", tau_ci_method = "BJ")
#      out2 = RTSA:::synthesize(synout, re_method = "DL_HKSJ", tau_ci_method = "BJ")
#      outpvalue[h, c(1,2,3)] = c(round(out1$peF[5],4), round(out1$peR[5],4), out2$peR[5])
#    }
#    outl[l,1] = sum(outpvalue[,1] <= 0.05)/nsim
#    outl[l,2] = sum(outpvalue[,2] <= 0.05)/nsim
#    outl[l,3] = sum(outpvalue[,3] <= 0.05)/nsim
#  }
#  
#  outl = cbind(outl, trial.out$nPax[1,], trial.out$nPax[2,])
#  colnames(outl) = c("Fixed-effect", "Random-effects DL", "Random-effects HKSJ",
#                     "Number of trials", "Participants per trial")
#  rownames(outl) = c("","","","")
#  save(outl, file = "random-effects.Rda")

## ----random, echo = FALSE-----------------------------------------------------
load("random-effects.Rda")
knitr::kable(outl[,c(4,5,1,2,3)], caption = "Power per model as a function of number of trials and number of participants per trial")

## -----------------------------------------------------------------------------
ma <- metaanalysis(data = perioOxy, outcome = "RR", mc = 0.8, beta = 0.1,
                   pC = 0.15)
ma$ris

## -----------------------------------------------------------------------------
ma$ris$NR_tau$NR_tau_ll

## -----------------------------------------------------------------------------
ma$ris$NR_tau$NR_tau_ul