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 FDR estimate 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

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

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].