Nothing
#### intern function
# method from hierarchical false discovery rate-controlling method
# fdr control of family: q
# fdr control of full tree: q * delta * 2 (delta < 1.44)
# fdr control of outer node: q * L * delta * 2 (delta < 1.44)
hierarchicalFDRTesting <- function(hierMat, group, grouplm, X, y, test = partialFtest) {
# root of the tree
indRoot <- findRoot2(hierMat)
# output container
pvalues <- rep(0, nrow(hierMat))
adjPvalues <- rep(0, nrow(hierMat))
# all Leaves in the hierarchical matrix
allLeaves <- which(rowSums(hierMat) == 1)
familyToTest <- list()
familyToTest[[1]] <- indRoot
continue <- TRUE
# hierarchical testing
while (continue) {
familyToTestTemp <- list()
compteur <- 1
continue <- FALSE
# current family at the level
for (i in seq_along(familyToTest))
{
# test the group of a family
for (gr in familyToTest[[i]])
{
# group included in gr
subGroup <- setdiff(which(hierMat[gr, ]), gr)
# if no subgroups, gr is a leaf
if (length(subGroup) == 0) {
subGroup <- gr
}
# only group corresponding to leaves
subLeaves <- intersect(subGroup, allLeaves)
toTest0 <- which(grouplm %in% group[subLeaves])
pvalues[gr] <- test(X, y, toTest0)
}
# BH correction for the family
adjPvalues[familyToTest[[i]]] <- p.adjust(pvalues[familyToTest[[i]]], "BH")
# if selected, we look for the children to test
for (j in seq_along(familyToTest[[i]]))
{
child <- children(familyToTest[[i]][j], hierMat)
if (length(child) > 0) {
familyToTestTemp[[compteur]] <- child
compteur <- compteur + 1
continue <- TRUE
}
} # end for selected group
} # end for family
familyToTest <- familyToTestTemp
} # end while hierarchy
return(list(pvalues = pvalues, adjPvalues = adjPvalues, groupId = as.numeric(colnames(hierMat))))
}
#'
#' Apply hierarchical test for each hierarchy, and test external variables for FDR control at level alpha
#'
#' @title Hierarchical testing with FDR control
#'
#' @param X original data
#' @param y associated response
#' @param group vector with index of groups. group[i] contains the index of the group of the variable var[i].
#' @param var vector with the variables contained in each group. group[i] contains the index of the group of the variable var[i].
#' @param test function for testing the nullity of a group of coefficients in linear regression.
#' The function has 3 arguments: \code{X}, the design matrix, \code{y}, response, and \code{varToTest},
#' a vector containing the indices of the variables to test. The function returns a p-value
#' @param addRoot If TRUE, add a common root containing all the groups
#'
#' @return a list containing:
#' \describe{
#' \item{pvalues}{pvalues of the different test (without correction)}
#' \item{adjPvalues}{adjusted pvalues}
#' \item{groupId}{Index of the group}
#' \item{hierMatrix}{Matrix describing the hierarchical tree.}
#' }
#'
#' @details
#' Version of the hierarchical testing procedure of Yekutieli for MLGL output. You can use th \link{selFDR} function
#' to select groups at a desired level alpha.
#'
#' @examples
#' set.seed(42)
#' X <- simuBlockGaussian(50, 12, 5, 0.7)
#' y <- X[, c(2, 7, 12)] %*% c(2, 2, -2) + rnorm(50, 0, 0.5)
#' res <- MLGL(X, y)
#' test <- hierarchicalFDR(X, y, res$group[[20]], res$var[[20]])
#' @references Yekutieli, Daniel. "Hierarchical False Discovery Rate-Controlling Methodology."
#' Journal of the American Statistical Association 103.481 (2008): 309-16.
#'
#' @seealso \link{selFDR}, \link{hierarchicalFWER}
#'
#' @export
hierarchicalFDR <- function(X, y, group, var, test = partialFtest, addRoot = FALSE) {
# check parameters
.checkhierarchicalFWER(X, y, group, var, test, TRUE)
# hierarchical matrix
hierInfo <- groupHier(group, var, addRoot)
# lm with leaves represented by their first principal component
reslm <- acpOLStest(X, y, hierInfo$grouplm, hierInfo$varlm)
# new hierMat with complementary group
# hierMatTot <- compHierMatTot(hierInfo)
hierMatTot <- hierInfo$hierTot
groupId <- colnames(hierMatTot)
# hierarchical testing on the tree of indRoot
out <- hierarchicalFDRTesting(hierMatTot, groupId, reslm$group, reslm$newdata, y, test = partialFtest)
out$group <- hierInfo$groupTot
out$var <- hierInfo$varTot
out$hierMatrix <- hierMatTot
return(out)
}
#'
#' Select groups from hierarchical testing procedure with FDR control (\link{hierarchicalFDR})
#'
#' @title Selection from hierarchical testing with FDR control
#'
#' @param out output of \link{hierarchicalFDR} function
#' @param alpha control level for test
#' @param global if FALSE the provided alpha is the desired level control for each family.
#' @param outer if TRUE, the FDR is controlled only on outer node (rejected groups without rejected children).
#' If FALSE, it is controlled on the full tree.
#'
#' @return a list containing:
#' \describe{
#' \item{toSel}{vector of boolean. TRUE if the group is selected}
#' \item{groupId}{Names of groups}
#' \item{local.alpha}{control level for each family of hypothesis}
#' \item{global.alpha}{control level for the tree (full tree or outer node)}
#' }
#'
#' @details
#' See the reference for mode details about the method.
#'
#' If each family is controlled at a level alpha, we have the following control:
#' FDR control of full tree: alpha * delta * 2 (delta = 1.44)
#' FDR control of outer node: alpha * L * delta * 2 (delta = 1.44)
#'
#' @examples
#' set.seed(42)
#' X <- simuBlockGaussian(50, 12, 5, 0.7)
#' y <- X[, c(2, 7, 12)] %*% c(2, 2, -2) + rnorm(50, 0, 0.5)
#' res <- MLGL(X, y)
#' test <- hierarchicalFDR(X, y, res$group[[20]], res$var[[20]])
#' sel <- selFDR(test, alpha = 0.05)
#' @references Yekutieli, Daniel. "Hierarchical False Discovery Rate-Controlling Methodology."
#' Journal of the American Statistical Association 103.481 (2008): 309-16.
#'
#' @seealso \link{hierarchicalFDR}
#'
#' @export
selFDR <- function(out, alpha = 0.05, global = TRUE, outer = TRUE) {
## check arguments
.checkselFDR(out$adjPvalues, out$hierMatrix, alpha, global, outer)
delta <- 1.44
# number of levels of the hierarchy
L <- ifelse(outer, numberLevels(out$hierMatrix), 1)
# compute local and global alpha
local.alpha <- ifelse(global, alpha / (delta * 2 * L), alpha)
global.alpha <- ifelse(global, alpha, min(alpha * (delta * 2 * L), 1))
# if global = 1, then use local = 1
local.alpha <- ifelse(global.alpha == 1, 1, local.alpha)
# if one level, simple BH adjust
local.alpha <- ifelse(outer & (L == 1), alpha, local.alpha)
global.alpha <- ifelse(outer & (L == 1), alpha, global.alpha)
#
toSel <- rep(FALSE, length(out$adjPvalues))
# indices of groups at the top of hierarchy
family <- findRoot2(out$hierMatrix)
continue <- TRUE
while (continue) {
continue <- FALSE
# select group with adjusted p-values <= local.alpha
toSel[family] <- (out$adjPvalues[family] <= local.alpha)
# find children of selected group
if (any(!is.na(toSel[family])) & any(toSel[family])) {
ind <- which(toSel[family])
newfamily <- c()
for (i in seq_along(ind)) {
newfamily <- c(newfamily, children(family[ind[i]], out$hierMatrix))
}
family <- newfamily
continue <- (length(family) > 0)
}
} # fin while
# we select only outer node
if (outer) {
toSel[which(toSel)[which(rowSums(out$hierMatrix[toSel, toSel, drop = FALSE]) > 1)]] <- FALSE
}
return(list(toSel = toSel, groupId = as.numeric(colnames(out$hierMatrix)), local.alpha = local.alpha,
global.alpha = global.alpha))
}
## check parameters of selFWER function
.checkselFDR <- function(adjPvalues, hierMatrix, alpha, global, outer) {
.checkselFWER(adjPvalues, hierMatrix, alpha)
# check if global is a boolean
if (length(global) != 1) {
stop("global must be a boolean.")
}
if (!is.logical(global)) {
stop("global must be a boolean.")
}
# check if outer is a boolean
if (length(outer) != 1) {
stop("outer must be a boolean.")
}
if (!is.logical(outer)) {
stop("outer must be a boolean.")
}
invisible(return(NULL))
}
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.