inst/doc/contingency-tables-vignette.R

In safestats: Safe Anytime-Valid Inference

## ----setup, include = FALSE---------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.height = 4,
  fig.width = 8
)

library(safestats)

## -----------------------------------------------------------------------------
ya <- c(1,1,1,1,1,1,0,1,1,1)
yb <- c(0,0,1,0,1,0,0,0,0,0)

## ---- echo = FALSE------------------------------------------------------------
set.seed(1082021)
# balancedSafeDesign <- designSafeTwoProportions(
#   na = 1,
#   nb = 1,
#   alpha = 0.05, #significance level for testing
#   beta = 1 - 0.8, #power 80% equals type II error (beta) of 0.2
#   delta = 0.3
# )
# save(balancedSafeDesign, file="safe2x2VignetteData/safe2x2DesignObject.Rdata")
load("safe2x2VignetteData/safe2x2DesignObject.Rdata")

## ---- eval=FALSE--------------------------------------------------------------
#  balancedSafeDesign <- designSafeTwoProportions(
#    na = 1,
#    nb = 1,
#    alpha = 0.05, #significance level for testing
#    beta = 1 - 0.8, #power 80% equals type II error (beta) of 0.2
#    delta = 0.3
#  )

## -----------------------------------------------------------------------------
print(balancedSafeDesign)

## -----------------------------------------------------------------------------
set.seed(19012022)
successProbabilityA <- 0.2
successProbabilityB <- 0.5
ya <- rbinom(n = balancedSafeDesign[["nPlan"]]["nBlocksPlan"], size = 1, prob = successProbabilityA)
yb <- rbinom(n = balancedSafeDesign[["nPlan"]]["nBlocksPlan"], size = 1, prob = successProbabilityB)

graphics::plot(x = seq_along(ya), y = cumsum(ya), type = "p", ylim = c(0, max(c(sum(ya), sum(yb)))),
     xlab = "block number", ylab = "# successes",
     main = "Total # successes per group over time")
graphics::points(x = seq_along(yb), y = cumsum(yb), type = "p", col = "grey")
graphics::legend(x = "topleft", legend = c("A", "B"), col = c("black", "grey"), pch = 1)

## -----------------------------------------------------------------------------
safe.prop.test(ya = ya, yb = yb, designObj = balancedSafeDesign)

## ---- echo = FALSE------------------------------------------------------------
# #make the plot voor standard parameter settings
# standardPrior <- list(betaA1 = 0.18,
#                       betaA2 = 0.18,
#                       betaB1 = 0.18,
#                       betaB2 = 0.18)
# #set a seed for the simulations
# set.seed(1082021)
# sampleSizePerMinimalDifference <- simulateTwoProportions(list(standard = standardPrior),
#                                            alternativeRestriction = "none",
#                                            alpha = 0.05,
#                                            beta = 0.2,
#                                            na = 1,
#                                            nb = 1,
#                                            deltamin = 0.3,
#                                            deltamax = 0.9
#                                          )
# save(sampleSizePerMinimalDifference, file = "safe2x2VignetteData/safe2x2SimObject.Rdata")
load("safe2x2VignetteData/safe2x2SimObject.Rdata")

## ---- eval = FALSE------------------------------------------------------------
#  #make the plot voor standard parameter settings
#  standardPrior <- list(betaA1 = 0.18,
#                        betaA2 = 0.18,
#                        betaB1 = 0.18,
#                        betaB2 = 0.18)
#  #set a seed for the simulations
#  set.seed(1082021)
#  sampleSizePerMinimalDifference <- simulateTwoProportions(list(standard = standardPrior),
#                                             alternativeRestriction = "none",
#                                             alpha = 0.05,
#                                             beta = 0.2,
#                                             na = 1,
#                                             nb = 1,
#                                             deltamin = 0.3,
#                                             deltamax = 0.9
#                                           )

