Confidence curves for structured hypotheses

G. Durand, G. Blanchard, P. Neuvial, E. Roquain

2022-11-11

Source: vignettes/confidenceCurves_localized.Rmd

Introduction

The goal of this vignette is to illustrate the post hoc bounds on the number of true/false positives proposed in Durand et al. (2020) for localized signals. More specifically, we reproduce one of the plots of Figure 12 in Durand et al. (2020) using the R package sanssouci. We will explicitly quote Durand et al. (2020) when relevant.

Objective

We consider \(m\) null ordered hypotheses partitioned in intervals of size \(s\). For simplicity we set the number of intervals to be a power of 2: \(m = s 2^q\) for some integer \(q\):

s <- 100
q <- 7
m <- s*2^q

Therefore, we have \(m = 1.28\times 10^{4}\).Our goal is to compare three post hoc bounds. These bounds are obtained by interpolation from a reference family where the amount of signal is estimated by probabilistic inequalities, following the general principle laid down by Blanchard, Neuvial, and Roquain (2020), and they differ by the choice of the reference family:

  • “Simes”: the bound derived from the Simes (1986) inequality as proposed by Goeman and Solari (2011) and further studied by Blanchard, Neuvial, and Roquain (2020). This bound was introduced in a context where the signal is not localized.

  • “Tree” and “Partition”: two bounds derived from the DKWM inequality (Dvoretzky, Kiefer, and Wolfowitz (1956), Massart (1990)), as proposed in Durand et al. (2020). For the “Partition” bound, the reference family is the original partition \((P_k)_k\) of the \(m\) null hypotheses into \(K=2^q\) intervals. For the “Tree” bound, the reference family is the perfect binary tree whose leaves are the elements of the original partition.

Settings

We define the following numerical parameters, which characterize the true/false null hypothesis configuration considered in Section 5 of Durand et al. (2020):

K1 <- 8
r <- 0.9
m1 <-  r*K1*s
barmu <- 3

More precisely, quoting Durand et al. (2020):

the false null hypotheses are contained in \(P_k\) for \(1 \leq k \leq K_1\), for some fixed value of \(K_1\). The quantity \(r\) is defined similarly as in (20), as the fraction of false null hypotheses in those \(P_k\), and is set to \(r =0.9\). All of the other partition pieces only contain true null hypotheses.

We start by creating a binary tree structure and generating the signal:

obj <- SansSouciDyadic(m, leaf_size = s, direction = "top-down")
mu <- gen.mu.leaves(m = m, K1 = K1, d = r, grouped = TRUE, setting = "const", 
                    barmu = barmu, leaf_list = obj$input$leaves)

This construction is illustrated in the figure below (Figure 11 in Durand et al. (2020)) in the case where \(q=3\).

Figure 11 in Durand et al.

Figure 11 in Durand et al.

We generate \(p\)-values according to the simulation setup in Durand et al. (2020):

The true null \(p\)-values are distributed as i.i.d. \(\mathcal{N}(0,1)\), and false null \(p\)-values are distributed as i.i.d. \(\mathcal{N}(\bar{\mu}, 1)\), where [\(\bar{\mu}= 3\)].

pvalues <- gen.p.values(m = m, mu = mu, rho = 0)

Calculate confidence curves

The confidence level for post hoc inference is set to \(\alpha = 0.05\).

alpha <- 0.05

Below, we will be considering confidence curves of the form \((k, V(S_k))_{1 \leq k \leq m}\), where \(S_k\) is the set of the \(k\) smallest \(p\)-values (regardless of the ordering given by the partition). Note that focusing on such sets is a priori favorable to the Simes bound, for which the elements of the reference family are among the \(S_k\).

ord <- order(pvalues)
idxs <- seq_len(nHyp(obj))
res <- list()

1- True number of false positives

The true number of false positives will be called “Oracle” bound in the plots below.

obj$input$truth <- as.numeric(mu != 0)
res_Oracle <- fit(obj, alpha = alpha, p.values = pvalues, family = "Oracle")
res[["Oracle"]] <- predict(res_Oracle, what = "FP", all = TRUE)
res[["Oracle"]]$method <- "Oracle"

2- Simes-based confidence curve

Here we use the Simes (1986) inequality to bound the number of false positives in each node of the tree, as proposed by Goeman and Solari (2011) and further studied by Blanchard, Neuvial, and Roquain (2020), both in a context where the signal is not localized.

res_Simes <- fit(obj, alpha = alpha, p.values = pvalues, family = "Simes")
res[["Simes"]] <- predict(res_Simes, what = "FP", all = TRUE)
res[["Simes"]]$method <- "Simes"

3- DKWM-based confidence curve

Here we use the DKWM inequality (Dvoretzky, Kiefer, and Wolfowitz (1956), Massart (1990)) to bound the number of false positives in each node of the tree, as suggested in Durand et al. (2020).

res_DKWM <- fit(obj, alpha, p.values = pvalues, family = "DKWM")
res[["Tree"]] <- predict(res_DKWM, what = "FP", n_out = 40)
res[["Tree"]]$method <- "Tree"
res_DKWM_part <- fit(obj, alpha, p.values = pvalues, 
                family = "DKWM", flavor = "partition")
res[["Partition"]] <- predict(res_DKWM_part, what = "FP", n_out = 40)
res[["Partition"]]$method <- "Partition"

Plot confidence curves 

