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

## ----eval=FALSE---------------------------------------------------------------
# library("safestats")
# freqColours <- c("#E31A1CE6", "#FB9A9980")
# bayesColours <- c("#B15928E6", "#FFFF9980")
# eColours <- c("#1F78B4E6", "#A6CEE380")

## ----eval=FALSE---------------------------------------------------------------
# n1 <- n2 <- 100
# sigma <- 12
# alpha <- 0.05
# 
# muGlobal <- 112
# 
# set.seed(2)
# dataGroup1 <- rnorm(n1, mean=muGlobal, sd=sigma)
# set.seed(3)
# dataGroup2 <- rnorm(n2, mean=muGlobal, sd=sigma)

## ----eval=FALSE---------------------------------------------------------------
# # # PSEUDO CODE
# # #
# # #     THIS CODE DOES NOT RUN
# # #
# # # Initiate
# # n <- 1
# # pValue <- 1
# #
# # while (pValue > alpha) {
# #   pValue <- computePValueZTest(x=x[1:n], y=y[1:n])
# #
# #   if (pValue < alpha) {
# #     "Reject the null"
# #     stop()
# #   } else {
# #     "Increase sample size and test again
# #           at the start the start of the while loop"
# #     n <- n + 1
# #   }
# # }

## ----eval=FALSE---------------------------------------------------------------
# pValueVector <- vector("numeric", length=n1)
# 
# for (i in 1:n1) {
#   pValueVector[i] <- pValueZTest(x=dataGroup1[1:i],
#                                  y=dataGroup2[1:i],
#                                  sigma=sigma)$pValue
# }
# 
# oldPar <- setSafeStatsPlotOptionsAndReturnOldOnes();
# 
# plot(1:n1, pValueVector, type="l", lwd=2, xlab="n",
#      ylab="p-value",
#      col=freqColours[1])
# abline(h=0.05, lty=2, lwd=2)

## ----eval=FALSE---------------------------------------------------------------
# mIter <- 1000
# 
# allData <- generateNormalData(c(n1, n2), muGlobal=muGlobal,
#                               nSim=mIter,
#                               meanDiffTrue=0, seed=1,
#                               sigma=sigma)
# allPValues <- matrix(nrow=mIter, ncol=n1)
# 
# # This indicates whether a simulation yielded a "significant" p-value
# pValueUnderAlpha <- vector("integer", length=mIter)
# 
# # This indicates the first time an experiment yielded a "significant" p-value
# # Default is Inf, which indicates that the p-value didn't dip below alpha
# firstPassageTime <- rep(Inf, times=mIter)
# 
# # Used to vectorise the computations for the the z-statistic
# n1Vector <- 1:n1
# n2Vector <- 1:n2
# nEffVector <- (1/n1Vector+1/n2Vector)^(-1)
# 
# 
# for (sim in 1:mIter) {
#   dataGroup1 <- allData$dataGroup1[sim, ]
#   dataGroup2 <- allData$dataGroup2[sim, ]
# 
#   x1BarVector <- 1/n1Vector*cumsum(dataGroup1)
#   x2BarVector <- 1/n2Vector*cumsum(dataGroup2)
#   zVector <- sqrt(nEffVector)*(x1BarVector - x2BarVector)/sigma
# 
#   for (i in 1:n1) {
#     currentPValue <- pValueFromZStat(zVector[i])
#     allPValues[sim, i] <- currentPValue
# 
#     if (currentPValue < alpha && pValueUnderAlpha[sim]!=1) {
#       pValueUnderAlpha[sim] <- 1
#       firstPassageTime[sim] <- i
#     }
#   }
# }

## ----eval=FALSE---------------------------------------------------------------
# numberOfDippingExperimentsAtTimeN <- integer(n1)
# 
# for (i in 1:n1) {
#   numberOfDippingExperimentsAtTimeN[i] <- sum(firstPassageTime <= i)
# }
# 
# pValueFalseRejects <-numberOfDippingExperimentsAtTimeN/mIter
# 
# oldPar <- setSafeStatsPlotOptionsAndReturnOldOnes();
# plot(1:n1, 100*pValueFalseRejects, type="l", xlab="n",
#      ylab="False positive rate (%)",
#      lwd=2, col=freqColours[1])
# lines(c(1, n1), c(5, 5), lwd=2, lty=2)

## ----eval=FALSE---------------------------------------------------------------
# # Null hyper parameters
# a0 <- 107
# b0 <- 2
# 
# # Alternative model
# #  Hyper parameters in the control group
# #
# a2 <- 107
# b2 <- 2
# 
# #  Hyper parameters in the treatment group
# #
# a1 <- 117
# b1 <- 4

## ----eval=FALSE---------------------------------------------------------------
# # # PSEUDO CODE
# # #
# # #     THIS CODE DOES NOT RUN
# # #
# # # Initiate
# # n <- 1
# # bf10 <- 1
# #
# # while (bf10 < 30) {
# #   bf10 <- computeBfZTest(x=x[1:n], y=y[1:n])
# #
# #   if (bf10 > 30) {
# #     "Reject the null"
# #     stop()
# #   } else {
# #     "Increase sample size and test again
# #           at the start of the while loop"
# #     n <- n + 1
# #   }
# # }