## ---- echo = FALSE------------------------------------------------------------
someIndex <- max(which(sampleSizePerMinimalDifference$simdata$worstCaseQuantile <= 50))
minDetectableEffect <- round(sampleSizePerMinimalDifference$simdata$delta[someIndex], 2)

## -----------------------------------------------------------------------------
graphics::plot(sampleSizePerMinimalDifference)

## ---- echo = FALSE------------------------------------------------------------
# unbalancedSafeDesign <- designSafeTwoProportions(
#   na = 2,
#   nb = 1,
#   alpha = 0.05,
#   beta = 0.2,
#   delta = 0.3
# )
# save(unbalancedSafeDesign, file = "safe2x2VignetteData/safe2x2UnbalancedDesignObject.Rdata")
load("safe2x2VignetteData/safe2x2UnbalancedDesignObject.Rdata")

## ---- eval = FALSE------------------------------------------------------------
#  unbalancedSafeDesign <- designSafeTwoProportions(
#    na = 2,
#    nb = 1,
#    alpha = 0.05,
#    beta = 0.2,
#    delta = 0.3
#  )

## -----------------------------------------------------------------------------
set.seed(692021)
yaUnbalanced <- rbinom(n = unbalancedSafeDesign[["nPlan"]]["nBlocksPlan"], size = 2, prob = successProbabilityA)
ybUnbalanced <- rbinom(n = unbalancedSafeDesign[["nPlan"]]["nBlocksPlan"], size = 1, prob = successProbabilityB)

print(yaUnbalanced)
print(ybUnbalanced)

## -----------------------------------------------------------------------------
safeTwoProportionsTest(ya = yaUnbalanced, yb = ybUnbalanced, designObj = unbalancedSafeDesign)

## -----------------------------------------------------------------------------
print(unbalancedSafeDesign)

## ---- echo = FALSE------------------------------------------------------------
# #make the plot voor standard parameter settings
# standardPrior <- unbalancedSafeDesign[["betaPriorParameterValues"]]
# uniformLikePrior <- list(betaA1 = 2, # pretend we "have seen" 2 success and
#                       betaA2 = 2,# 2 failure in group A before starting
#                       betaB1 = 1, # and 1 success and failure in group B
#                       betaB2 = 1)
# #set a seed for the simulations
# set.seed(1082021)
# sampleSizePerMinimalDifference <- simulateTwoProportions(
#    hyperparameterList = list(standard = standardPrior, uniform = uniformLikePrior),
#    alternativeRestriction = "none",
#    alpha = 0.05,
#    beta = 0.2,
#    na = 2,
#    nb = 1,
#    deltamin = 0.3,
#    deltamax = 0.9
#   )
# save(sampleSizePerMinimalDifference, file = "safe2x2VignetteData/safe2x2UnbalancedSimObject.Rdata")
load("safe2x2VignetteData/safe2x2UnbalancedSimObject.Rdata")

## ---- eval = FALSE------------------------------------------------------------
#  #make the plot voor standard parameter settings
#  standardPrior <- unbalancedSafeDesign[["betaPriorParameterValues"]]
#  uniformLikePrior <- list(betaA1 = 2, # pretend we "have seen" 2 success and
#                        betaA2 = 2,# 2 failure in group A before starting
#                        betaB1 = 1, # and 1 success and failure in group B
#                        betaB2 = 1)
#  #set a seed for the simulations
#  set.seed(1082021)
#  sampleSizePerMinimalDifference <- simulateTwoProportions(
#     hyperparameterList = list(standard = standardPrior, uniform = uniformLikePrior),
#     alternativeRestriction = "none",
#     alpha = 0.05,
#     beta = 0.2,
#     na = 2,
#     nb = 1,
#     deltamin = 0.3,
#     deltamax = 0.9
#    )

## -----------------------------------------------------------------------------
graphics::plot(sampleSizePerMinimalDifference)

