## ----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)