Simulations

Empirical Power and False Discovery Rate for different target thresholds

We take the colon cancer (CRC) dataset which encompasses five distinct studies each originating from different countries. We randomly select roughly half of the samples from each metagenomic study, which are pooled together for further analysis. For the sake of this example, we retain the top 100 most abundant species according to the average proportions. Among these species, we randomly selected 30 biomarkers and simulated both continuous and binary outcomes using the symmetric reformulation of the log-contrast model:

$$g(E(y_i))=\sum_{j=1}^{p} \log(\bar{X}_{ij})\beta_j \quad \text{subject to} \quad \sum_{j=1}^{p} \beta_j=0$$ where g(.) is the link function and $\bar{X}_{ij}$ is the observed proportions based on the count matrix with zero replaced by 0.5 to avoid taking log on zero values. For simulating the continuous responses, the link function g(.) is the identity function, and for simulating the binary responses, the link function g(.) is the logit link function. We adopt the signal strengths in such a way that the zero-sum constraint is satisfied. We replicate a single iteration of this setting for both continuous and binary outcomes.

load("/Users/Patron/Documents/zinLDA research/count.Rdata")
load("/Users/Patron/Documents/zinLDA research/meta.RData")
library(rstan)

Loading required package: StanHeaders

Loading required package: ggplot2

rstan (Version 2.21.8, GitRev: 2e1f913d3ca3)

For execution on a local, multicore CPU with excess RAM we recommend calling
options(mc.cores = parallel::detectCores()).
To avoid recompilation of unchanged Stan programs, we recommend calling
rstan_options(auto_write = TRUE)

library(glmnet)

Loading required package: Matrix

Loaded glmnet 4.1-7

library(knockoff)
library(zinck)
library(kosel)
library(dplyr)

Attaching package: 'dplyr'

The following objects are masked from 'package:stats':

    filter, lag

The following objects are masked from 'package:base':

    intersect, setdiff, setequal, union

library(reshape2)

generate_data <- function(p,seed){
  dcount <- count[,order(decreasing=T,colSums(count,na.rm=T),
                         apply(count,2L,paste,collapse=''))] ## ordering the columns w/ decreasing abundance
  ####### Randomly sampling patients from 574 observations #######
  set.seed(seed)
  sel_index <- rbinom(nrow(meta),size=1,prob=0.5)
  selected_samples <- which(sel_index==1)
  meta_selected <- meta[selected_samples,]
  X <- dcount[selected_samples,]
  if(p == 100){
    X1 <- X[,1:p]
    n = nrow(X1)
    
    Five1 = c(-3,3,2.5,-1, -1.5)
    Five2 = c(3,3,-2,-2,-2)
    Five3 = c(1,-1,3,-2,-1)
    Five4 = c(-1,1,2,-1,-1)
    Five5 = c(3,3,-3,-2,-1)
    Five6 = c(-1,1,-2,1,1)
    Five_all <- c(Five1,Five2,Five3,Five4,Five5,Five6)  ## Signals satisfy sum-to-zero constraint ##
    randBeta <- rep(0,p)
    set.seed(1)
    rand_indices <- sample(1:p,size=30,replace=FALSE)  ## Randomly injecting 30 signals out of p=100 ##
    set.seed(1)
    randBeta[rand_indices] <- sample(Five_all,size=30,replace=FALSE)
    
    W1 <- log_normalize(X1) ## log-normalizing counts to make it compositional
    
    ##### Generating continuous responses ######
    set.seed(1)
    eps=rnorm(n,mean = 0, sd=1)
    Y1 <- W1 %*% randBeta + eps
    
    ##### Generating binary responses #####
    set.seed(1)
    pr = 1/(1+exp(-W1 %*% randBeta))
    Y1_bin = rbinom(n,1,pr)
  } else if (p %in% c(200,300,400)) {
    X1 <- X[,1:p]
    n = nrow(X1)
    
    Five1 = c(-3,3,2.5,-1, -1.5)
    Five2 = c(3,3,-2,-2,-2)
    Five3 = c(1,-1,3,-2,-1)
    Five4 = c(-1,1,2,-1,-1)
    Five5 = c(3,3,-3,-2,-1)
    Five6 = c(-1,1,-2,1,1)
    Five_all <- c(Five1,Five2,Five3,Five4,Five5,Five6)  ## Signals satisfy sum-to-zero constraint ##
    randBeta <- rep(0,p)
    set.seed(1)
    rand_indices <- sample(1:200,size=30,replace=FALSE)  ## Randomly injecting 30 signals out of p=200 ##
    set.seed(1)
    randBeta[rand_indices] <- sample(Five_all,size=30,replace=FALSE)
    
    W1 <- log_normalize(X1) ## log-normalizing counts to make it compositional
    
    ##### Generating continuous responses ######
    set.seed(1)
    eps=rnorm(n,mean = 0, sd=1)
    Y1 <- W1 %*% randBeta + eps
    
    ##### Generating binary responses #####
    set.seed(1)
    pr = 1/(1+exp(-W1 %*% randBeta))
    Y1_bin = rbinom(n,1,pr)
  }

  return(list(Y1 = Y1, X1 = X1,  W1 = W1, Y1_bin = Y1_bin, index = rand_indices))
}

