Last updated: 2025-02-11
Checks: 6 1
Knit directory: zinck-website/
This reproducible R Markdown analysis was created with workflowr (version 1.7.1). The Checks tab describes the reproducibility checks that were applied when the results were created. The Past versions tab lists the development history.
The R Markdown file has unstaged changes. To know which version of
the R Markdown file created these results, you’ll want to first commit
it to the Git repo. If you’re still working on the analysis, you can
ignore this warning. When you’re finished, you can run
wflow_publish to commit the R Markdown file and build the
HTML.
Great job! The global environment was empty. Objects defined in the global environment can affect the analysis in your R Markdown file in unknown ways. For reproduciblity it’s best to always run the code in an empty environment.
The command set.seed(20240617) was run prior to running
the code in the R Markdown file. Setting a seed ensures that any results
that rely on randomness, e.g. subsampling or permutations, are
reproducible.
Great job! Recording the operating system, R version, and package versions is critical for reproducibility.
Nice! There were no cached chunks for this analysis, so you can be confident that you successfully produced the results during this run.
Great job! Using relative paths to the files within your workflowr project makes it easier to run your code on other machines.
Great! You are using Git for version control. Tracking code development and connecting the code version to the results is critical for reproducibility.
The results in this page were generated with repository version d4297f5. See the Past versions tab to see a history of the changes made to the R Markdown and HTML files.
Note that you need to be careful to ensure that all relevant files for
the analysis have been committed to Git prior to generating the results
(you can use wflow_publish or
wflow_git_commit). workflowr only checks the R Markdown
file, but you know if there are other scripts or data files that it
depends on. Below is the status of the Git repository when the results
were generated:
Ignored files:
Ignored: .DS_Store
Ignored: analysis/.DS_Store
Unstaged changes:
Modified: .gitignore
Modified: analysis/CRC.Rmd
Deleted: analysis/CRC.html
Modified: analysis/Heatmaps.Rmd
Modified: analysis/IBD.Rmd
Modified: analysis/_site.yml
Modified: analysis/index.Rmd
Modified: analysis/simulation.Rmd
Note that any generated files, e.g. HTML, png, CSS, etc., are not included in this status report because it is ok for generated content to have uncommitted changes.
These are the previous versions of the repository in which changes were
made to the R Markdown (analysis/simulation.Rmd) and HTML
(docs/simulation.html) files. If you’ve configured a remote
Git repository (see ?wflow_git_remote), click on the
hyperlinks in the table below to view the files as they were in that
past version.
| File | Version | Author | Date | Message |
|---|---|---|---|---|
| html | d4297f5 | ghoshstats | 2024-12-12 | Manually deploy website with existing HTML files |
| Rmd | d39920f | Patron | 2024-06-18 | Changed the simulation settings |
| html | d39920f | Patron | 2024-06-18 | Changed the simulation settings |
| html | ab6400d | Patron | 2024-06-18 | Build and publish the website |
| Rmd | a6c38f8 | Patron | 2024-06-18 | Add home, experiment, and simulation pages |
We take the colon cancer (CRC) dataset which encompasses five distinct studies each originating from different countries. We calculated the average proportions of species across all \(574\) samples from five metagenomic studies and retained the top \(p = 200\) most abundant species. Then, we randomly sampled \(D = 500\) samples from the cohort of \(574\) observations and normalized their raw OTU counts row-wise, ensuring each row sums to \(1\). To avoid zero probabilities during multinomial sampling, a small constant (\(0.5\)) was added to the raw OTU counts before normalization.
Next, we paired the calculated true proportions \(\Pi\) with sequencing depths resampled from the original dataset. Each sequencing depth was independently resampled and paired with a proportion vector from a different subject. Using these paired proportions and sequencing depths, we generated the taxa count matrix \(\mathbf{X}^{D \times p}\) via multinomial sampling. For each sample \(d = 1, 2, \ldots, 500\), the count vector corresponding to the \(d^{\text{th}}\) row of \(\mathbf{X}\) was generated as: \[ \mathbf{X}_d \sim \text{Multinomial}(N_d, \Pi_d), \] where \(N_d\) is the sequencing depth for sample \(d\), and \(\Pi_d\) is the vector of true proportions for sample \(d\) with \(\Pi_{dj}\) (j = 1…p) as its element.
Given the simulated count, among the first \(200\) taxa, \(30\) were randomly assumed to be associated with the outcome, while the remaining \(p-30\) taxa had no effect. To simulate outcomes reflecting realistic microbiome data, signal strengths \(\boldsymbol{s}\) for the associated \(30\) biomarkers were defined as: \[ s_j \overset{d}{=} \begin{cases} (2Z_j - 1) \frac{U_1}{\sqrt{C_j}}, & \text{for continuous outcomes}, \\ (2Z_j - 1) \frac{U_2}{\sqrt{C_j}}, & \text{for binary outcomes}, \end{cases} \] where \(U_1 \sim \text{Uniform}(6, 12)\) and \(U_2 \sim \text{Uniform}(30, 50)\) represent the signal magnitudes, \(Z_j \sim \text{Bernoulli}(0.5)\) determines the sign of \(s_j\) (\(+1\) or \(-1\)), and \(C_j\) is the column-wise average abundance of the \(j^{\text{th}}\) bio-marker in \(\Pi\).
Outcomes were generated using the following model: \[ g(\mathbb{E}(y_d)) = \sum_{j=1}^{p} \left \{ \Pi_{dj}^2 \times \left(\frac{s_j}{2}\right) + \Pi_{dj} \times s_j \right \} \]
where \(g(\cdot)\) is the identity link for continuous outcomes and the logit link for binary outcomes. We replicate a single iteration of this setting for both continuous and binary outcomes.
generate_data <- function(p, seed) {
dcount <- count[, order(decreasing = TRUE, colSums(count, na.rm = TRUE), apply(count, 2L, paste, collapse = ''))] # Ordering the columns with decreasing abundance
# Randomly sampling patients from 574 observations
set.seed(seed)
norm_count <- count / rowSums(count)
col_means <- colMeans(norm_count > 0)
indices <- which(col_means > 0.2)
sorted_indices <- indices[order(col_means[indices], decreasing = TRUE)]
if (p %in% c(200, 300, 400)) {
dcount <- count[, sorted_indices][, 1:p]
sel_index <- sort(sample(1:nrow(dcount), 500))
dcount <- dcount[sel_index, ]
original_OTU <- dcount + 0.5
seq_depths <- rowSums(original_OTU)
Pi <- sweep(original_OTU, 1, seq_depths, "/")
n <- nrow(Pi)
col_abundances <- colMeans(Pi)
##### Generating continuous responses ######
set.seed(1)
signal_indices <- sample(1:min(p, 200), 30, replace = FALSE) # Randomly selecting 30 indices for signal injection
signals <- (2 * rbinom(30, 1, 0.5) - 1) * runif(30, 6, 12)
kBeta <- numeric(p)
kBeta[signal_indices] <- signals / sqrt(col_abundances[signal_indices])
eps <- rnorm(n, mean = 0, sd = 1)
Y <- Pi^2 %*% (kBeta / 2) + Pi %*% kBeta + eps
##### Generating binary responses #####
set.seed(1)
signals <- (2 * rbinom(30, 1, 0.5) - 1) * runif(30, 30, 50)
kBeta[signal_indices] <- signals / sqrt(col_abundances[signal_indices])
pr <- 1 / (1 + exp(-(Pi^2 %*% (kBeta / 2) + Pi %*% kBeta)))
Y_bin <- rbinom(n, 1, pr)
######### Generate a copy of X #########
X <- matrix(0, nrow = nrow(Pi), ncol = ncol(Pi))
nSeq <- seq_depths
# Loop over each row to generate the new counts based on the multinomial distribution
set.seed(1)
for (i in 1:nrow(Pi)) {
X[i, ] <- rmultinom(1, size = nSeq[i], prob = Pi[i, ])
}
} else {
print("Enter p within 200 to 400")
}
colnames(X) <- colnames(Pi)
return(list(Y = Y, X = X, Y_bin = Y_bin, signal_indices = signal_indices))
}
ntaxa = 200 # Change to p = 200, 300, 400 accordingly.
X <- generate_data(p=ntaxa, seed=1)$X
Y1 <- generate_data(p=ntaxa, seed=1)$YOn extracting the sample taxa count matrix \(X\) and continuous outcomes \(Y_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}\).
# K =16, seed=19
fit <- fit.zinck(X,num_clusters = 18,method="ADVI",seed=1,boundary_correction = TRUE, prior_ZIGD = TRUE)
theta <- fit$theta
beta <- fit$beta
X_tilde <- zinck::generateKnockoff(X,theta,beta,seed=1) ## getting the knockoff copyNow that we have generated the knockoff copy, we will fit a model associating the response \(Y_1\) with the augmented set of covariates \([X,\tilde{X}]^{D \times 2p}\). We fit a Random Forest model and compute the importance statistics for each feature. Then, we detect the non-zero features for the FDR threshold of \(0.2\). Then, we compute the Power (True Positive Rate) and the empirical FDR by comparing the selected set of taxa with the true set of taxa.
index <- generate_data(p=ntaxa,seed=1)$signal_indices
################### Continuous Outcomes #########################
############## Target FDR = 0.2 ###############
cts_rf <- suppressWarnings(zinck.filter(X,X_tilde,Y1,model="Random Forest",fdr=0.2,ntrees=5000,offset=1,mtry=400,seed=15,metric = "Accuracy", rftuning = TRUE))
index_est <- cts_rf[["selected"]]
### 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_zinck_cts <- FP / (FP + TP) # Evaluating the empirical False Discovery Rate
estimated_power_zinck_cts <- TP / (TP + FN) # Evaluating the empirical Power or TPR
print(paste("Estimated FDR for Target 0.2 (Continuous case):", estimated_FDR_zinck_cts))[1] "Estimated FDR for Target 0.2 (Continuous case): 0.166666666666667"print(paste("Estimated Power for Target 0.2 (Continuous case):", estimated_power_zinck_cts))[1] "Estimated Power for Target 0.2 (Continuous case): 0.5"It is to be noted that zinck is outcome agnostic! That
is, it performs feature selection irrespective of the outcome type.
Thus, we demonstrate its power and empirical FDR for binary outcomes
\(Y_2\).
Y2 <- generate_data(p=ntaxa, seed=1)$Y_bin
index <- generate_data(p=ntaxa,seed=1)$signal_indices
################### Binary Outcomes #########################
############## Target FDR = 0.2 ###############
bin_rf <- suppressWarnings(zinck.filter(X,X_tilde,as.factor(Y2),model="Random Forest",fdr=0.2,offset=0,mtry=45,seed=68,metric="Gini",rftuning = TRUE))
index_est <- bin_rf[["selected"]]
### 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_zinck_bin <- FP / (FP + TP) # Evaluating the empirical False Discovery Rate
estimated_power_zinck_bin <- TP / (TP + FN) # Evaluating the empirical Power or TPR
print(paste("Estimated FDR for Target 0.2 (Binary case):", estimated_FDR_zinck_bin))[1] "Estimated FDR for Target 0.2 (Binary case): 0.181818181818182"print(paste("Estimated Power for Target 0.2 (Binary case):", estimated_power_zinck_bin))[1] "Estimated Power for Target 0.2 (Binary case): 0.3"We now compare the performance of zinck with two
standard knockoff filters namely, Model-X Knockoff Filter (MX-KF) and
the standard LDA based knockoff filter (LDA-KF) keeping the same FDR
threshold of \(0.2\). It is to be noted
that for fitting MX-KF, we need to generate the second-order knockoff
copies for the log-normalized version of \(\mathbf{X}\).
################ Continuous Outcomes (MX-KF) ##################
set.seed(1)
Xlog = log_normalize(X)
Xlog_tilde = create.second_order(Xlog)
index_est <- zinck.filter(Xlog,Xlog_tilde,Y1,model="Random Forest",fdr=0.2,seed=4)$selected
### 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_KF_cts <- FP / (FP + TP) # Evaluating the empirical False Discovery Rate
estimated_power_KF_cts <- TP / (TP + FN) # Evaluating the empirical Power or TPR
cat("MX-KF (Continuous) metrics\n")MX-KF (Continuous) metricscat("Estimated FDR for Target 0.2:", estimated_FDR_KF_cts, "\n")Estimated FDR for Target 0.2: 0.1428571 cat("Estimated Power for Target 0.2:", estimated_power_KF_cts, "\n")Estimated Power for Target 0.2: 0.2 ############### Binary Outcomes (MX-KF) ###################
index_est <- zinck.filter(Xlog,Xlog_tilde,as.factor(Y2),model="Random Forest",fdr=0.2,seed=1)$selected
### 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_KF_bin <- FP / (FP + TP) # Evaluating the empirical False Discovery Rate
estimated_power_KF_bin <- TP / (TP + FN) # Evaluating the empirical Power or TPR
cat("MX-KF (Binary) metrics\n")MX-KF (Binary) metricscat("Estimated FDR for Target 0.2:", estimated_FDR_KF_bin, "\n")Estimated FDR for Target 0.2: 0.4545455 cat("Estimated Power for Target 0.2:", estimated_power_KF_bin, "\n")Estimated Power for Target 0.2: 0.2 ################ Continuous Outcomes (LDA-KF) ##################
df.LDA = as(as.matrix(X),"dgCMatrix")
set.seed(1)
vanilla.LDA <- LDA(df.LDA,k=8,method="VEM")
theta.LDA <- vanilla.LDA@gamma
beta.LDA <- vanilla.LDA@beta
beta.LDA <- t(apply(beta.LDA,1,function(row) row/sum(row)))
set.seed(1)
X_tilde.LDA <- zinck::generateKnockoff(X,theta.LDA,beta.LDA,seed=1) ## Generating vanilla LDA knockoff copy
index_est <- zinck.filter(X,X_tilde.LDA,Y1,model="Random Forest",fdr=0.2,offset=0,seed=2)$selected
### 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_LDA_cts <- FP / (FP + TP) # Evaluating the empirical False Discovery Rate
estimated_power_LDA_cts <- TP / (TP + FN) # Evaluating the empirical Power or TPR
cat("LDA-KF (Continuous) metrics\n")LDA-KF (Continuous) metricscat("Estimated FDR for Target 0.2:", estimated_FDR_LDA_cts, "\n")Estimated FDR for Target 0.2: 0.44 cat("Estimated Power for Target 0.2:", estimated_power_LDA_cts, "\n")Estimated Power for Target 0.2: 0.4666667 ############### Binary Outcomes (LDA-KF) ###################
index_est <- zinck.filter(X,X_tilde.LDA,as.factor(Y2),model="Random Forest",fdr=0.2,offset=1,seed=4,rftuning = TRUE,metric="Gini",mtry=67)$selected
### 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_LDA_bin <- FP / (FP + TP) # Evaluating the empirical False Discovery Rate
estimated_power_LDA_bin <- TP / (TP + FN) # Evaluating the empirical Power or TPR
cat("LDA-KF (Binary) metrics\n")LDA-KF (Binary) metricscat("Estimated FDR for Target 0.2:", estimated_FDR_LDA_bin, "\n")Estimated FDR for Target 0.2: 0.3076923 cat("Estimated Power for Target 0.2:", estimated_power_LDA_bin, "\n")Estimated Power for Target 0.2: 0.3 To summarize the metrics for the three methods, we create the following table.
# Example collected results for demonstration purposes
results <- data.frame(
Method = c("Zinck", "MX-KF", "LDA-KF", "Zinck", "MX-KF", "LDA-KF"),
Outcome = c("Continuous", "Continuous", "Continuous", "Binary", "Binary", "Binary"),
Power = c(estimated_power_zinck_cts, estimated_power_KF_cts, estimated_power_LDA_cts, estimated_power_zinck_bin, estimated_power_KF_bin, estimated_power_LDA_bin),
FDR = c(estimated_FDR_zinck_cts, estimated_FDR_KF_cts, estimated_FDR_LDA_cts, estimated_FDR_zinck_bin, estimated_FDR_KF_bin, estimated_FDR_LDA_bin)
)
kable(results, format = "html", caption = "Power and Empirical FDR Summary") %>%
kable_styling(bootstrap_options = c("striped", "hover"))| Method | Outcome | Power | FDR |
|---|---|---|---|
| Zinck | Continuous | 0.5000000 | 0.1666667 |
| MX-KF | Continuous | 0.2000000 | 0.1428571 |
| LDA-KF | Continuous | 0.4666667 | 0.4400000 |
| Zinck | Binary | 0.3000000 | 0.1818182 |
| MX-KF | Binary | 0.2000000 | 0.4545455 |
| LDA-KF | Binary | 0.3000000 | 0.3076923 |
This shows that zinck has decent power and controlled
FDR for both binary and continuous outcome scenarios compared to the
other knockoff filters.
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] kableExtra_1.4.0 knitr_1.43 topicmodels_0.2-14
[4] glmnet_4.1-7 Matrix_1.5-1 dplyr_1.1.2
[7] kosel_0.0.1 phyloseq_1.38.0 rstan_2.21.8
[10] StanHeaders_2.21.0-7 ggplot2_3.4.2 knockoff_0.3.6
[13] reshape2_1.4.4 zinck_0.0.0.9000
loaded via a namespace (and not attached):
[1] colorspace_2.1-0 modeltools_0.2-23 rprojroot_2.0.3
[4] XVector_0.34.0 fs_1.6.2 rstudioapi_0.14
[7] RSpectra_0.16-1 fansi_1.0.4 ranger_0.15.1
[10] xml2_1.3.4 codetools_0.2-19 splines_4.1.3
[13] cachem_1.0.8 ade4_1.7-22 jsonlite_1.8.5
[16] workflowr_1.7.1 cluster_2.1.4 compiler_4.1.3
[19] fastmap_1.1.1 cli_3.6.1 later_1.3.1
[22] htmltools_0.5.5 prettyunits_1.1.1 tools_4.1.3
[25] igraph_1.4.2 NLP_0.2-1 gtable_0.3.3
[28] glue_1.6.2 GenomeInfoDbData_1.2.7 Rcpp_1.0.10
[31] slam_0.1-50 Biobase_2.54.0 jquerylib_0.1.4
[34] vctrs_0.6.5 Biostrings_2.62.0 rhdf5filters_1.6.0
[37] multtest_2.50.0 ape_5.7-1 svglite_2.1.1
[40] nlme_3.1-162 iterators_1.0.14 ordinalNet_2.12
[43] xfun_0.39 stringr_1.5.0 ps_1.7.5
[46] lifecycle_1.0.3 zlibbioc_1.40.0 MASS_7.3-60
[49] scales_1.2.1 promises_1.2.0.1 parallel_4.1.3
[52] biomformat_1.22.0 rhdf5_2.38.1 inline_0.3.19
[55] yaml_2.3.7 gridExtra_2.3 loo_2.6.0
[58] sass_0.4.6 stringi_1.7.12 highr_0.10
[61] S4Vectors_0.32.4 foreach_1.5.2 randomForest_4.7-1.1
[64] permute_0.9-7 BiocGenerics_0.40.0 pkgbuild_1.4.2
[67] shape_1.4.6 GenomeInfoDb_1.30.1 rlang_1.1.1
[70] pkgconfig_2.0.3 systemfonts_1.0.4 matrixStats_0.63.0
[73] bitops_1.0-7 evaluate_0.21 lattice_0.21-8
[76] Rhdf5lib_1.16.0 Rdsdp_1.0.5.2 processx_3.8.1
[79] tidyselect_1.2.0 plyr_1.8.8 magrittr_2.0.3
[82] R6_2.5.1 IRanges_2.28.0 generics_0.1.3
[85] DBI_1.1.3 pillar_1.9.0 whisker_0.4.1
[88] withr_2.5.0 mgcv_1.8-42 survival_3.5-5
[91] RCurl_1.98-1.12 tibble_3.2.1 crayon_1.5.2
[94] utf8_1.2.3 rmarkdown_2.22 grid_4.1.3
[97] data.table_1.14.8 callr_3.7.3 git2r_0.32.0
[100] vegan_2.6-4 digest_0.6.31 tm_0.7-8
[103] httpuv_1.6.11 RcppParallel_5.1.7 stats4_4.1.3
[106] munsell_0.5.0 viridisLite_0.4.2 bslib_0.5.0