Screening before selecting

A genomics biobank holds hundreds of thousands of variables and hundreds of phenotypes. Running the full T-Rex calibration on every phenotype is costly, and for most phenotypes the answer is that nothing is selected. Screen-T-Rex [1] is a single-pass version of the random experiments without the L-loop and the T-loop: \(K\) experiments with \(p\) dummies each, stopped at the first dummy, and an estimate of the false discovery rate (FDR) read off the vote. It reports an estimated FDR instead of controlling a target, so it is a fast first pass. The biobank screening wraps it in a decision rule per phenotype: when the estimated FDR of the screening falls into an acceptance window, the screening result is kept, otherwise the full T-Rex selector runs as a fallback.

Screen-T-Rex

The screening selector takes no target FDR.

set.seed(2024)
n <- 100
p <- 300
true_support <- c(3, 50, 120, 200, 280)
beta <- numeric(p)
beta[true_support] <- 2
X <- matrix(rnorm(n * p), n, p)
y <- as.numeric(X %*% beta + rnorm(n))

scr <- TRexScreeningSelector$new(X, y, seed = 3, verbose = FALSE,
                                 control = trex_control(K = 20))
scr$select()
scr$selected_indices
#> [1]   3  25  50 120 200 231 280
scr$estimated_FDR
#> [1] 0.1428571

beta_mat holds the coefficients of the original variables at the stop of each experiment (one column per experiment) and dummy_betas those of the dummies, the material the estimate is read off. The L-loop strategy is forced to "STANDARD" here, fresh dummies per experiment.

dim(scr$beta_mat)
#> [1] 300   1
dim(scr$dummy_betas)
#> [1] 20  1

A bootstrap confidence interval

use_bootstrap_CI = TRUE resamples the experiments R_boot times and reports the upper confidence bound of the estimated FDR instead of the point estimate, at confidence_level. ci_grid_step is the resolution of the bound.

boot <- TRexScreeningSelector$new(X, y, seed = 3, verbose = FALSE,
                                  screen_control = trex_screen_control(use_bootstrap_CI = TRUE,
                                                                       R_boot = 200),
                                  control = trex_control(K = 20))
boot$select()
boot$used_bootstrap
#> [1] TRUE
boot$confidence_level
#> [1] 0.389
boot$estimated_FDR
#> [1] 0.1428571
boot$selected_indices
#> [1]   3  25  50 120 200 231 280

Dependent variables

trex_method switches the screening to dependency-aware dummies: "TREX_DA_AR1", "TREX_DA_EQUI" and "TREX_DA_BLOCK_EQUI" draw autoregressive, equi-correlated and block equi-correlated dummies. estimated_correlation reports the correlation the dummies were drawn with.

ar <- TRexScreeningSelector$new(X, y, seed = 3, verbose = FALSE,
                                screen_control = trex_screen_control(trex_method = "TREX_DA_AR1"),
                                control = trex_control(K = 20))
ar$select()
round(ar$estimated_correlation, 3)
#> [1] 0.008

The biobank screening

TRexBiobankScreeningSelector screens one or many phenotypes against one genotype matrix. Y is a vector for one phenotype or a matrix with one column per phenotype. The acceptance window of the estimated FDR is trex_biobank_control(lower_bound_FDR, upper_bound_FDR), and tFDR is the target of the fallback selector.

Y <- cbind(y, as.numeric(X %*% (0.5 * beta) + rnorm(n)), rnorm(n))
colnames(Y) <- c("strong", "weak", "null")
bio <- TRexBiobankScreeningSelector$new(X, Y, tFDR = 0.1, seed = 3, verbose = FALSE,
                                        biobank_control = trex_biobank_control(0.05, 0.15),
                                        control = trex_control(K = 20))
res <- bio$select()
length(res)
#> [1] 3

Each record names the phenotype, the method that produced the result and whether the fallback ran.

summary_tab <- do.call(rbind, lapply(res, function(r)
  data.frame(phenotype = r$phenotype_index, selected = length(r$selected_indices),
             estimated_FDR = round(r$estimated_FDR, 3), method = r$method_used,
             fallback = r$used_fallback_trex)))
summary_tab
#>   phenotype selected estimated_FDR                     method fallback
#> 1         1        7         0.143 Screen-TRex (bootstrap-CI)    FALSE
#> 2         2        3         0.100           T-Rex (fallback)     TRUE
#> 3         3        0         0.100           T-Rex (fallback)     TRUE
res[[1]]$selected_indices
#> [1]   3  25  50 120 200 231 280

The record also keeps the ordinary and the bootstrap screening results (selected_indices_screen_ordinary, estimated_FDR_screen_bootstrap, …), so the decision can be audited. For a single phenotype select() returns one record instead of a list.

Large designs

Both screening selectors accept a memory-mapped design (mem_matrix()), which the memory-mapping page describes. The combination of the screening and a mapped design is what solved sparse regression problems with five million variables on a laptop [2].

Examples

The TRexSelector_Examples repository holds runnable demos and Monte Carlo studies of Screen-T-Rex and the biobank screening, in memory and memory-mapped, with their results and a README per demo, in R, Python and C++.

Bibliography

[1]

J. Machkour, M. Muma, and D. P. Palomar, “False discovery rate control for fast screening of large-scale genomics biobanks,” in Proc. 22nd IEEE statist. Signal process. Workshop (SSP), 2023, pp. 666–670. doi: 10.1109/SSP53291.2023.10207957.

[2]

F. Scheidt, J. Machkour, and M. Muma, “Solving FDR-controlled sparse regression problems with five million variables on a laptop,” in Proc. IEEE 9th int. Workshop comput. Adv. Multi-sensor adapt. Process. (CAMSAP), 2023, pp. 116–120. doi: 10.1109/CAMSAP58249.2023.10403478.

Screening: Screen-T-Rex and the biobank screening