## ---- echo = FALSE------------------------------------------------------------
# differenceRestrictedSafeDesign <- designSafeTwoProportions(
#   na = 1,
#   nb = 1,
#   alpha = 0.05, #significance level for testing
#   beta = 1 - 0.8, #power 80% equals type II error (beta) of 0.2
#   delta = 0.3,
#   alternativeRestriction = "difference" #also available: logOddsRatio
# )
# save(differenceRestrictedSafeDesign, file = "safe2x2VignetteData/safe2x2DiffRestrictedDesign.Rdata")
load("safe2x2VignetteData/safe2x2DiffRestrictedDesign.Rdata")

## ---- eval = FALSE------------------------------------------------------------
#  differenceRestrictedSafeDesign <- designSafeTwoProportions(
#    na = 1,
#    nb = 1,
#    alpha = 0.05, #significance level for testing
#    beta = 1 - 0.8, #power 80% equals type II error (beta) of 0.2
#    delta = 0.3,
#    alternativeRestriction = "difference" #also available: logOddsRatio
#  )

## -----------------------------------------------------------------------------
safe.prop.test(ya = ya, yb = yb, designObj = differenceRestrictedSafeDesign)

## -----------------------------------------------------------------------------
#set a seed for the simulations
set.seed(1082021)
optionalStoppingResult <- simulateOptionalStoppingScenarioTwoProportions(
    safeDesign = balancedSafeDesign,
    M = 1e3,
    thetaA = 0.2,
    thetaB = 0.5
)

## -----------------------------------------------------------------------------
plotHistogramDistributionStoppingTimes(safeSim = optionalStoppingResult, 
                                       nPlan = balancedSafeDesign[["nPlan"]]["nBlocksPlan"],
                                       deltaTrue = 0.3)

## -----------------------------------------------------------------------------
optionalStoppingResult[["powerOptioStop"]]

## -----------------------------------------------------------------------------
#set a seed for the simulations
set.seed(1082021)
optionalStoppingResultTrueDifferenceBigger <-
  simulateOptionalStoppingScenarioTwoProportions(
    safeDesign = balancedSafeDesign,
    M = 1e3,
    thetaA = 0.2,
    thetaB = 0.7
  )

## -----------------------------------------------------------------------------
plotHistogramDistributionStoppingTimes(safeSim = optionalStoppingResultTrueDifferenceBigger, 
                                       nPlan = balancedSafeDesign[["nPlan"]]["nBlocksPlan"],
                                       deltaTrue = 0.5)

## -----------------------------------------------------------------------------
optionalStoppingResultTrueDifferenceBigger[["powerOptioStop"]]

## -----------------------------------------------------------------------------
#set a seed for the simulations
set.seed(1082021)
optionalStoppingResultTrueDifferenceNull <- simulateOptionalStoppingScenarioTwoProportions(safeDesign = balancedSafeDesign, 
                                                    M = 1e3, 
                                                    thetaA = 0.5, 
                                                    thetaB = 0.5)

## -----------------------------------------------------------------------------
plotHistogramDistributionStoppingTimes(safeSim = optionalStoppingResultTrueDifferenceNull, 
                                       nPlan = balancedSafeDesign[["nPlan"]]["nBlocksPlan"],
                                       deltaTrue = 0)

## -----------------------------------------------------------------------------
optionalStoppingResultTrueDifferenceNull[["powerOptioStop"]]

## ---- echo = FALSE------------------------------------------------------------
# optionalStoppingWrongFisher <- simulateIncorrectStoppingTimesFisher(thetaA = 0.5,
#                                      thetaB = 0.5,
#                                      alpha = 0.05,
#                                      na = balancedSafeDesign[["nPlan"]][["na"]],
#                                      nb = balancedSafeDesign[["nPlan"]][["nb"]],
#                                      maxSimStoptime = balancedSafeDesign[["nPlan"]][["nBlocksPlan"]],
#                                      M = 1e3,
#                                      numberForSeed = 1082021)
# save(optionalStoppingWrongFisher, file = "safe2x2VignetteData/safe2x2SimFisher.Rdata")
load("safe2x2VignetteData/safe2x2SimFisher.Rdata")