## ----eval=FALSE---------------------------------------------------------------
# # Single sample path data set as in the p-value example
# set.seed(2)
# dataGroup1 <- rnorm(n1, mean=muGlobal, sd=sigma)
# set.seed(3)
# dataGroup2 <- rnorm(n2, mean=muGlobal, sd=sigma)
# 
# bfVector <- vector("numeric", length=n1)
# 
# for (i in 1:n1) {
#   x1 <- mean(dataGroup1[1:i])
#   s21 <- var(dataGroup1[1:i])
# 
#   x2 <- mean(dataGroup2[1:i])
#   s22 <- var(dataGroup2[1:i])
# 
#   bfVector[i] <- subjectiveBfZStat(x1=x1, s21=s21, n1=i,
#                                    x2=x2, s22=s22, n2=i,
#                                    sigma=sigma,
#                                    a1=a1, b1=b1,
#                                    a2=a2, b2=b2,
#                                    a0=a0, b0=b0)
# }
# 
# oldPar <- setSafeStatsPlotOptionsAndReturnOldOnes();
# 
# plot(1:n1, bfVector, type="l", lwd=2, xlab="n", ylab="BF10",
#      col=bayesColours[1])
# lines(c(1, n1), c(30, 30), lwd=2, lty=3)
# lines(c(1, n1), c(1/alpha, 1/alpha), lwd=2, lty=2)

## ----eval=FALSE---------------------------------------------------------------
# mIter <- 1000
# 
# allData <- generateNormalData(c(n1, n2), muGlobal=muGlobal,
#                               nSim=mIter,
#                               meanDiffTrue=0, seed=1,
#                               sigma=sigma)
# allBfValues <- matrix(nrow=mIter, ncol=n1)
# 
# # This indicates whether a simulation yielded a "significant" Bayes factor values
# bfValueOver20 <- bfValueOver30 <- vector("integer", length=mIter)
# 
# # This indicates the first time an experiment yielded a "significant" Bayes factor
# # Default is Inf, which indicates that the bf didn't cross 1/alpha=20, or 30 respectively
# firstPassageTime20 <- firstPassageTime30 <- rep(Inf, times=mIter)
# 
# for (sim in 1:mIter) {
#   dataGroup1 <- allData$dataGroup1[sim, ]
#   dataGroup2 <- allData$dataGroup2[sim, ]
# 
#   for (i in 1:n1) {
#     x1 <- mean(dataGroup1[1:i])
#     s21 <- var(dataGroup1[1:i])
# 
#     x2 <- mean(dataGroup2[1:i])
#     s22 <- var(dataGroup2[1:i])
# 
#     currentBf10 <- subjectiveBfZStat(x1=x1, s21=s21, n1=i,
#                                      x2=x2, s22=s22, n2=i,
#                                      sigma=sigma,
#                                      a1=a1, b1=b1,
#                                      a2=a2, b2=b2,
#                                      a0=a0, b0=b0)
# 
#     allBfValues[sim, i] <- currentBf10
# 
#     if (currentBf10 > 30 && bfValueOver30[sim]!=1) {
#       bfValueOver30[sim] <- 1
#       firstPassageTime30[sim] <- i
#     }
# 
#     # Prepare for the comparison with the e-value procedure
#     if (currentBf10 > 1/alpha && bfValueOver20[sim]!=1) {
#       bfValueOver20[sim] <- 1
#       firstPassageTime20[sim] <- i
#     }
#   }
# }

## ----eval=FALSE---------------------------------------------------------------
# numberOfCrossingExperimentsAtTimeN20 <- numberOfCrossingExperimentsAtTimeN30 <- integer(n1)
# 
# for (i in 1:n1) {
#   numberOfCrossingExperimentsAtTimeN30[i] <- sum(firstPassageTime30 <= i)
#     numberOfCrossingExperimentsAtTimeN20[i] <- sum(firstPassageTime20 <= i)
# }
# 
# subjectiveBfValueFalseRejects30 <-numberOfCrossingExperimentsAtTimeN30/mIter
# 
# subjectiveBfValueFalseRejects <-numberOfCrossingExperimentsAtTimeN20/mIter
# 
# oldPar <- setSafeStatsPlotOptionsAndReturnOldOnes();
# plot(1:n1, 100*subjectiveBfValueFalseRejects, type="l",
#      xlab="n", ylab="False positive rate (%)",
#      lwd=2, col=bayesColours[1], ylim=c(0, 20))
# 
# lines(1:n1, 100*subjectiveBfValueFalseRejects30, lwd=2,
#       col=bayesColours[1], lty=2)
# lines(c(1, n1), c(5, 5), lwd=2, lty=2)