Now, let’s look at the case when $p = 100$. On extracting the sample taxa count matrix $X_1$, we train the zinck model using ADVI. Plugging in the learnt parameters into the generative framework of zinck, we generate the knockoff matrix $\tilde{X}_1$.

X1 <- generate_data(p=100, seed=2)$X1
fit1 <- zinck::fit.zinck(X1,num_clusters = 6, method="ADVI", seed = 1)
theta <- fit1$theta
beta <- fit1$beta
X1_tilde <- generateKnockoff(X1,theta,beta,seed=1)

Now that we have generated the knockoff copy, we will fit a model associating the response $Y_1$ (both continuous and binary) with the augmented set of covariates $[X_1,\tilde{X}_1]^{D \times 2p}$. We fit a glmnet model separately for each outcome type and compute the importance statistics for each feature. Then, we detect the non-zero features for four different FDR thresholds 0.05, 0.1, 0.15, and 0.2.

W1 <- generate_data(p=100, seed=2)$W1 ## Log-Normalized Sample Taxa matrix 
W_tilde1 <- log_normalize(X1_tilde) ## Log-Normalized Knockoff Matrix
Y1 <- generate_data(p=100, seed=2)$Y1 ## Continuous Responses
Y1_bin <- generate_data(p = 100, seed = 2)$Y1_bin ## Binary Responses
index <- generate_data(p = 100, seed = 2)$index ## True signals

################### Continuous Outcomes #########################
############## FDR = 0.05 ###############
index_est <- zinck.filter(W1,W_tilde1,Y1,model="glmnet",fdr=0.05)

Selected variables:
 [1]  1  7 14 20 21 28 33 34 35 37 38 39 42 43 44 51 54 59 68 70 73 74 79 82 83
[26] 84 85 87 96 97

# ### Evaluation metrics ###
FN <- sum(index %in% index_est == FALSE) ## False Negatives
FP <- sum(index_est %in% index == FALSE) ## False Positives
TP <- sum(index_est %in% index == TRUE) ## True Positives

estimated_FDR <- FP / (FP + TP) # Evaluating the empirical False Discovery Rate
estimated_power <- TP / (TP + FN) # Evaluating the empirical Power or TPR

print(paste("Estimated FDR for Target 0.05 (Continuous case):", estimated_FDR))

[1] "Estimated FDR for Target 0.05 (Continuous case): 0"

print(paste("Estimated Power for Target 0.05 (Continuous case):", estimated_power))

[1] "Estimated Power for Target 0.05 (Continuous case): 1"

############## FDR = 0.1 ###############
index_est <- zinck.filter(W1,W_tilde1,Y1,model="glmnet",fdr=0.1)

Selected variables:
 [1]  1  7 14 20 21 28 33 34 35 37 38 39 42 43 44 51 54 59 68 70 73 74 79 82 83
[26] 84 85 87 96 97

# ### Evaluation metrics ###
FN <- sum(index %in% index_est == FALSE) ## False Negatives
FP <- sum(index_est %in% index == FALSE) ## False Positives
TP <- sum(index_est %in% index == TRUE) ## True Positives

estimated_FDR <- FP / (FP + TP) # Evaluating the empirical False Discovery Rate
estimated_power <- TP / (TP + FN) # Evaluating the empirical Power or TPR

print(paste("Estimated FDR for Target 0.1 (Continuous case):", estimated_FDR))

[1] "Estimated FDR for Target 0.1 (Continuous case): 0"

print(paste("Estimated Power for Target 0.1 (Continuous case):", estimated_power))

[1] "Estimated Power for Target 0.1 (Continuous case): 1"

############## FDR = 0.15 ###############
index_est <- zinck.filter(W1,W_tilde1,Y1,model="glmnet",fdr=0.15)

Selected variables:
 [1]  1  7 14 20 21 28 33 34 35 37 38 39 42 43 44 51 54 59 68 70 73 74 79 82 83