## ---- eval = FALSE------------------------------------------------------------
#  optionalStoppingWrongFisher <-
#    simulateIncorrectStoppingTimesFisher(
#      thetaA = 0.5,
#      thetaB = 0.5,
#      alpha = 0.05,
#      na = balancedSafeDesign[["nPlan"]][["na"]],
#      nb = balancedSafeDesign[["nPlan"]][["nb"]],
#      maxSimStoptime = balancedSafeDesign[["nPlan"]][["nBlocksPlan"]],
#      M = 1e3,
#      numberForSeed = 1082021
#    )

## -----------------------------------------------------------------------------
mean(optionalStoppingWrongFisher[["rejections"]] == 1)

## -----------------------------------------------------------------------------
eValuesNotRejected <-
  optionalStoppingResult[["eValues"]][!optionalStoppingResult[["allSafeDecisions"]]]
eValuesNotRejectedNull <-
  optionalStoppingResultTrueDifferenceNull[["eValues"]][!optionalStoppingResultTrueDifferenceNull[["allSafeDecisions"]]]

## ---- echo = FALSE------------------------------------------------------------
trueHist <- graphics::hist(x = eValuesNotRejected, plot = FALSE)
nullHist <- graphics::hist(x = eValuesNotRejectedNull, plot = FALSE)
oldPar <- graphics::par(cex.main=1.5, mar=c(5, 6, 4, 4)+0.1, mgp=c(3.5, 1, 0), cex.lab=1.5,
              font.lab=2, cex.axis=1.3, bty="n", las=1)
graphics::plot(nullHist, xlim = c(0, max(eValuesNotRejected, eValuesNotRejectedNull)),
               freq = FALSE, col = "blue", density = 20, angle = 45, xlab = "e-values",
               main = "Histogram of e-values where null not rejected")
graphics::plot(trueHist, add = TRUE, freq = FALSE, col = "red", density = 20,
               angle = -45)
graphics::legend(x = "topright", legend = c("True delta: null", "True delta: design"), fill = c("blue", "red"))
par(oldPar)

## -----------------------------------------------------------------------------
promisingEValues <- eValuesNotRejected[eValuesNotRejected >= 10]
followUpDesign <- designSafeTwoProportions(na = 1, nb = 1, nBlocksPlan = 30)
followUpResults <- simulateOptionalStoppingScenarioTwoProportions(followUpDesign, 
                                                                  M = length(promisingEValues), 
                                                                  thetaA = 0.2, thetaB = 0.5)
newEValues <- followUpResults[["eValues"]] * promisingEValues

## ---- echo = FALSE, cache = FALSE---------------------------------------------
trueHist <- graphics::hist(x = newEValues, plot = FALSE)
oldPar <- graphics::par(cex.main=1.5, mar=c(5, 6, 4, 4)+0.1, mgp=c(3.5, 1, 0), cex.lab=1.5,
                        font.lab=2, cex.axis=1.3, bty="n", las=1)
graphics::plot(trueHist, xlim = c(0, max(newEValues)),
               freq = FALSE, col = "blue", density = 20, angle = 45, xlab = "e-values",
               main = "Histogram of e-values after continuing data collection")
par(oldPar)

## ---- cache = FALSE-----------------------------------------------------------
sum(newEValues >= 20)/length(promisingEValues)

## -----------------------------------------------------------------------------
promisingEValuesNull <-
  eValuesNotRejectedNull[eValuesNotRejectedNull >= 1]
followUpDesign <-
  designSafeTwoProportions(na = 1,
                           nb = 1,
                           nBlocksPlan = 30)