## ----eval=FALSE---------------------------------------------------------------
# failedSimIndeces <- which(bfValueOver20==1)
# 
# # Take an arbitrary index where it fails
# someIndex <- failedSimIndeces[3]
# 
# oldPar <- setSafeStatsPlotOptionsAndReturnOldOnes();
# plot(1:n1, allBfValues[someIndex, ], type="l", lwd=2,
#      xlab="n", ylab="BF10", col=bayesColours[1])
# lines(c(1, n1), c(1/alpha, 1/alpha), lwd=2, lty=2)

## ----eval=FALSE---------------------------------------------------------------
# designObj <- designSaviZ(meanDiffMin=10,
#                          testType="twoSample",
#                          sigma=sigma)
# designObj

## ----eval=FALSE---------------------------------------------------------------
# # # PSEUDO CODE
# # #
# # #     THIS CODE DOES NOT RUN
# # #
# # # Initiate
# # n <- 1
# # eValue <- 1
# #
# # while (eValue < 1/alpha) {
# #   eValue <- saviZTest(x=x[1:n], y=y[1:n],
# #                       designObj=designObj)
# #
# #   if (eValue > 1/alpha) {
# #     "Reject the null"
# #     stop()
# #   } else {
# #     "Increase sample size and test again
# #           at the start of the while loop"
# #     n <- n + 1
# #   }
# # }

## ----eval=FALSE---------------------------------------------------------------
# # Single sample path data set as in the p-value example
# set.seed(2)
# dataGroup1 <- rnorm(n1, mean=muGlobal, sd=sigma)
# set.seed(3)
# dataGroup2 <- rnorm(n2, mean=muGlobal, sd=sigma)
# 
# result <- saviZTest(dataGroup1, dataGroup2[1:i],
#                     designObj=designObj)
# plot(result)

## ----eval=FALSE---------------------------------------------------------------
# mIter <- 1000
# 
# allData <- generateNormalData(c(n1, n2),
#                               muGlobal=muGlobal,
#                               nSim=mIter,
#                               meanDiffTrue=0, seed=1, sigma=sigma)
# allEValues <- matrix(nrow=mIter, ncol=n1)
# 
# # This indicates whether a simulation yielded e > 1/alpha
# eValueOver <- vector("integer", length=mIter)
# 
# # This indicates the first time an experiment yielded e > 1/alpha
# # Default is Inf, which indicates that the e didn't cross 1/alpha
# firstPassageTime <- rep(Inf, times=mIter)
# 
# # Used to vectorise the computations for the the z-statistic
# n1Vector <- 1:n1
# n2Vector <- 1:n2
# nEffVector <- (1/n1Vector+1/n2Vector)^(-1)
# 
# for (sim in 1:mIter) {
#   dataGroup1 <- allData$dataGroup1[sim, ]
#   dataGroup2 <- allData$dataGroup2[sim, ]
# 
#   x1BarVector <- 1/n1Vector*cumsum(dataGroup1)
#   x2BarVector <- 1/n2Vector*cumsum(dataGroup2)
#   zVector <- sqrt(nEffVector)*(x1BarVector - x2BarVector)/sigma
# 
#   for (i in 1:n1) {
#     currentEValue <- saviZTestStat(zVector[i],
#                                    parameter=designObj$parameter,
#                                    n1=n1Vector[i], n2=n2Vector[i],
#                                    sigma=sigma,
#                                    eType=designObj$eType)$eValue
# 
#     allEValues[sim, i] <- currentEValue
# 
#     if (currentEValue > 1/alpha && eValueOver[sim]!=1) {
#       eValueOver[sim] <- 1
#       firstPassageTime[sim] <- i
#     }
#   }
# }

## ----eval=FALSE---------------------------------------------------------------
# numberOfCrossingExperimentsAtTimeN <- integer(n1)
# 
# for (i in 1:n1) {
#   numberOfCrossingExperimentsAtTimeN[i] <- sum(firstPassageTime <= i)
# }
# 
# eValueFalseRejects <- numberOfCrossingExperimentsAtTimeN/mIter
# 
# oldPar <- setSafeStatsPlotOptionsAndReturnOldOnes();
# plot(1:n1, 100*eValueFalseRejects, type="l",
#      xlab="n", ylab="False positive rate (%)", lwd=2,
#      col=eColours[1], ylim=c(0, 5))
# lines(c(1, n1), c(5, 5), lwd=2, lty=2)

## ----eval=FALSE---------------------------------------------------------------
# oldPar <- setSafeStatsPlotOptionsAndReturnOldOnes();
# plot(1:n1, 100*pValueFalseRejects, type="l", xlab="n",
#      ylab="False positive rate (%)",
#      bty="n", lwd=2, ylim=c(0, 40), col=freqColours[1])
# lines(c(1, n1), c(5, 5), lwd=2, lty=2)
# lines(1:n1, 100*subjectiveBfValueFalseRejects,
#       col=bayesColours[1], lwd=2)
# lines(1:n1, 100*eValueFalseRejects, col=eColours[1],
#       lwd=2)

