Nothing
######################################################################################################################
# Function: MultipleSequenceGatekeepingAdj.
# Argument: rawp, Raw p-value.
# par, List of procedure parameters: vector of family (1 x m) Vector of component procedure labels ('BonferroniAdj.global' or 'HolmAdj.global' or 'HochbergAdj.global' or 'HommelAdj.global') (1 x nfam) Vector of truncation parameters for component procedures used in individual families (1 x nfam)
# Description: Computation of adjusted p-values for gatekeeping procedures based on the modified mixture methods (ref Dmitrienko et al. (2014))
MultipleSequenceGatekeepingAdj = function(rawp, par) {
# Determine the function call, either to generate the p-value or to return description
call = (par[[1]] == "Description")
if (any(call == FALSE) | any(is.na(call))) {
# Error check
if (is.null(par[[2]]$family)) stop("Analysis model: Multiple sequence gatekeeping procedure: Hypothesis families must be specified.")
if (is.null(par[[2]]$proc)) stop("Analysis model: Multiple sequence gatekeeping procedure: Procedures must be specified.")
if (is.null(par[[2]]$gamma)) stop("Analysis model: Multiple sequence gatekeeping procedure: Gamma must be specified.")
# Number of p-values
nhyp = length(rawp)
# Extract the vector of family (1 x m)
family = par[[2]]$family
# Number of families in the multiplicity problem
nfam = length(family)
# Number of null hypotheses per family
nperfam = nhyp/nfam
# Extract the vector of procedures (1 x m)
proc = paste(unlist(par[[2]]$proc), ".global", sep = "")
# Extract the vector of truncation parameters (1 x m)
gamma = unlist(par[[2]]$gamma)
# Simple error checks
if (nhyp != length(unlist(family)))
stop("Multiple-sequence gatekeeping adjustment: Length of the p-value vector must be equal to the number of hypothesis.")
if (length(proc) != nfam)
stop("Multiple-sequence gatekeeping adjustment: Length of the procedure vector must be equal to the number of families.") else {
for (i in 1:nfam) {
if (proc[i] %in% c("BonferroniAdj.global", "HolmAdj.global", "HochbergAdj.global", "HommelAdj.global") == FALSE)
stop("Multiple-sequence gatekeeping adjustment: Only Bonferroni (BonferroniAdj), Holm (HolmAdj), Hochberg (HochbergAdj) and Hommel (HommelAdj) component procedures are supported.")
}
}
if (length(gamma) != nfam)
stop("Multiple-sequence gatekeeping adjustment: Length of the gamma vector must be equal to the number of families.") else {
for (i in 1:nfam) {
if (gamma[i] < 0 | gamma[i] > 1)
stop("Multiple-sequence gatekeeping adjustment: Gamma must be between 0 (included) and 1 (included).") else if (proc[i] == "bonferroni.global" & gamma[i] != 0)
stop("Multiple-sequence gatekeeping adjustment: Gamma must be set to 0 for the global Bonferroni procedure.")
}
}
# Number of intersection hypotheses in the closed family
nint = 2^nhyp - 1
# Construct the intersection index sets (int_orig) before the logical restrictions are applied. Each row is a vector of binary indicators (1 if the hypothesis is
# included in the original index set and 0 otherwise)
int_orig = matrix(0, nint, nhyp)
for (i in 1:nhyp) {
for (j in 0:(nint - 1)) {
k = floor(j/2^(nhyp - i))
if (k/2 == floor(k/2))
int_orig[j + 1, i] = 1
}
}
# Construct the intersection index sets (int_rest) and family index sets (fam_rest) after the logical restrictions are applied. Each row is a vector of binary
# indicators (1 if the hypothesis is included in the restricted index set and 0 otherwise)
int_rest = int_orig
fam_rest = matrix(1, nint, nhyp)
for (i in 1:nint) {
for (j in 1:(nfam - 1)) {
for (k in 1:nperfam) {
# Index of the current null hypothesis in Family j
m = (j - 1) * nperfam + k
# If this null hypothesis is included in the intersection hypothesis all dependent null hypotheses must be removed from the intersection hypothesis
if (int_orig[i, m] == 1) {
for (l in 1:(nfam - j)) {
int_rest[i, m + l * nperfam] = 0
fam_rest[i, m + l * nperfam] = 0
}
}
}
}
}
# Number of null hypotheses from each family included in each intersection before the logical restrictions are applied
korig = matrix(0, nint, nfam)
# Number of null hypotheses from each family included in the current intersection after the logical restrictions are applied
krest = matrix(0, nint, nfam)
# Number of null hypotheses from each family after the logical restrictions are applied
nrest = matrix(0, nint, nfam)
# Compute korig, krest and nrest
for (j in 1:nfam) {
# Index vector in the current family
# index = which(family == j)
index = family[[j]]
korig[, j] = apply(as.matrix(int_orig[, index]), 1, sum)
krest[, j] = apply(as.matrix(int_rest[, index]), 1, sum)
nrest[, j] = apply(as.matrix(fam_rest[, index]), 1, sum)
}
# Vector of intersection p-values
pint = rep(1, nint)
# Matrix of component p-values within each intersection
pcomp = matrix(0, nint, nfam)
# Matrix of family weights within each intersection
c = matrix(0, nint, nfam)
# P-value for each hypothesis within each intersection
p = matrix(0, nint, nhyp)
# Compute the intersection p-value for each intersection hypothesis
for (i in 1:nint) {
# Compute component p-values
for (j in 1:nfam) {
# Consider non-empty restricted index sets
if (krest[i, j] > 0) {
# Restricted index set in the current family
int = int_rest[i, family[[j]]]
# Set of p-values in the current family
pv = rawp[family[[j]]]
# Select raw p-values included in the restricted index set
pselected = pv[int == 1]
# Total number of hypotheses used in the computation of the component p-value
tot = nrest[i, j]
pcomp[i, j] = do.call(proc[j], list(pselected, tot, gamma[j]))
} else if (krest[i, j] == 0)
pcomp[i, j] = 1
}
# Compute family weights
c[i, 1] = 1
for (j in 2:nfam) {
c[i, j] = c[i, j - 1] * (1 - errorfrac(krest[i, j - 1], nrest[i, j - 1], gamma[j - 1]))
}
# Compute the intersection p-value for the current intersection hypothesis
pint[i] = pmin(1, min(ifelse(c[i,]>0, pcomp[i, ]/c[i, ], NA), na.rm = TRUE))
# Compute the p-value for each hypothesis within the current intersection
p[i, ] = int_orig[i, ] * pint[i]
}
# Compute adjusted p-values
adjustedp = apply(p, 2, max)
result = adjustedp
}
else if (call == TRUE) {
family = par[[2]]$family
nfam = length(family)
proc = unlist(par[[2]]$proc)
gamma = unlist(par[[2]]$gamma)
test.id=unlist(par[[3]])
proc.par = data.frame(nrow = nfam, ncol = 4)
for (i in 1:nfam){
proc.par[i,1] = i
proc.par[i,2] = paste0("{",paste(test.id[family[[i]]], collapse = ", "),"}")
proc.par[i,3] = proc[i]
proc.par[i,4] = gamma[i]
}
colnames(proc.par) = c("Family", "Tests", "Component procedure", "Truncation parameter")
result=list(list("Multiple-sequence gatekeeping"),list(proc.par))
}
return(result)
}
# End of MultipleSequenceGatekeepingAdj
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.