[26] 84 85 87 96 97

# ### Evaluation metrics ###
FN <- sum(index %in% index_est == FALSE) ## False Negatives
FP <- sum(index_est %in% index == FALSE) ## False Positives
TP <- sum(index_est %in% index == TRUE) ## True Positives

estimated_FDR <- FP / (FP + TP) # Evaluating the empirical False Discovery Rate
estimated_power <- TP / (TP + FN) # Evaluating the empirical Power or TPR

print(paste("Estimated FDR for Target 0.15 (Continuous case):", estimated_FDR))

[1] "Estimated FDR for Target 0.15 (Continuous case): 0"

print(paste("Estimated Power for Target 0.15 (Continuous case):", estimated_power))

[1] "Estimated Power for Target 0.15 (Continuous case): 1"

############## FDR = 0.2 ###############
index_est <- zinck.filter(W1,W_tilde1,Y1,model="glmnet",fdr=0.2)

Selected variables:
 [1]  1  7 14 20 21 28 33 34 35 37 38 39 42 43 44 51 54 59 68 70 73 74 79 82 83
[26] 84 85 87 96 97

# ### Evaluation metrics ###
FN <- sum(index %in% index_est == FALSE) ## False Negatives
FP <- sum(index_est %in% index == FALSE) ## False Positives
TP <- sum(index_est %in% index == TRUE) ## True Positives

estimated_FDR <- FP / (FP + TP) # Evaluating the empirical False Discovery Rate
estimated_power <- TP / (TP + FN) # Evaluating the empirical Power or TPR

print(paste("Estimated FDR for Target 0.2 (Continuous case):", estimated_FDR))

[1] "Estimated FDR for Target 0.2 (Continuous case): 0"

print(paste("Estimated Power for Target 0.2 (Continuous case):", estimated_power))

[1] "Estimated Power for Target 0.2 (Continuous case): 1"

################### Binary Outcomes #########################
############## FDR = 0.05 ###############
index_est <- zinck.filter(W1,W_tilde1,Y1_bin,model="glmnet",fdr=0.05)

Selected variables:
 [1]  7 20 28 33 34 35 37 38 43 44 54 59 68 70 73 74 82 87 97

# ### Evaluation metrics ###
FN <- sum(index %in% index_est == FALSE) ## False Negatives
FP <- sum(index_est %in% index == FALSE) ## False Positives
TP <- sum(index_est %in% index == TRUE) ## True Positives

estimated_FDR <- FP / (FP + TP) # Evaluating the empirical False Discovery Rate
estimated_power <- TP / (TP + FN) # Evaluating the empirical Power or TPR

print(paste("Estimated FDR for Target 0.05 (Binary case):", estimated_FDR))

[1] "Estimated FDR for Target 0.05 (Binary case): 0"

print(paste("Estimated Power for Target 0.05 (Binary case):", estimated_power))

[1] "Estimated Power for Target 0.05 (Binary case): 0.633333333333333"

############## FDR = 0.1 ###############
index_est <- zinck.filter(W1,W_tilde1,Y1_bin,model="glmnet",fdr=0.1)

Selected variables:
 [1]  1  7 20 28 33 34 35 37 38 43 44 54 59 68 70 73 74 82 87 97

# ### Evaluation metrics ###
FN <- sum(index %in% index_est == FALSE) ## False Negatives
FP <- sum(index_est %in% index == FALSE) ## False Positives
TP <- sum(index_est %in% index == TRUE) ## True Positives

estimated_FDR <- FP / (FP + TP) # Evaluating the empirical False Discovery Rate
estimated_power <- TP / (TP + FN) # Evaluating the empirical Power or TPR

print(paste("Estimated FDR for Target 0.1 (Binary case):", estimated_FDR))

[1] "Estimated FDR for Target 0.1 (Binary case): 0"

print(paste("Estimated Power for Target 0.1 (Binary case):", estimated_power))

[1] "Estimated Power for Target 0.1 (Binary case): 0.666666666666667"

############## FDR = 0.15 ###############
index_est <- zinck.filter(W1,W_tilde1,Y1_bin,model="glmnet",fdr=0.15)

Selected variables:
 [1]  1  7 20 21 28 33 34 35 37 38 43 44 54 59 63 68 70 73 74 82 84 87 97

# ### Evaluation metrics ###
FN <- sum(index %in% index_est == FALSE) ## False Negatives
FP <- sum(index_est %in% index == FALSE) ## False Positives
TP <- sum(index_est %in% index == TRUE) ## True Positives

