threshold_values {rfvimptest}R Documentation

Threshold values of a SPRT and SAPT, and calculation of sequential Monte Carlo p-values

Description

Function to produce threshold values for the number of cumulative successes a sequential permutation test of type="SRPT" or type="SAPT" needs to show in order to "keep H0" or "accept H1" early in the sequence of tests.

Usage

threshold_values(
  p0 = 0.06,
  p1 = 0.04,
  alpha = 0.05,
  beta = 0.2,
  A = 0.1,
  B = 10,
  Mmax = 500,
  type = c("SPRT", "SAPT")[1]
)

Arguments

p0

The value of the p-value in the null hypothesis (H0: p = p0) of SPRT and SAPT. Default is 0.06.

p1

The value of the p-value in the alternative hypothesis (H1: p = p1) of SPRT and SAPT. Default is 0.04.

alpha

The significance level of SPRT when p = p0. Also known as type I error. Default is 0.05.

beta

One minus the power of SPRT when p = p1. Also known as type II error. Default is 0.2.

A

The quantity A in the formula of SAPT. Default is 0.1 for a type I error of 0.05. Usually not changed by the user.

B

The quantity B in the formula of SAPT. Default is 10 (1/A) for a type I error of 0.05. Usually not changed by the user.

Mmax

Maximum number of permutations used in each permutation test. Default is 500.

type

Type of the sequential test method. The choices are "SPRT" and "SAPT".

Details

An example of how to calculate sequential Monte Carlo p-values is also given.

The function and examples are intended to support the use of SPRT, SAPT and sequential p-values for any type of modelling approach, even outside the application in rfvimptest::rfvimptest().

Value

List of two vectors, each containing Mmax threshold values for the number of cumulative successes a sequential permutation test of type="SRPT" or type="SAPT" needs to show in order to "keep H0" or "accept H1" early in the sequence of tests.

Author(s)

Alexander Hapfelmeier, Roman Hornung

References

Examples



## Load package:
library("rfvimptest")

## Set seed to obtain reproducible results:
set.seed(1234)


# Load example data
data(hearth2)


# Calculation of a SAPT

# Create critical values for a SAPT
myvals <- threshold_values(p0 = 0.06, p1 = 0.04, A = 0.1, B = 10, Mmax = 200, type = "SAPT")
myvals

# Fit a model to the original data
mod <- glm(I(Class == 2) ~ ., data = hearth2)
summary(mod)

# Derive a statistic of interest from the model. Here, the negative AIC is used as a statistic
# to be maximised.
stat <- -1 * mod$aic

# Perform the permutations, extract the statistics and stop early according to the critical values
myresult <- sapply(setdiff(names(hearth2), "Class"), function(j) {
  permdat <- hearth2
  permstats <- c()
  for (i in 1:length(myvals[[1]])) {
    permdat[, j] <- sample(permdat[, j])
    permstats <- c(permstats, -1 * glm(I(Class == 2) ~ ., data = permdat)$aic)
    if (sum(permstats >= stat) >= myvals[[1]][i] | sum(permstats >= stat) <= myvals[[2]][i]) break
  }
  permstats
})

# Statistics obtained after permutation of the data
myresult

# Number of permutations performed until (early) stopping
sapply(myresult, length)

# Derive the SAPT decisions from the statistics
sapply(myresult, function(x) {
  if (sum(x >= stat) >= myvals[[1]][length(x)]) return("keep H0")
  else if (sum(x >= stat) <= myvals[[2]][length(x)]) return("accept H1")
  else ifelse(sum(x >= stat) / length(myvals[[1]]) <= 0.05, "accept H1", "keep H0")
})


# Calculation of sequential Monte Carlo p-values

h <- 10 # Parameter h of the sequential Monte Carlo p-value calculation.

# Perform the permutations, extract the statistics and stop early when h successes are reached
mypvals <- sapply(setdiff(names(hearth2), "Class"), function(j) {
  permdat <- hearth2
  permstats <- c()
  Mmax <- length(myvals[[1]])
  for (i in 1:Mmax) {
    permdat[, j] <- sample(permdat[, j])
    permstats <- c(permstats, -1 * glm(I(Class == 2) ~ ., data = permdat)$aic)
    d <- sum(permstats >= stat)
    if (d == h) {pval <- d/length(permstats); m <- i; break}
    else if (i == Mmax) {pval <- (d + 1)/(Mmax + 1); m <- i}
  }
  c(pval, m)
})

# p-values and number of permutations performed until (early) stopping
rownames(mypvals) <- c("p-value", "m")
t(mypvals)




[Package rfvimptest version 0.1.4 Index]