library("ggplot2")
dat <- Reduce(rbind, res)
lvls <- c("Oracle", "Partition", "Simes", "Tree", "Hybrid")
cols <- RColorBrewer::brewer.pal(length(lvls), "Set1")
names(cols) <- lvls

Upper bound on the number of false positives 

xymax <- m1;
pV <- ggplot(dat, aes(x, bound, colour = method)) + 
    geom_line() +
    ylab("Upper bound on the number of false positives") +
    xlab("sorted hypotheses") +
    scale_colour_manual(values = cols)
pV

The “Tree” and “Partition” bounds are sharper than the “Simes” bound as soon as we are considering “large” sets of hypotheses. The fact that the “Tree” and “Partition” bounds are not as sharp as the “Simes” bound for the first hundred of hypotheses can be explained by our choice of the ordering of the null hypotheses in the sets \(S_k\), which as discussed above is favorable to the “Simes” bound.

Zooming on the first 720 null hypotheses (in the order of the \(p\)-values):

pV + coord_cartesian(xlim = c(1, xymax),
                     ylim = c(0, xymax))

Lower bound on the number of true positives

The same information can be displayed as a lower bound on the number of true positives, defined for any \(S \subset \{1 \dots m\}\) by \(|S| - V(S)\):

dat$S <- dat$x - dat$bound
xmax <- m1;
ymax <- max(dat$S);
pS <- ggplot(dat, aes(x, S, colour = method)) + 
    geom_line() +
    ylab("Lower bound on the number of true positives") +
    xlab("sorted hypotheses") +
    scale_colour_manual(values = cols)
pS

Zooming in the first 720 null hypotheses (in the order of the \(p\)-values), we recover the middle plot in Figure 12 of Durand et al. (2020).

pS + xlim(1, xmax) + ylim(0, ymax)

Session information 

## R version 4.2.2 (2022-10-31)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 20.04.5 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.9.0
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.9.0
## 
## locale:
##  [1] LC_CTYPE=C.UTF-8       LC_NUMERIC=C           LC_TIME=C.UTF-8       
##  [4] LC_COLLATE=C.UTF-8     LC_MONETARY=C.UTF-8    LC_MESSAGES=C.UTF-8   
##  [7] LC_PAPER=C.UTF-8       LC_NAME=C              LC_ADDRESS=C          
## [10] LC_TELEPHONE=C         LC_MEASUREMENT=C.UTF-8 LC_IDENTIFICATION=C   
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] ggplot2_3.4.0    sanssouci_0.12.8
## 
## loaded via a namespace (and not attached):
##  [1] Rcpp_1.0.9          RColorBrewer_1.1-3  pillar_1.8.1       
##  [4] bslib_0.4.1         compiler_4.2.2      jquerylib_0.1.4    
##  [7] highr_0.9           tools_4.2.2         digest_0.6.30      
## [10] tibble_3.1.8        jsonlite_1.8.3      evaluate_0.18      
## [13] memoise_2.0.1       lifecycle_1.0.3     gtable_0.3.1       
## [16] lattice_0.20-45     pkgconfig_2.0.3     rlang_1.0.6        
## [19] Matrix_1.5-3        cli_3.4.1           yaml_2.3.6         
## [22] pkgdown_2.0.6       xfun_0.34           fastmap_1.1.0      
## [25] withr_2.5.0         stringr_1.4.1       knitr_1.40         
## [28] vctrs_0.5.0         desc_1.4.2          generics_0.1.3     
## [31] fs_1.5.2            sass_0.4.2          systemfonts_1.0.4  
## [34] rprojroot_2.0.3     grid_4.2.2          glue_1.6.2         
## [37] R6_2.5.1            textshaping_0.3.6   fansi_1.0.3        
## [40] rmarkdown_2.18      farver_2.1.1        purrr_0.3.5        
## [43] magrittr_2.0.3      scales_1.2.1        htmltools_0.5.3    
## [46] matrixStats_0.62.0  matrixTests_0.1.9.1 colorspace_2.0-3   
## [49] labeling_0.4.2      ragg_1.2.4          utf8_1.2.2         
## [52] stringi_1.7.8       munsell_0.5.0       cachem_1.0.6

References

Blanchard, Gilles, Pierre Neuvial, and Etienne Roquain. 2020. “Post Hoc Confidence Bounds on False Positives Using Reference Families.” Annals of Statistics 48 (3): 1281–1303. https://projecteuclid.org/euclid.aos/1594972818.

Durand, Guillermo, Gilles Blanchard, Pierre Neuvial, and Etienne Roquain. 2020. “Post hoc false positive control for structured hypotheses.” Scandinavian Journal of Statistics. https://hal.archives-ouvertes.fr/hal-01829037.

Dvoretzky, Aryeh, Jack Kiefer, and Jacob Wolfowitz. 1956. “Asymptotic Minimax Character of the Sample Distribution Function and of the Classical Multinomial Estimator.” The Annals of Mathematical Statistics 27 (3): 642–69.

Goeman, Jelle J, and Aldo Solari. 2011. “Multiple Testing for Exploratory Research.” Statistical Science 26 (4): 584–97.

Massart, Pascal. 1990. “The tight constant in the Dvoretzky-Kiefer-Wolfowitz inequality.” The Annals of Probability 18 (3): 1269–83.

Simes, R J. 1986. “An improved Bonferroni procedure for multiple tests of significance.” Biometrika 73 (3): 751–54.