followUpResults <-
  simulateOptionalStoppingScenarioTwoProportions(
    followUpDesign,
    M = length(promisingEValuesNull),
    thetaA = 0.5,
    thetaB = 0.5
  )
newEValuesNull <-
  followUpResults[["eValues"]] * promisingEValuesNull

## ---- echo = FALSE, cache = FALSE---------------------------------------------
nullHist <- graphics::hist(x = newEValuesNull, plot = FALSE)
oldPar <- graphics::par(cex.main=1.5, mar=c(5, 6, 4, 4)+0.1, mgp=c(3.5, 1, 0), cex.lab=1.5,
                        font.lab=2, cex.axis=1.3, bty="n", las=1)
graphics::plot(nullHist, xlim = c(0, max(newEValuesNull)),
               freq = FALSE, col = "blue", density = 20, angle = 45, xlab = "e-values",
               main = "Histogram of e-values after continuing data collection under the null")
par(oldPar)

## ---- eval = FALSE------------------------------------------------------------
#  plotConfidenceSequenceTwoProportions(ya = ya,
#                                       yb = yb,
#                                       safeDesign = balancedSafeDesign,
#                                       differenceMeasure = "difference",
#                                       precision = 100,
#                                       trueDifference = 0.3)

## ---- echo = FALSE------------------------------------------------------------
# coverageSimResult <- simulateCoverageDifferenceTwoProportions(successProbabilityA = 0.2,
#                                            trueDelta = 0.3,
#                                            safeDesign = balancedSafeDesign,
#                                            numberForSeed = 1082021)
# save(coverageSimResult, file = "safe2x2VignetteData/safe2x2Coverage.Rdata")
load("safe2x2VignetteData/safe2x2Coverage.Rdata")

## ---- eval = FALSE------------------------------------------------------------
#  coverageSimResult <- simulateCoverageDifferenceTwoProportions(successProbabilityA = 0.2,
#                                             trueDelta = 0.3,
#                                             safeDesign = balancedSafeDesign,
#                                             numberForSeed = 1082021)

## -----------------------------------------------------------------------------
print(coverageSimResult)

## -----------------------------------------------------------------------------
log(0.5/0.5 * 0.8/0.2)

## ---- echo = FALSE------------------------------------------------------------
# confidenceBoundOdds <- computeConfidenceBoundForLogOddsTwoProportions(ya = ya,
#                                            yb = yb,
#                                            safeDesign = balancedSafeDesign,
#                                            bound = "lower",
#                                            precision = 100,
#                                            deltaStart = 0.001,
#                                            deltaStop = 2)
# save(confidenceBoundOdds, file = "safe2x2VignetteData/safe2x2ConfidenceOdds.Rdata")
load("safe2x2VignetteData/safe2x2ConfidenceOdds.Rdata")

## ---- eval = FALSE------------------------------------------------------------
#  confidenceBoundOdds <- computeConfidenceBoundForLogOddsTwoProportions(ya = ya,
#                                             yb = yb,
#                                             safeDesign = balancedSafeDesign,
#                                             bound = "lower",
#                                             precision = 100,
#                                             deltaStart = 0.001,
#                                             deltaStop = 2)

## -----------------------------------------------------------------------------
print(confidenceBoundOdds)

## -----------------------------------------------------------------------------
computeConfidenceBoundForLogOddsTwoProportions(ya = ya,
                                           yb = yb, 
                                           safeDesign = balancedSafeDesign,
                                           bound = "upper", 
                                           precision = 10, 
                                           deltaStart = -0.01, 
                                           deltaStop = -2)

## ---- eval = FALSE------------------------------------------------------------
#  plotConfidenceSequenceTwoProportions(ya = ya,
#                                       yb = yb,
#                                       safeDesign = balancedSafeDesign,
#                                       differenceMeasure = "odds",
#                                       precision = 100,
#                                       trueDifference = log(0.5/0.5 * 0.8/0.2),
#                                       deltaStart = 0.001,
#                                       deltaStop = 2)