From cb9403bc8fe304c6f75c517eaf3cf51a49b95ecc Mon Sep 17 00:00:00 2001 From: elior rahmani Date: Mon, 20 May 2019 12:21:52 -0700 Subject: [PATCH] adding more refactor tests --- manual.pdf | Bin 168013 -> 168013 bytes .../refactor.R_est.adjusted.C.remove.txt | 200 +++++++++ tests/assets/refactor.R_est.adjusted.txt | 200 +++++++++ tests/assets/refactor.R_est.txt | 400 +++++++++--------- .../assets/refactor.ranked_list.adjusted.txt | 1 + tests/assets/refactor.ranked_list.txt | 2 +- tests/testthat/test_refactor.R | 69 +-- tests/testthat/test_tca.R | 382 ++++++++--------- tests/testthat/test_tcasub.R | 26 +- tests/testthat/test_test_data.R | 44 +- 10 files changed, 870 insertions(+), 454 deletions(-) create mode 100644 tests/assets/refactor.R_est.adjusted.C.remove.txt create mode 100644 tests/assets/refactor.R_est.adjusted.txt create mode 100644 tests/assets/refactor.ranked_list.adjusted.txt diff --git a/manual.pdf b/manual.pdf index cc206f17f3e4c3d3bee963231c48a04d5f95700a..6f9dd5b9fb6d22eb54badb964695cad11c416d2a 100644 GIT binary patch delta 135 zcmX@Rfa~l6u7(!IEll@1H4Tjn%uS4S4b0UI4AeEb^nLSFToOxC6*OF|j0_Ad4NQR~ zSjl$wE+#K#XD3&4S4(pPLsu6=GgBv1b2k$=M{{EnAji_k#nQ>mPQiwdlI>G^nRYS( E0ITjI;xSF{-8@RbSI=i?!8#_5!xH`K!I$9c-I~f~Tn%F7W5K^*zN-xt+ FCIHkkBf0

0.99, TRUE) -# } -# -# }) +test_that("Compare the result of refactor with those of the matlab version of refactor", { + + skip_on_cran() + + basedir <- "../assets/" + + # load data and the matlab results + X <- as.matrix(read.table(paste(basedir,"X.txt",sep=""), header = FALSE, sep=",")) + rownames(X) <- 1:nrow(X) + C <- as.matrix(read.table(paste(basedir,"C2.txt",sep=""), header = FALSE, sep=",")) + # add intercept + C <- cbind(matrix(1,nrow(X),1),C) + t <- 50 + k <- 3 + matlab.ref_comp <- as.matrix(read.table(paste(basedir,"refactor.R_est.txt",sep=""), header = FALSE, sep=",")) + matlab.ranked_list <- as.matrix(read.table(paste(basedir,"refactor.ranked_list.txt",sep=""), header = FALSE, sep=",")) + matlab.ref_comp.adj <- as.matrix(read.table(paste(basedir,"refactor.R_est.adjusted.txt",sep=""), header = FALSE, sep=",")) + matlab.ref_comp.adj.C.remove <- as.matrix(read.table(paste(basedir,"refactor.R_est.adjusted.C.remove.txt",sep=""), header = FALSE, sep=",")) + matlab.ranked_list.adj <- as.matrix(read.table(paste(basedir,"refactor.ranked_list.adjusted.txt",sep=""), header = FALSE, sep=",")) + + # run refactor + ref1 <- refactor(t(X), k, sparsity = t, C = NULL, sd_threshold = 0, num_comp = ncol(matlab.ref_comp), debug = TRUE, log_file = NULL) + ref2 <- refactor(t(X), k, sparsity = t, C = C, C.remove = FALSE, sd_threshold = 0, num_comp = ncol(matlab.ref_comp), debug = TRUE, log_file = NULL) + ref3 <- refactor(t(X), k, sparsity = t, C = C, C.remove = TRUE, sd_threshold = 0, num_comp = ncol(matlab.ref_comp), debug = TRUE, log_file = NULL) + + # evaluating the feature ranking + expect_equal(sum(as.numeric(substring(ref1[["ranked_list"]], 2)) == matlab.ranked_list), ncol(X)) + expect_equal(sum(as.numeric(substring(ref2[["ranked_list"]], 2)) == matlab.ranked_list.adj), ncol(X)) + expect_equal(sum(as.numeric(substring(ref3[["ranked_list"]], 2)) == matlab.ranked_list.adj), ncol(X)) + + # evaluate correlation of the refactor components + for (h in 1:ncol(matlab.ref_comp)){ + expect_equal(abs(cor(matlab.ref_comp[,h],ref1[["scores"]][,h])) > 0.99, TRUE) + + expect_equal(abs(cor(matlab.ref_comp.adj.C.remove[,h],ref3[["scores"]][,h])) > 0.99, TRUE) + } + # the version with C.remove=TRUE has more differences compared with the matlab version - evaluate the the first two components only. + for (h in 1:2){ + expect_equal(abs(cor(matlab.ref_comp.adj[,h],ref2[["scores"]][,h])) > 0.99, TRUE) + } + +}) diff --git a/tests/testthat/test_tca.R b/tests/testthat/test_tca.R index cf2fee9..456c42e 100644 --- a/tests/testthat/test_tca.R +++ b/tests/testthat/test_tca.R @@ -8,194 +8,194 @@ context("Test TCA") # expect_is(data$W, "matrix") # }) -# -# test_that("Comparison with the results given by the matlab version", { -# -# skip_on_cran() -# -# basedir <- "../assets/" -# # load data -# X <- as.matrix(read.table(paste(basedir,"X.txt",sep=""), header = FALSE, sep=",")) -# W <- as.matrix(read.table(paste(basedir,"W.txt",sep=""), header = FALSE, sep=",")) -# C1 <- as.matrix(read.table(paste(basedir,"C1.txt",sep=""), header = FALSE, sep=",")) -# C2 <- as.matrix(read.table(paste(basedir,"C2.txt",sep=""), header = FALSE, sep=",")) -# -# # load real parameters -# gammas <- as.matrix(read.table(paste(basedir,"gammas.txt",sep=""), header = FALSE, sep=",")) -# deltas <- as.matrix(read.table(paste(basedir,"deltas.txt",sep=""), header = FALSE, sep=",")) -# mus <- as.matrix(read.table(paste(basedir,"mus.txt",sep=""), header = FALSE, sep=",")) -# sigmas <- as.matrix(read.table(paste(basedir,"sigmas.txt",sep=""), header = FALSE, sep=",")) -# -# # load matlab estimates -# matlab.mus_hat <- as.matrix(read.table(paste(basedir,"mus_hat.txt",sep=""), header = FALSE, sep=",")) -# matlab.sigmas_hat <- as.matrix(read.table(paste(basedir,"sigmas_hat.txt",sep=""), header = FALSE, sep=",")) -# matlab.gammas_hat <- as.matrix(read.table(paste(basedir,"gammas_hat.txt",sep=""), header = FALSE, sep=",")) -# matlab.deltas_hat <- as.matrix(read.table(paste(basedir,"deltas_hat.txt",sep=""), header = FALSE, sep=",")) -# matlab.W_hat <- as.matrix(read.table(paste(basedir,"W_hat.txt",sep=""), header = FALSE, sep=",")) -# -# # run tca -# res <- tca(t(X), W, C1 = C1, C2 = C2, refit_W = TRUE, refit_W.sparsity = ncol(X), parallel = FALSE, log_file = NULL) -# -# # compare parameter estimates -# for (h in 1:ncol(res$mus_hat)){ -# expect_equal(cor(res$mus_hat[,h],matlab.mus_hat[,h]) > 0.99, TRUE) -# expect_equal(cor(res$sigmas_hat[,h],matlab.sigmas_hat[,h]) > 0.95, TRUE) -# } -# for (h in 1:ncol(res$gammas_hat)){ -# expect_equal(cor(res$gammas_hat[,h],matlab.gammas_hat[,h]) > 0.95, TRUE) -# } -# for (h in 1:ncol(res$deltas_hat)){ -# expect_equal(cor(res$deltas_hat[,h],matlab.deltas_hat[,h]) > 0.99, TRUE) -# } -# for (h in 1:ncol(res$W)){ -# expect_equal(cor(res$W[,h],matlab.W_hat[,h]) > 0.99, TRUE) -# } -# -# # use noisy initial estimates of W -# W_noisy <- as.matrix(read.table(paste(basedir,"W_noisy.txt",sep=""), header = FALSE, sep=",")) -# res2 <- tca(t(X), W_noisy, C1 = C1, C2 = C2, refit_W = TRUE, refit_W.sparsity = ncol(X), parallel = FALSE) -# matlab.W_noisy_hat <- as.matrix(read.table(paste(basedir,"W_noisy_hat.txt",sep=""), header = FALSE, sep=",")) -# for (h in 1:ncol(res2$W)){ -# expect_equal(cor(res2$W[,h],matlab.W_noisy_hat[,h]) > 0.99, TRUE) -# } -# -# }) -# -# -# test_that("Comparison with the results given by the matlab version for tensor", { -# -# skip_on_cran() -# -# basedir <- "../assets/" -# # load data -# X <- as.matrix(read.table(paste(basedir,"X.txt",sep=""), header = FALSE, sep=",")) -# tca.mdl <- list() -# tca.mdl[["W"]] <- as.matrix(read.table(paste(basedir,"W.txt",sep=""), header = FALSE, sep=",")) -# tca.mdl[["mus_hat"]] <- t(as.matrix(read.table(paste(basedir,"mus.txt",sep=""), header = FALSE, sep=","))) -# tca.mdl[["sigmas_hat"]] <- t(as.matrix(read.table(paste(basedir,"sigmas.txt",sep=""), header = FALSE, sep=","))) -# tca.mdl[["gammas_hat"]] <- as.matrix(read.table(paste(basedir,"gammas.txt",sep=""), header = FALSE, sep=",")) -# tca.mdl[["deltas_hat"]] <- as.matrix(read.table(paste(basedir,"deltas.txt",sep=""), header = FALSE, sep=",")) -# tca.mdl[["tau_hat"]] <- 0.01 -# tca.mdl[["C1"]] <- as.matrix(read.table(paste(basedir,"C1.txt",sep=""), header = FALSE, sep=",")) -# tca.mdl[["C2"]] <- as.matrix(read.table(paste(basedir,"C2.txt",sep=""), header = FALSE, sep=",")) -# -# # load matlab's estimates of Z -# matlab.Z1_hat <- t(as.matrix(read.table(paste(basedir,"Z1_hat.txt",sep=""), header = FALSE, sep=","))) -# matlab.Z2_hat <- t(as.matrix(read.table(paste(basedir,"Z2_hat.txt",sep=""), header = FALSE, sep=","))) -# matlab.Z3_hat <- t(as.matrix(read.table(paste(basedir,"Z3_hat.txt",sep=""), header = FALSE, sep=","))) -# -# # run tensor -# Z_hat <- tensor(t(X), tca.mdl, log_file = NULL) -# -# # compare Z estimates -# th = 0.001 -# expect_equal(mean((Z_hat[[1]]-matlab.Z1_hat)**2) < th, TRUE) -# expect_equal(mean((Z_hat[[2]]-matlab.Z2_hat)**2) < th, TRUE) -# expect_equal(mean((Z_hat[[3]]-matlab.Z3_hat)**2) < th, TRUE) -# for (h in 1:ncol(matlab.Z1_hat)){ -# expect_equal(cor(matlab.Z1_hat[,h],Z_hat[[1]][,h]) > 0.99, TRUE) -# } -# for (h in 1:ncol(matlab.Z1_hat)){ -# expect_equal(cor(matlab.Z2_hat[,h],Z_hat[[2]][,h]) > 0.99, TRUE) -# } -# for (h in 1:ncol(matlab.Z1_hat)){ -# expect_equal(cor(matlab.Z3_hat[,h],Z_hat[[3]][,h]) > 0.99, TRUE) -# } -# -# }) -# -# -# -# run_ewas <- function(X, W, tca.mdl, basedir, yfile, test){ -# Y <- as.matrix(read.table(paste(basedir,yfile,sep=""), header = FALSE, sep=",")) -# m <- ncol(Y) -# k <- ncol(W) -# tca.mdl2 <- tca.mdl -# pvals <- matrix(0,m,1) -# if (test == "marginal" | test == "marginal_conditional" ) pvals <- matrix(0,m,k) -# for (j in 1:m){ -# tca.mdl2[["mus_hat"]] <- t(as.matrix(tca.mdl[["mus_hat"]][j,])) -# tca.mdl2[["sigmas_hat"]] <- t(as.matrix(tca.mdl[["sigmas_hat"]][j,])) -# tca.mdl2[["deltas_hat"]] <- t(as.matrix(tca.mdl[["deltas_hat"]][j,])) -# tca.mdl2[["gammas_hat"]] <- t(as.matrix(tca.mdl[["gammas_hat"]][j,])) -# if (test == "custom"){ -# res <- tcareg(X = t(X[,j]), tca.mdl = tca.mdl2, y = as.matrix(Y[,j]), test = test, save_results = FALSE, null_model = c("V3"), alternative_model = c("V1","V2","V3"), log_file = NULL) -# }else{ -# res <- tcareg(X = t(X[,j]), tca.mdl = tca.mdl2, y = as.matrix(Y[,j]), test = test, save_results = FALSE, log_file = NULL) -# } -# if (test == "single_effect" | test == "custom" | test == "joint"){ -# pvals[j,] <- res$pvals[1] -# } -# if (test == "marginal" | test == "marginal_conditional"){ -# for (h in 1:k){ -# pvals[j,h] <- res[[h]]$pvals[1] -# } -# } -# } -# return(pvals) -# } -# -# -# test_that("Evaluate tcareg with the results of the matlab version", { -# -# skip_on_cran() -# -# basedir <- "../assets/ewas/" -# # load data -# X <- as.matrix(read.table(paste(basedir,"TCA_methylation_exp2_sim1_k_6_n_500_data.txt",sep=""), header = FALSE, sep=",")) -# colnames(X) <- 1:ncol(X) -# rownames(X) <- 1:nrow(X) -# W <- as.matrix(read.table(paste(basedir,"TCA_methylation_exp2_sim1_k_6_n_500_W.txt",sep=""), header = FALSE, sep=",")) -# m <- ncol(X) -# k <- ncol(W) -# -# tca.mdl <- tca(t(X), W, refit_W = FALSE, parallel = FALSE, log_file = NULL) -# -# # (1) joint test -# # power -# pvals <- run_ewas(X, W, tca.mdl, basedir, "TCA_methylation_exp2_sim1_k_6_n_500_Y_joint_effect_size_10_W_noise_level_0.txt", test="joint") -# expect_equal(sum(pvals < 0.05/m)/m > 0.5, TRUE) -# # fps rate -# pvals <- run_ewas(X, W, tca.mdl, basedir, "TCA_methylation_exp2_sim1_k_6_n_500_Y_fps_effect_size_1_W_noise_level_0.txt", test="joint") -# expect_equal(sum(pvals < 0.05/m)/m < 0.05, TRUE) -# -# # (2) marginal test -# # power -# pvals <- run_ewas(X, W, tca.mdl, basedir, "TCA_methylation_exp2_sim1_k_6_n_500_Y_marginal_ct_1_effect_size_5_W_noise_level_0.txt", test="marginal") -# expect_equal(sum(pvals[,1] < 0.05/m)/m > 0.9, TRUE) -# # fps rate -# pvals <- run_ewas(X, W, tca.mdl, basedir, "TCA_methylation_exp2_sim1_k_6_n_500_Y_fps_effect_size_1_W_noise_level_0.txt", test="marginal") -# for (h in 1:k){ -# expect_equal(sum(pvals[,h] < 0.05)/m < 0.1, TRUE) -# } -# -# # (3) conditional marginal test -# pvals <- run_ewas(X, W, tca.mdl, basedir, "TCA_methylation_exp2_sim1_k_6_n_500_Y_marginal_ct_1_effect_size_5_W_noise_level_0.txt", test="marginal_conditional") -# expect_equal(sum(pvals[,1] < 0.05)/m > 0.9, TRUE) -# for (h in 2:k){ -# expect_equal(sum(pvals[,h] < 0.05)/m > 0.1, FALSE) -# } -# pvals <- run_ewas(X, W, tca.mdl, basedir, "TCA_methylation_exp2_sim1_k_6_n_500_Y_marginal_ct_3_effect_size_10_W_noise_level_0.txt", test="marginal_conditional") -# l <- 3 -# expect_equal(sum(pvals[,l] < 0.05)/m > 0.1, TRUE) -# for (h in setdiff(1:k,l)){ -# expect_equal(sum(pvals[,h] < 0.05)/m > 0.1, FALSE) -# } -# -# # (4) joint test - single beta -# # power -# pvals <- run_ewas(X, W, tca.mdl, basedir, "TCA_methylation_exp2_sim1_k_6_n_500_Y_joint_single_beta_effect_size_10_W_noise_level_0.txt", test="single_effect") -# expect_equal(sum(pvals < 0.05/m)/m > 0.8, TRUE) -# # fps rate -# pvals <- run_ewas(X, W, tca.mdl, basedir, "TCA_methylation_exp2_sim1_k_6_n_500_Y_fps_effect_size_1_W_noise_level_0.txt", test="single_effect") -# expect_equal(sum(pvals < 0.05/m)/m > 0.05, FALSE) -# -# # (5) custom test -# # power -# pvals <- run_ewas(X, W, tca.mdl, basedir, "TCA_methylation_exp2_sim1_k_6_n_500_Y_joint_effect_size_10_W_noise_level_0.txt", test="custom") -# expect_equal(sum(pvals < 0.05/m)/m > 0.4, TRUE) -# # fps rate -# pvals <- run_ewas(X, W, tca.mdl, basedir, "TCA_methylation_exp2_sim1_k_6_n_500_Y_fps_effect_size_1_W_noise_level_0.txt", test="custom") -# expect_equal(sum(pvals < 0.05/m)/m > 0.1, FALSE) -# -# }) + +test_that("Comparison with the results given by the matlab version", { + + skip_on_cran() + + basedir <- "../assets/" + # load data + X <- as.matrix(read.table(paste(basedir,"X.txt",sep=""), header = FALSE, sep=",")) + W <- as.matrix(read.table(paste(basedir,"W.txt",sep=""), header = FALSE, sep=",")) + C1 <- as.matrix(read.table(paste(basedir,"C1.txt",sep=""), header = FALSE, sep=",")) + C2 <- as.matrix(read.table(paste(basedir,"C2.txt",sep=""), header = FALSE, sep=",")) + + # load real parameters + gammas <- as.matrix(read.table(paste(basedir,"gammas.txt",sep=""), header = FALSE, sep=",")) + deltas <- as.matrix(read.table(paste(basedir,"deltas.txt",sep=""), header = FALSE, sep=",")) + mus <- as.matrix(read.table(paste(basedir,"mus.txt",sep=""), header = FALSE, sep=",")) + sigmas <- as.matrix(read.table(paste(basedir,"sigmas.txt",sep=""), header = FALSE, sep=",")) + + # load matlab estimates + matlab.mus_hat <- as.matrix(read.table(paste(basedir,"mus_hat.txt",sep=""), header = FALSE, sep=",")) + matlab.sigmas_hat <- as.matrix(read.table(paste(basedir,"sigmas_hat.txt",sep=""), header = FALSE, sep=",")) + matlab.gammas_hat <- as.matrix(read.table(paste(basedir,"gammas_hat.txt",sep=""), header = FALSE, sep=",")) + matlab.deltas_hat <- as.matrix(read.table(paste(basedir,"deltas_hat.txt",sep=""), header = FALSE, sep=",")) + matlab.W_hat <- as.matrix(read.table(paste(basedir,"W_hat.txt",sep=""), header = FALSE, sep=",")) + + # run tca + res <- tca(t(X), W, C1 = C1, C2 = C2, refit_W = TRUE, refit_W.sparsity = ncol(X), parallel = FALSE, log_file = NULL) + + # compare parameter estimates + for (h in 1:ncol(res$mus_hat)){ + expect_equal(cor(res$mus_hat[,h],matlab.mus_hat[,h]) > 0.99, TRUE) + expect_equal(cor(res$sigmas_hat[,h],matlab.sigmas_hat[,h]) > 0.95, TRUE) + } + for (h in 1:ncol(res$gammas_hat)){ + expect_equal(cor(res$gammas_hat[,h],matlab.gammas_hat[,h]) > 0.95, TRUE) + } + for (h in 1:ncol(res$deltas_hat)){ + expect_equal(cor(res$deltas_hat[,h],matlab.deltas_hat[,h]) > 0.99, TRUE) + } + for (h in 1:ncol(res$W)){ + expect_equal(cor(res$W[,h],matlab.W_hat[,h]) > 0.99, TRUE) + } + + # use noisy initial estimates of W + W_noisy <- as.matrix(read.table(paste(basedir,"W_noisy.txt",sep=""), header = FALSE, sep=",")) + res2 <- tca(t(X), W_noisy, C1 = C1, C2 = C2, refit_W = TRUE, refit_W.sparsity = ncol(X), parallel = FALSE, log_file = NULL) + matlab.W_noisy_hat <- as.matrix(read.table(paste(basedir,"W_noisy_hat.txt",sep=""), header = FALSE, sep=",")) + for (h in 1:ncol(res2$W)){ + expect_equal(cor(res2$W[,h],matlab.W_noisy_hat[,h]) > 0.99, TRUE) + } + +}) + + +test_that("Comparison with the results given by the matlab version for tensor", { + + skip_on_cran() + + basedir <- "../assets/" + # load data + X <- as.matrix(read.table(paste(basedir,"X.txt",sep=""), header = FALSE, sep=",")) + tca.mdl <- list() + tca.mdl[["W"]] <- as.matrix(read.table(paste(basedir,"W.txt",sep=""), header = FALSE, sep=",")) + tca.mdl[["mus_hat"]] <- t(as.matrix(read.table(paste(basedir,"mus.txt",sep=""), header = FALSE, sep=","))) + tca.mdl[["sigmas_hat"]] <- t(as.matrix(read.table(paste(basedir,"sigmas.txt",sep=""), header = FALSE, sep=","))) + tca.mdl[["gammas_hat"]] <- as.matrix(read.table(paste(basedir,"gammas.txt",sep=""), header = FALSE, sep=",")) + tca.mdl[["deltas_hat"]] <- as.matrix(read.table(paste(basedir,"deltas.txt",sep=""), header = FALSE, sep=",")) + tca.mdl[["tau_hat"]] <- 0.01 + tca.mdl[["C1"]] <- as.matrix(read.table(paste(basedir,"C1.txt",sep=""), header = FALSE, sep=",")) + tca.mdl[["C2"]] <- as.matrix(read.table(paste(basedir,"C2.txt",sep=""), header = FALSE, sep=",")) + + # load matlab's estimates of Z + matlab.Z1_hat <- t(as.matrix(read.table(paste(basedir,"Z1_hat.txt",sep=""), header = FALSE, sep=","))) + matlab.Z2_hat <- t(as.matrix(read.table(paste(basedir,"Z2_hat.txt",sep=""), header = FALSE, sep=","))) + matlab.Z3_hat <- t(as.matrix(read.table(paste(basedir,"Z3_hat.txt",sep=""), header = FALSE, sep=","))) + + # run tensor + Z_hat <- tensor(t(X), tca.mdl, log_file = NULL) + + # compare Z estimates + th = 0.001 + expect_equal(mean((Z_hat[[1]]-matlab.Z1_hat)**2) < th, TRUE) + expect_equal(mean((Z_hat[[2]]-matlab.Z2_hat)**2) < th, TRUE) + expect_equal(mean((Z_hat[[3]]-matlab.Z3_hat)**2) < th, TRUE) + for (h in 1:ncol(matlab.Z1_hat)){ + expect_equal(cor(matlab.Z1_hat[,h],Z_hat[[1]][,h]) > 0.99, TRUE) + } + for (h in 1:ncol(matlab.Z1_hat)){ + expect_equal(cor(matlab.Z2_hat[,h],Z_hat[[2]][,h]) > 0.99, TRUE) + } + for (h in 1:ncol(matlab.Z1_hat)){ + expect_equal(cor(matlab.Z3_hat[,h],Z_hat[[3]][,h]) > 0.99, TRUE) + } + +}) + + + +run_ewas <- function(X, W, tca.mdl, basedir, yfile, test){ + Y <- as.matrix(read.table(paste(basedir,yfile,sep=""), header = FALSE, sep=",")) + m <- ncol(Y) + k <- ncol(W) + tca.mdl2 <- tca.mdl + pvals <- matrix(0,m,1) + if (test == "marginal" | test == "marginal_conditional" ) pvals <- matrix(0,m,k) + for (j in 1:m){ + tca.mdl2[["mus_hat"]] <- t(as.matrix(tca.mdl[["mus_hat"]][j,])) + tca.mdl2[["sigmas_hat"]] <- t(as.matrix(tca.mdl[["sigmas_hat"]][j,])) + tca.mdl2[["deltas_hat"]] <- t(as.matrix(tca.mdl[["deltas_hat"]][j,])) + tca.mdl2[["gammas_hat"]] <- t(as.matrix(tca.mdl[["gammas_hat"]][j,])) + if (test == "custom"){ + res <- tcareg(X = t(X[,j]), tca.mdl = tca.mdl2, y = as.matrix(Y[,j]), test = test, save_results = FALSE, null_model = c("V3"), alternative_model = c("V1","V2","V3"), log_file = NULL) + }else{ + res <- tcareg(X = t(X[,j]), tca.mdl = tca.mdl2, y = as.matrix(Y[,j]), test = test, save_results = FALSE, log_file = NULL) + } + if (test == "single_effect" | test == "custom" | test == "joint"){ + pvals[j,] <- res$pvals[1] + } + if (test == "marginal" | test == "marginal_conditional"){ + for (h in 1:k){ + pvals[j,h] <- res[[h]]$pvals[1] + } + } + } + return(pvals) +} + + +test_that("Evaluate tcareg with the results of the matlab version", { + + skip_on_cran() + + basedir <- "../assets/ewas/" + # load data + X <- as.matrix(read.table(paste(basedir,"TCA_methylation_exp2_sim1_k_6_n_500_data.txt",sep=""), header = FALSE, sep=",")) + colnames(X) <- 1:ncol(X) + rownames(X) <- 1:nrow(X) + W <- as.matrix(read.table(paste(basedir,"TCA_methylation_exp2_sim1_k_6_n_500_W.txt",sep=""), header = FALSE, sep=",")) + m <- ncol(X) + k <- ncol(W) + + tca.mdl <- tca(t(X), W, refit_W = FALSE, parallel = FALSE, log_file = NULL) + + # (1) joint test + # power + pvals <- run_ewas(X, W, tca.mdl, basedir, "TCA_methylation_exp2_sim1_k_6_n_500_Y_joint_effect_size_10_W_noise_level_0.txt", test="joint") + expect_equal(sum(pvals < 0.05/m)/m > 0.5, TRUE) + # fps rate + pvals <- run_ewas(X, W, tca.mdl, basedir, "TCA_methylation_exp2_sim1_k_6_n_500_Y_fps_effect_size_1_W_noise_level_0.txt", test="joint") + expect_equal(sum(pvals < 0.05/m)/m < 0.05, TRUE) + + # (2) marginal test + # power + pvals <- run_ewas(X, W, tca.mdl, basedir, "TCA_methylation_exp2_sim1_k_6_n_500_Y_marginal_ct_1_effect_size_5_W_noise_level_0.txt", test="marginal") + expect_equal(sum(pvals[,1] < 0.05/m)/m > 0.9, TRUE) + # fps rate + pvals <- run_ewas(X, W, tca.mdl, basedir, "TCA_methylation_exp2_sim1_k_6_n_500_Y_fps_effect_size_1_W_noise_level_0.txt", test="marginal") + for (h in 1:k){ + expect_equal(sum(pvals[,h] < 0.05)/m < 0.1, TRUE) + } + + # (3) conditional marginal test + pvals <- run_ewas(X, W, tca.mdl, basedir, "TCA_methylation_exp2_sim1_k_6_n_500_Y_marginal_ct_1_effect_size_5_W_noise_level_0.txt", test="marginal_conditional") + expect_equal(sum(pvals[,1] < 0.05)/m > 0.9, TRUE) + for (h in 2:k){ + expect_equal(sum(pvals[,h] < 0.05)/m > 0.1, FALSE) + } + pvals <- run_ewas(X, W, tca.mdl, basedir, "TCA_methylation_exp2_sim1_k_6_n_500_Y_marginal_ct_3_effect_size_10_W_noise_level_0.txt", test="marginal_conditional") + l <- 3 + expect_equal(sum(pvals[,l] < 0.05)/m > 0.1, TRUE) + for (h in setdiff(1:k,l)){ + expect_equal(sum(pvals[,h] < 0.05)/m > 0.1, FALSE) + } + + # (4) joint test - single beta + # power + pvals <- run_ewas(X, W, tca.mdl, basedir, "TCA_methylation_exp2_sim1_k_6_n_500_Y_joint_single_beta_effect_size_10_W_noise_level_0.txt", test="single_effect") + expect_equal(sum(pvals < 0.05/m)/m > 0.8, TRUE) + # fps rate + pvals <- run_ewas(X, W, tca.mdl, basedir, "TCA_methylation_exp2_sim1_k_6_n_500_Y_fps_effect_size_1_W_noise_level_0.txt", test="single_effect") + expect_equal(sum(pvals < 0.05/m)/m > 0.05, FALSE) + + # (5) custom test + # power + pvals <- run_ewas(X, W, tca.mdl, basedir, "TCA_methylation_exp2_sim1_k_6_n_500_Y_joint_effect_size_10_W_noise_level_0.txt", test="custom") + expect_equal(sum(pvals < 0.05/m)/m > 0.4, TRUE) + # fps rate + pvals <- run_ewas(X, W, tca.mdl, basedir, "TCA_methylation_exp2_sim1_k_6_n_500_Y_fps_effect_size_1_W_noise_level_0.txt", test="custom") + expect_equal(sum(pvals < 0.05/m)/m > 0.1, FALSE) + +}) diff --git a/tests/testthat/test_tcasub.R b/tests/testthat/test_tcasub.R index 0a48474..b291785 100644 --- a/tests/testthat/test_tcasub.R +++ b/tests/testthat/test_tcasub.R @@ -2,16 +2,16 @@ library("TCA") context("Test tcasub") -# test_that("subset the results of tca", { -# n = 20; m = 30; k = 3; tau = 0.01; p1 = 0; p2 = 0; -# res <- test_data(n, m, k, p1, p2, tau, log_file = NULL) -# X <- res$X -# W <- res$W -# mus <- res$mus -# sigmas <- res$sigmas -# tca.mdl <- tca(X, W, refit_W = FALSE, parallel = FALSE, log_file = NULL) -# remove <- c("1","3","5") -# features <- setdiff(as.character(1:m),remove) -# tca.mdl.subset <- tcasub(tca.mdl, features, log_file = NULL) -# expect_equal(sum(rownames(tca.mdl.subset$mus_hat) == features) == length(features), TRUE) -# }) +test_that("subset the results of tca", { + n = 20; m = 30; k = 3; tau = 0.01; p1 = 0; p2 = 0; + res <- test_data(n, m, k, p1, p2, tau, log_file = NULL) + X <- res$X + W <- res$W + mus <- res$mus + sigmas <- res$sigmas + tca.mdl <- tca(X, W, refit_W = FALSE, parallel = FALSE, log_file = NULL) + remove <- c("1","3","5") + features <- setdiff(as.character(1:m),remove) + tca.mdl.subset <- tcasub(tca.mdl, features, log_file = NULL) + expect_equal(sum(rownames(tca.mdl.subset$mus_hat) == features) == length(features), TRUE) +}) diff --git a/tests/testthat/test_test_data.R b/tests/testthat/test_test_data.R index 2456a4a..b4ae1fa 100644 --- a/tests/testthat/test_test_data.R +++ b/tests/testthat/test_test_data.R @@ -1,25 +1,25 @@ library("TCA") context("Test test_data") -# -# test_that("generate test data", { -# -# n = 200; m = 100; k = 3; tau = 0.01; p1 = 2; p2 = 2; -# res <- test_data(n, m, k, p1, p2, tau) -# -# X <- res$X -# W <- res$W -# mus <- res$mus -# sigmas <- res$sigmas -# C1 <- res$C1 -# C2 <- res$C2 -# gammas <- res$gammas -# deltas <- res$deltas -# -# tca.mdl <- tca(X, W, C1 = C1, C2 = C2, refit_W = TRUE, refit_W.sparsity = nrow(X), parallel = FALSE, log_file = NULL) -# Z_hat <- tensor(X, tca.mdl) -# expect_is(Z_hat[[1]], "matrix") -# expect_is(Z_hat[[2]], "matrix") -# expect_is(Z_hat[[3]], "matrix") -# -# }) + +test_that("generate test data", { + + n = 200; m = 100; k = 3; tau = 0.01; p1 = 2; p2 = 2; + res <- test_data(n, m, k, p1, p2, tau, log_file = NULL) + + X <- res$X + W <- res$W + mus <- res$mus + sigmas <- res$sigmas + C1 <- res$C1 + C2 <- res$C2 + gammas <- res$gammas + deltas <- res$deltas + + tca.mdl <- tca(X, W, C1 = C1, C2 = C2, refit_W = TRUE, refit_W.sparsity = nrow(X), parallel = FALSE, log_file = NULL) + Z_hat <- tensor(X, tca.mdl, log_file = NULL) + expect_is(Z_hat[[1]], "matrix") + expect_is(Z_hat[[2]], "matrix") + expect_is(Z_hat[[3]], "matrix") + +})