Nothing
paired.reject.region <-
function(N, alternative=c("two.sided", "less", "greater"), alpha=0.05,
npNumbers=100, np.interval=FALSE, beta=0.001,
method=c("uam", "ucm", "uamcc", "csm", "cm", "am", "amcc"),
tsmethod=c("square", "central"), delta=0, convexity=TRUE, useStoredCSM=TRUE){
alternative <- match.arg(tolower(alternative), c("two.sided", "less", "greater"))
method <- convertMethod(method)
tsmethod <- match.arg(tolower(tsmethod), c("square", "central"))
#Perform several checks
checkPairedParam(N=N, alternative=alternative, alpha=alpha, npNumbers=npNumbers, np.interval=np.interval, beta=beta,
method=method, tsmethod=tsmethod, delta=delta, convexity=convexity, useStoredCSM=useStoredCSM)
# if two-sided and central, easiest to simply find rejection regions when greater and less and dividing alpha by 2
if (alternative == "two.sided" && tsmethod == "central") {
rejectRegionUpper <- paired.reject.region(N=N, alternative="less", alpha=alpha/2,
npNumbers=npNumbers, np.interval=np.interval, beta=beta,
method=method, tsmethod=tsmethod, delta=delta, convexity=convexity,
useStoredCSM=useStoredCSM)
rejectRegionLower <- paired.reject.region(N=N, alternative="greater", alpha=alpha/2,
npNumbers=npNumbers, np.interval=np.interval, beta=beta,
method=method, tsmethod=tsmethod, delta=delta, convexity=convexity,
useStoredCSM=useStoredCSM)
rejectRegion <- rejectRegionUpper + rejectRegionLower
if (!all(rejectRegion %in% c(0,1,NA))) { stop("Check rejection region logic") }
return(rejectRegion)
}
# Instead of trying to solve both "less" and "greater", just have code only work for "less"
# Will transpose at the end
alternativeTemp <- alternative
swapFlg <- (alternative == "greater" || (alternative=="two.sided" && delta < 0))
if (swapFlg) {
delta <- -delta
if (alternative == "greater") { alternative <- "less" }
}
# Each test has a different way of minimizing the computation time to calculate the rejection region.
# Most tests have a test statistic which shows the order of the tables considered. This can be utilized to greatly speed up time.
# The CSM test is essentially just calculating the p-value except stopping at alpha
# The conditional mcnemar test doesn't have an ordered test statistic, so calculate sequentially with convexity property
# The np.interval approach is a special case since the NP interval depends on the number of successes, which makes it not completely
# convex. The precise way is to calculate a p-value for all tables, but in almost all cases the convexity property holds
if (!np.interval && !(method %in% c("csm", "cm"))) {
if (delta == 0) { int <- seq(0.00001,.49999,length=npNumbers)
} else if (delta > 0) { int <- seq(0.00001, (1-delta)/2 - 0.00001, length=npNumbers)
} else if (delta < 0) { int <- seq(-delta + 0.00001, (1 - delta)/2 - 0.00001, length=npNumbers)}
lookupArray <- trinomCalc(N, int, delta)
# For the tests with a test statistic:
TX <- switch(method,
"uam" = mcnemar_TX(NULL, N, delta=delta, CC=FALSE),
"ucm" = mcnemar.2x2(NULL, N, alternative=alternative, pval1ties = TRUE),
"uamcc" = mcnemar_TX(NULL, N, delta=delta, CC=TRUE),
"am" = mcnemar_TX(NULL, N, delta=delta, CC=FALSE),
"amcc" = mcnemar_TX(NULL, N, delta=delta, CC=TRUE))
# Instead of starting with the most extreme table and adding sequentially, start in the middle and move up or down
# depending on whether we have a more extreme table. This can greatly reduce the number of p-values calculated
# If two-sided test, then first table will be the largest absolute value of the test statistic.
# Take the negative value to match the "less" alternative
if (alternative == "two.sided" && method != "ucm") {TX[, 3] <- -abs(TX[,3])}
TX <- signif(TX[order(TX[,3]), ], 12)
TX <- cbind(TX, NA)
# We can just cross out all test statistics >=0 since this p-value would be >= 0.5.
if (method %in% c("uam", "uamcc", "am", "amcc") && delta==0) { TX[TX[,3] >= 0, 4] <- FALSE }
# For ucm, we know significant if Conditional McNemar's p-value < alpha
if (method == "ucm" && delta==0) { TX[TX[,3] <= alpha, 4] <- TRUE }
moreExtremeTbls <- searchExtremePaired(TX = TX, N = N, alternative = alternative, method = method, int = int, delta = delta,
alpha = alpha, lookupArray = lookupArray)
rejectRegion <- matrix(0, N+1, N+1)
rejectRegion <- lowerMatVal(rejectRegion, NA)
if (nrow(moreExtremeTbls) > 0) { #If at least one more extreme table
for (i in 1:nrow(moreExtremeTbls)) {
rejectRegion[moreExtremeTbls[i, 1] + 1, moreExtremeTbls[i, 2] + 1] <- 1
}
}
rownames(rejectRegion) <- 0:N
colnames(rejectRegion) <- 0:N
if (swapFlg){ rejectRegion <- t(rejectRegion) }
return(rejectRegion)
} else if (method == "csm") {
# Use stored ordering matrix if available
if ( N <= 200 && delta == 0 && useStoredCSM) {
if (!requireNamespace("ExactData", quietly = TRUE)) {
stop(paste("ExactData R package must be installed when useStoredCSM=TRUE. To install ExactData R package, run:",
"`install.packages('ExactData', repos='https://pcalhoun1.github.io/drat/', type='source')`"))
}
# Note: delta=0
int <- seq(0.00001, .49999, length=npNumbers)
lookupArray <- trinomCalc(N, int, delta)
if (alternative == "two.sided") { orderMat <- ExactData::orderCSMPairedMatTwoSided[[N]]
} else { orderMat <- ExactData::orderCSMPairedMatOneSided[[N]] }
orderVec <- as.vector(orderMat)[!is.na(orderMat)]
TX <- matrix(cbind(which(!is.na(orderMat), arr.ind = TRUE) - 1, orderVec), ncol=3, dimnames = NULL)
TX <- signif(TX[order(TX[,3]), ], 12)
TX <- cbind(TX, NA)
moreExtremeTbls <- searchExtremePaired(TX = TX, N = N, alternative = alternative, method = method, int = int, delta = delta,
alpha = alpha, lookupArray = lookupArray)
rejectRegion <- matrix(0, N+1, N+1)
rejectRegion <- lowerMatVal(rejectRegion, NA)
if (nrow(moreExtremeTbls) > 0) { #If at least one more extreme table
for (i in 1:nrow(moreExtremeTbls)) {
rejectRegion[moreExtremeTbls[i, 1] + 1, moreExtremeTbls[i, 2] + 1] <- 1
}
}
rownames(rejectRegion) <- 0:N
colnames(rejectRegion) <- 0:N
} else {
if (delta == 0) { int <- seq(0.00001,.49999,length=npNumbers)
} else if (delta > 0) { int <- seq(0.00001, (1-delta)/2 - 0.00001, length=npNumbers)
} else if (delta < 0) { int <- seq(-delta + 0.00001, (1 - delta)/2 - 0.00001, length=npNumbers)}
rejectRegion <- moreExtremePairedCSM(data = NULL, N = N, alternative = alternative,
int = int, doublePvalue = FALSE, delta = delta, reject.alpha = alpha)$moreExtremeMat
}
if (swapFlg) { rejectRegion <- t(rejectRegion) }
return(rejectRegion)
} else if (method == "cm" || np.interval) {
rejectRegion <- matrix(NA, N+1, N+1)
# Lower diagonal always starts at 0
rejectRegion[round((row(rejectRegion)-1)/(nrow(rejectRegion)-1), digits=10) >= round((col(rejectRegion)-1)/(ncol(rejectRegion)-1), digits=10)] <- 0
rejectRegion <- lowerMatVal(rejectRegion)
# Conditional McNemar's test has the convexity property, so can determine the entire rejection region through the boundary.
# Essentially, start at first row and move through each column until table rejects. Then move to next row
# For a given alpha and sample size, the np.interval approach almost always has the convexity property. However, there are
# very rare instances when the convexity property doesn't hold and one has to consider all possible tables. Add an input parameter
# for whether user is wanting to assume convexity (very fast), or do it correctly by going through all possible tables (very slow)
for (i in 0:N) { #Go through each row
# Find first column with NA
startJ <- which(is.na(rejectRegion[i+1, ]))[1] - 1
if (!is.na(startJ)) {
for (j in startJ:(N-i)) {
tables <- matrix(c(N-i-j, i, j, 0), 2, 2, byrow=TRUE)
if (method=="cm") {
rejectRegionTemp <- (mcnemar.2x2(tables, N, alternative=alternative, pval1ties = FALSE)[,3] <= alpha)
} else {
rejectRegionTemp <- pairedCode(tables, alternative=alternative, npNumbers=npNumbers,
np.interval=np.interval, beta=beta, method=method, tsmethod=tsmethod, to.plot=FALSE,
ref.pvalue=FALSE, delta=delta, reject.alpha=alpha, useStoredCSM=useStoredCSM)
}
if (rejectRegionTemp) {
if (method=="cm" || convexity) {
# If p-value <= alpha, then know the remaining columns in the row is more extreme
rejectRegion[i+1, (j+1):(N+1)] <- 1
if (alternative=="two.sided") { rejectRegion[(j+1):(N+1), i+1] <- 1 }
break
} else {
rejectRegion[i+1, j+1] <- 1
if (alternative=="two.sided") { rejectRegion[j+1, i+1] <- 1 }
}
} else {
if (method=="cm" || convexity) {
# If p-value > alpha, then know the remaining rows in the column is less extreme
rejectRegion[is.na(rejectRegion[ , j+1]), j+1] <- 0
} else { rejectRegion[i+1, j+1] <- 0 }
}
}
}
}
rownames(rejectRegion) <- 0:N
colnames(rejectRegion) <- 0:N
rejectRegion <- lowerMatVal(rejectRegion, NA)
if (swapFlg){ rejectRegion <- t(rejectRegion) }
return(rejectRegion)
}
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.