estimated_FDR <- FP / (FP + TP) # Evaluating the empirical False Discovery Rate
estimated_power <- TP / (TP + FN) # Evaluating the empirical Power or TPR

print(paste("Estimated FDR for Target 0.15 (Binary case):", estimated_FDR))

[1] "Estimated FDR for Target 0.15 (Binary case): 0.0434782608695652"

print(paste("Estimated Power for Target 0.15 (Binary case):", estimated_power))

[1] "Estimated Power for Target 0.15 (Binary case): 0.733333333333333"

############## FDR = 0.2 ###############
index_est <- zinck.filter(W1,W_tilde1,Y1_bin,model="glmnet",fdr=0.2)

Selected variables:
 [1]  1  7 20 21 28 33 34 35 37 38 43 44 51 54 59 63 68 70 73 74 82 84 87 97

# ### Evaluation metrics ###
FN <- sum(index %in% index_est == FALSE) ## False Negatives
FP <- sum(index_est %in% index == FALSE) ## False Positives
TP <- sum(index_est %in% index == TRUE) ## True Positives

estimated_FDR <- FP / (FP + TP) # Evaluating the empirical False Discovery Rate
estimated_power <- TP / (TP + FN) # Evaluating the empirical Power or TPR

print(paste("Estimated FDR for Target 0.2 (Binary case):", estimated_FDR))

[1] "Estimated FDR for Target 0.2 (Binary case): 0.0416666666666667"

print(paste("Estimated Power for Target 0.2 (Binary case):", estimated_power))

[1] "Estimated Power for Target 0.2 (Binary case): 0.766666666666667"

The users are encouraged to play with other values of $p$, such as $p = 200, 300, 400$.

Session information

sessionInfo()

R version 4.1.3 (2022-03-10)
Platform: x86_64-apple-darwin17.0 (64-bit)
Running under: macOS Big Sur/Monterey 10.16

Matrix products: default
BLAS:   /Library/Frameworks/R.framework/Versions/4.1/Resources/lib/libRblas.0.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/4.1/Resources/lib/libRlapack.dylib

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] reshape2_1.4.4       dplyr_1.1.2          kosel_0.0.1         
 [4] zinck_0.0.0.9000     knockoff_0.3.6       glmnet_4.1-7        
 [7] Matrix_1.5-1         rstan_2.21.8         ggplot2_3.4.2       
[10] StanHeaders_2.21.0-7 workflowr_1.7.1     

loaded via a namespace (and not attached):
 [1] Rcpp_1.0.10          lattice_0.21-8       prettyunits_1.1.1   
 [4] getPass_0.2-2        ps_1.7.5             rprojroot_2.0.3     
 [7] digest_0.6.31        foreach_1.5.2        utf8_1.2.3          
[10] plyr_1.8.8           R6_2.5.1             stats4_4.1.3        
[13] evaluate_0.21        httr_1.4.6           pillar_1.9.0        
[16] rlang_1.1.1          rstudioapi_0.14      whisker_0.4.1       
[19] callr_3.7.3          jquerylib_0.1.4      ordinalNet_2.12     
[22] rmarkdown_2.22       splines_4.1.3        stringr_1.5.0       
[25] loo_2.6.0            munsell_0.5.0        compiler_4.1.3      
[28] httpuv_1.6.11        xfun_0.39            pkgconfig_2.0.3     
[31] pkgbuild_1.4.2       shape_1.4.6          htmltools_0.5.5     
[34] tidyselect_1.2.0     tibble_3.2.1         gridExtra_2.3       
[37] codetools_0.2-19     matrixStats_0.63.0   randomForest_4.7-1.1
[40] fansi_1.0.4          crayon_1.5.2         withr_2.5.0         
[43] later_1.3.1          grid_4.1.3           jsonlite_1.8.5      
[46] gtable_0.3.3         lifecycle_1.0.3      DBI_1.1.3           
[49] git2r_0.32.0         magrittr_2.0.3       scales_1.2.1        
[52] RcppParallel_5.1.7   cli_3.6.1            stringi_1.7.12      
[55] cachem_1.0.8         fs_1.6.2             promises_1.2.0.1    
[58] doParallel_1.0.17    bslib_0.5.0          generics_0.1.3      
[61] vctrs_0.6.5          iterators_1.0.14     tools_4.1.3         
[64] glue_1.6.2           processx_3.8.1       parallel_4.1.3      
[67] fastmap_1.1.1        survival_3.5-5       yaml_2.3.7          
[70] inline_0.3.19        colorspace_2.1-0     knitr_1.43          
[73] sass_0.4.6