From c2e9bfad0d07661514a1f44bf67bdd243386460e Mon Sep 17 00:00:00 2001 From: suzannejin Date: Thu, 12 Sep 2024 14:32:53 +0200 Subject: [PATCH 1/5] add more tests to theta --- tests/testthat/test-PROPD-theta.R | 71 ++++++++++++++++++++++++++----- 1 file changed, 60 insertions(+), 11 deletions(-) diff --git a/tests/testthat/test-PROPD-theta.R b/tests/testthat/test-PROPD-theta.R index 875fa84..95e3568 100644 --- a/tests/testthat/test-PROPD-theta.R +++ b/tests/testthat/test-PROPD-theta.R @@ -1,30 +1,79 @@ library(testthat) library(propr) -# RDA with 3 groups -x <- iris[,1:4] -y <- iris[,5] -v <- vegan::rda(log(x[,1]/x[,2]) ~ y) -pd <- propd(x, as.character(y)) - test_that("3-group RDA agrees with theta", { + # RDA with 3 groups + x <- iris[,1:4] + y <- iris[,5] + v <- vegan::rda(log(x[,1]/x[,2]) ~ y) + pd <- propd(x, as.character(y)) + expect_equal( 1 - sum(v$CCA$eig) / v$tot.chi, pd@results[1,"theta"] ) }) -# RDA with 2 groups -x <- iris[1:100,1:4] -y <- iris[1:100,5] -v <- vegan::rda(log(x[,1]/x[,2]) ~ y) -pd <- propd(x, as.character(y)) test_that("2-group RDA agrees with theta", { + # RDA with 2 groups + x <- iris[1:100,1:4] + y <- iris[1:100,5] + v <- vegan::rda(log(x[,1]/x[,2]) ~ y) + pd <- propd(x, as.character(y)) + expect_equal( 1 - sum(v$CCA$eig) / v$tot.chi, pd@results[1,"theta"] ) }) + + +# data +keep <- iris$Species %in% c("setosa", "versicolor") +counts <- iris[keep, 1:4] * 10 +group <- ifelse(iris[keep, "Species"] == "setosa", "A", "B") + +# calculate propd +pd <- propd(counts, group, p = 5) +pd_w <- propd(counts, group, p = 5, weighted = TRUE) + +test_that("active theta_e matches calculation using theta_d", { + + n1 <- 50 + n2 <- 50 + + expect_equal( + setActive(pd, what = "theta_e")@results$theta, + 1 - pd@results$theta + pmin((n1-1) * pd@results$lrv1, (n2-1) * pd@results$lrv2) / ((n1+n2-1) * pd@results$lrv) + ) + + # when weighted + groups <- lapply(unique(group), function(g) g == group) + ngrp <- length(unique(group)) + W <- pd_w@weights + ps <- lapply(groups, function(g) propr:::omega(W[g,])) + names(ps) <- paste0("p", 1:ngrp) + p <- propr:::omega(W) + expect_equal( + setActive(pd_w, what = "theta_e")@results$theta, + 1 - pd_w@results$theta + pmin(ps[[1]] * pd_w@results$lrv1, ps[[2]] * pd_w@results$lrv2) / (p * pd_w@results$lrv) + ) +}) + +test_that("active theta_f matches calculation using theta_e", { + + expect_equal( + setActive(pd, what = "theta_f")@results$theta, + 1 - setActive(pd, what = "theta_e")@results$theta + ) + + expect_equal( + setActive(pd_w, what = "theta_f")@results$theta, + 1 - setActive(pd_w, what = "theta_e")@results$theta + ) +}) + + From 7a9cb0ca742f912a48b9b97c9c9d2f365817243a Mon Sep 17 00:00:00 2001 From: Suzanne Jin Date: Fri, 13 Sep 2024 14:49:32 +0200 Subject: [PATCH 2/5] add more test graflex --- src/graflex.cpp | 1 + tests/testthat/test-GRAFLEX.R | 12 ++++++++++++ 2 files changed, 13 insertions(+) diff --git a/src/graflex.cpp b/src/graflex.cpp index 255e6e3..236c2bc 100644 --- a/src/graflex.cpp +++ b/src/graflex.cpp @@ -32,6 +32,7 @@ IntegerVector shuffle_and_get_lower_triangle(IntegerMatrix& mat) { IntegerVector shuffled_triangle(n); IntegerVector index = sample(ncol, ncol, false) - 1; + // TODO check if this is correct int k = 0; for (int i = 0; i < ncol; ++i) { int index_i = index[i]; diff --git a/tests/testthat/test-GRAFLEX.R b/tests/testthat/test-GRAFLEX.R index e06fd55..f3126d3 100644 --- a/tests/testthat/test-GRAFLEX.R +++ b/tests/testthat/test-GRAFLEX.R @@ -83,6 +83,18 @@ runGraflex_old <- function(A, K, p = 500) { # run tests # # ===================== # +test_that("check if get_lower_triangle works correctly", { + + n <- 50 + # Create a square and symmetric matrix + mat <- matrix(sample(0:1, n*n, replace = TRUE), nrow = n, byrow = TRUE) + # Expected lower triangle + expected <- mat[lower.tri(mat)] + # Compute the lower triangle + result <- propr:::get_lower_triangle(mat) + # Check if the result matches the expected lower triangle + expect_equal(result, expected) +}) test_that("check if odds ratio is computed correctly", { From bb95ae6d22b05abe91897ac6d502859752e060ea Mon Sep 17 00:00:00 2001 From: suzannejin Date: Fri, 13 Sep 2024 18:30:04 +0200 Subject: [PATCH 3/5] update graflex --- R/RcppExports.R | 20 +++---- src/RcppExports.cpp | 81 ++++++++++------------------ src/graflex.cpp | 125 ++++++++++++++++++++------------------------ 3 files changed, 90 insertions(+), 136 deletions(-) diff --git a/R/RcppExports.R b/R/RcppExports.R index 60f02dc..0a823c6 100644 --- a/R/RcppExports.R +++ b/R/RcppExports.R @@ -105,16 +105,12 @@ ctzRcpp <- function(X) { .Call(`_propr_ctzRcpp`, X) } -get_lower_triangle <- function(mat) { - .Call(`_propr_get_lower_triangle`, mat) +get_triangle <- function(mat) { + .Call(`_propr_get_triangle`, mat) } -shuffle_and_get_lower_triangle <- function(mat) { - .Call(`_propr_shuffle_and_get_lower_triangle`, mat) -} - -binTab <- function(A, G) { - .Call(`_propr_binTab`, A, G) +get_triangle_from_index <- function(mat, index) { + .Call(`_propr_get_triangle_from_index`, mat, index) } getOR <- function(A, G) { @@ -125,12 +121,8 @@ permuteOR <- function(A, Gstar, p = 100L) { .Call(`_propr_permuteOR`, A, Gstar, p) } -getFDR_over <- function(actual, permuted) { - .Call(`_propr_getFDR_over`, actual, permuted) -} - -getFDR_under <- function(actual, permuted) { - .Call(`_propr_getFDR_under`, actual, permuted) +getFDR <- function(actual, permuted) { + .Call(`_propr_getFDR`, actual, permuted) } graflex <- function(A, G, p = 100L) { diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index 1e4572b..ab8c424 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -317,97 +317,74 @@ BEGIN_RCPP return rcpp_result_gen; END_RCPP } -// get_lower_triangle -IntegerVector get_lower_triangle(IntegerMatrix& mat); -RcppExport SEXP _propr_get_lower_triangle(SEXP matSEXP) { +// get_triangle +IntegerVector get_triangle(const IntegerMatrix& mat); +RcppExport SEXP _propr_get_triangle(SEXP matSEXP) { BEGIN_RCPP Rcpp::RObject rcpp_result_gen; Rcpp::RNGScope rcpp_rngScope_gen; - Rcpp::traits::input_parameter< IntegerMatrix& >::type mat(matSEXP); - rcpp_result_gen = Rcpp::wrap(get_lower_triangle(mat)); + Rcpp::traits::input_parameter< const IntegerMatrix& >::type mat(matSEXP); + rcpp_result_gen = Rcpp::wrap(get_triangle(mat)); return rcpp_result_gen; END_RCPP } -// shuffle_and_get_lower_triangle -IntegerVector shuffle_and_get_lower_triangle(IntegerMatrix& mat); -RcppExport SEXP _propr_shuffle_and_get_lower_triangle(SEXP matSEXP) { +// get_triangle_from_index +IntegerVector get_triangle_from_index(const IntegerMatrix& mat, const IntegerVector& index); +RcppExport SEXP _propr_get_triangle_from_index(SEXP matSEXP, SEXP indexSEXP) { BEGIN_RCPP Rcpp::RObject rcpp_result_gen; Rcpp::RNGScope rcpp_rngScope_gen; - Rcpp::traits::input_parameter< IntegerMatrix& >::type mat(matSEXP); - rcpp_result_gen = Rcpp::wrap(shuffle_and_get_lower_triangle(mat)); - return rcpp_result_gen; -END_RCPP -} -// binTab -NumericMatrix binTab(IntegerVector& A, IntegerVector& G); -RcppExport SEXP _propr_binTab(SEXP ASEXP, SEXP GSEXP) { -BEGIN_RCPP - Rcpp::RObject rcpp_result_gen; - Rcpp::RNGScope rcpp_rngScope_gen; - Rcpp::traits::input_parameter< IntegerVector& >::type A(ASEXP); - Rcpp::traits::input_parameter< IntegerVector& >::type G(GSEXP); - rcpp_result_gen = Rcpp::wrap(binTab(A, G)); + Rcpp::traits::input_parameter< const IntegerMatrix& >::type mat(matSEXP); + Rcpp::traits::input_parameter< const IntegerVector& >::type index(indexSEXP); + rcpp_result_gen = Rcpp::wrap(get_triangle_from_index(mat, index)); return rcpp_result_gen; END_RCPP } // getOR -NumericMatrix getOR(IntegerVector& A, IntegerVector& G); +NumericVector getOR(const IntegerVector& A, const IntegerVector& G); RcppExport SEXP _propr_getOR(SEXP ASEXP, SEXP GSEXP) { BEGIN_RCPP Rcpp::RObject rcpp_result_gen; Rcpp::RNGScope rcpp_rngScope_gen; - Rcpp::traits::input_parameter< IntegerVector& >::type A(ASEXP); - Rcpp::traits::input_parameter< IntegerVector& >::type G(GSEXP); + Rcpp::traits::input_parameter< const IntegerVector& >::type A(ASEXP); + Rcpp::traits::input_parameter< const IntegerVector& >::type G(GSEXP); rcpp_result_gen = Rcpp::wrap(getOR(A, G)); return rcpp_result_gen; END_RCPP } // permuteOR -NumericMatrix permuteOR(IntegerMatrix& A, IntegerVector& Gstar, int p); +NumericMatrix permuteOR(const IntegerMatrix& A, const IntegerVector& Gstar, int p); RcppExport SEXP _propr_permuteOR(SEXP ASEXP, SEXP GstarSEXP, SEXP pSEXP) { BEGIN_RCPP Rcpp::RObject rcpp_result_gen; Rcpp::RNGScope rcpp_rngScope_gen; - Rcpp::traits::input_parameter< IntegerMatrix& >::type A(ASEXP); - Rcpp::traits::input_parameter< IntegerVector& >::type Gstar(GstarSEXP); + Rcpp::traits::input_parameter< const IntegerMatrix& >::type A(ASEXP); + Rcpp::traits::input_parameter< const IntegerVector& >::type Gstar(GstarSEXP); Rcpp::traits::input_parameter< int >::type p(pSEXP); rcpp_result_gen = Rcpp::wrap(permuteOR(A, Gstar, p)); return rcpp_result_gen; END_RCPP } -// getFDR_over -double getFDR_over(double actual, NumericVector permuted); -RcppExport SEXP _propr_getFDR_over(SEXP actualSEXP, SEXP permutedSEXP) { -BEGIN_RCPP - Rcpp::RObject rcpp_result_gen; - Rcpp::RNGScope rcpp_rngScope_gen; - Rcpp::traits::input_parameter< double >::type actual(actualSEXP); - Rcpp::traits::input_parameter< NumericVector >::type permuted(permutedSEXP); - rcpp_result_gen = Rcpp::wrap(getFDR_over(actual, permuted)); - return rcpp_result_gen; -END_RCPP -} -// getFDR_under -double getFDR_under(double actual, NumericVector permuted); -RcppExport SEXP _propr_getFDR_under(SEXP actualSEXP, SEXP permutedSEXP) { +// getFDR +List getFDR(double actual, const NumericVector& permuted); +RcppExport SEXP _propr_getFDR(SEXP actualSEXP, SEXP permutedSEXP) { BEGIN_RCPP Rcpp::RObject rcpp_result_gen; Rcpp::RNGScope rcpp_rngScope_gen; Rcpp::traits::input_parameter< double >::type actual(actualSEXP); - Rcpp::traits::input_parameter< NumericVector >::type permuted(permutedSEXP); - rcpp_result_gen = Rcpp::wrap(getFDR_under(actual, permuted)); + Rcpp::traits::input_parameter< const NumericVector& >::type permuted(permutedSEXP); + rcpp_result_gen = Rcpp::wrap(getFDR(actual, permuted)); return rcpp_result_gen; END_RCPP } // graflex -NumericMatrix graflex(IntegerMatrix& A, IntegerMatrix& G, int p); +NumericVector graflex(const IntegerMatrix& A, const IntegerMatrix& G, int p); RcppExport SEXP _propr_graflex(SEXP ASEXP, SEXP GSEXP, SEXP pSEXP) { BEGIN_RCPP Rcpp::RObject rcpp_result_gen; Rcpp::RNGScope rcpp_rngScope_gen; - Rcpp::traits::input_parameter< IntegerMatrix& >::type A(ASEXP); - Rcpp::traits::input_parameter< IntegerMatrix& >::type G(GSEXP); + Rcpp::traits::input_parameter< const IntegerMatrix& >::type A(ASEXP); + Rcpp::traits::input_parameter< const IntegerMatrix& >::type G(GSEXP); Rcpp::traits::input_parameter< int >::type p(pSEXP); rcpp_result_gen = Rcpp::wrap(graflex(A, G, p)); return rcpp_result_gen; @@ -539,13 +516,11 @@ static const R_CallMethodDef CallEntries[] = { {"_propr_count_less_equal_than", (DL_FUNC) &_propr_count_less_equal_than, 2}, {"_propr_count_greater_equal_than", (DL_FUNC) &_propr_count_greater_equal_than, 2}, {"_propr_ctzRcpp", (DL_FUNC) &_propr_ctzRcpp, 1}, - {"_propr_get_lower_triangle", (DL_FUNC) &_propr_get_lower_triangle, 1}, - {"_propr_shuffle_and_get_lower_triangle", (DL_FUNC) &_propr_shuffle_and_get_lower_triangle, 1}, - {"_propr_binTab", (DL_FUNC) &_propr_binTab, 2}, + {"_propr_get_triangle", (DL_FUNC) &_propr_get_triangle, 1}, + {"_propr_get_triangle_from_index", (DL_FUNC) &_propr_get_triangle_from_index, 2}, {"_propr_getOR", (DL_FUNC) &_propr_getOR, 2}, {"_propr_permuteOR", (DL_FUNC) &_propr_permuteOR, 3}, - {"_propr_getFDR_over", (DL_FUNC) &_propr_getFDR_over, 2}, - {"_propr_getFDR_under", (DL_FUNC) &_propr_getFDR_under, 2}, + {"_propr_getFDR", (DL_FUNC) &_propr_getFDR, 2}, {"_propr_graflex", (DL_FUNC) &_propr_graflex, 3}, {"_propr_lr2vlr", (DL_FUNC) &_propr_lr2vlr, 1}, {"_propr_lr2phi", (DL_FUNC) &_propr_lr2phi, 1}, diff --git a/src/graflex.cpp b/src/graflex.cpp index 236c2bc..04d0c28 100644 --- a/src/graflex.cpp +++ b/src/graflex.cpp @@ -5,8 +5,7 @@ using namespace Rcpp; // Function to extract the triangle of a square and symmetric IntegerMatrix // [[Rcpp::export]] -IntegerVector get_lower_triangle(IntegerMatrix& mat) { - int nrow = mat.nrow(); +IntegerVector get_triangle(const IntegerMatrix& mat) { int ncol = mat.ncol(); int n = ncol * (ncol - 1) / 2; @@ -14,7 +13,7 @@ IntegerVector get_lower_triangle(IntegerMatrix& mat) { int k = 0; for (int j = 0; j < ncol; ++j) { - for (int i = j+1; i < nrow; ++i) { + for (int i = j+1; i < ncol; ++i) { triangle[k++] = mat(i, j); } } @@ -22,70 +21,59 @@ IntegerVector get_lower_triangle(IntegerMatrix& mat) { return triangle; } -// Function to shuffle a square and symmetric IntegerMatrix, and get the lower triangle +// Function to get the triangle of a matrix based on the given indeces // [[Rcpp::export]] -IntegerVector shuffle_and_get_lower_triangle(IntegerMatrix& mat) { - int nrow = mat.nrow(); +IntegerVector get_triangle_from_index(const IntegerMatrix& mat, const IntegerVector& index) { int ncol = mat.ncol(); int n = ncol * (ncol - 1) / 2; - IntegerVector shuffled_triangle(n); - IntegerVector index = sample(ncol, ncol, false) - 1; + IntegerVector triangle(n); - // TODO check if this is correct int k = 0; - for (int i = 0; i < ncol; ++i) { - int index_i = index[i]; - for (int j = 0; j < i; ++j) { - int index_j = index[j]; - shuffled_triangle[k++] = mat(index_i, index_j); + for (int j = 0; j < ncol; ++j) { + for (int i = j+1; i < ncol; ++i) { + triangle[k++] = mat(index[i], index[j]); } } - return shuffled_triangle; + return triangle; } // Function to calculate the contingency table // [[Rcpp::export]] -NumericMatrix binTab(IntegerVector& A, IntegerVector& G) { - NumericMatrix tab(2, 2); +NumericVector getOR(const IntegerVector& A, const IntegerVector& G) { int n = A.size(); - tab(0, 0) = sum((1 - A) * (1 - G)); // not in A and not in G - tab(0, 1) = sum((1 - A) * G); // not in A but in G - tab(1, 0) = sum(A * (1 - G)); // in A but not G - tab(1, 1) = sum(A * G); // in A and G - - return tab; -} - -// Function to calculate the odds ratio, and parse all information into a matrix -// [[Rcpp::export]] -NumericMatrix getOR(IntegerVector& A, IntegerVector& G) { - - NumericMatrix tab = binTab(A, G); - double odds_ratio = (tab(0, 0) * tab(1, 1)) / (tab(0, 1) * tab(1, 0)); - - NumericMatrix or_table(1, 8); - or_table(0, 0) = tab(0, 0); - or_table(0, 1) = tab(0, 1); - or_table(0, 2) = tab(1, 0); - or_table(0, 3) = tab(1, 1); - or_table(0, 4) = odds_ratio; - or_table(0, 5) = log(odds_ratio); + // calculate the contingency table + int a = 0, b = 0, c = 0, d = 0; + for (int i = 0; i < n; ++i) { + if (A[i] == 0) { + if (G[i] == 0) ++a; // not in A and not in G + else ++b; // not in A but in G + } else { + if (G[i] == 0) ++c; // in A but not in G + else ++d; // in A and in G + } + } + + // calculate the odds ratio + double odds_ratio = static_cast(a * d) / (b * c); - return or_table; + return NumericVector::create( + a, b, c, d, odds_ratio, std::log(odds_ratio), R_NaN, R_NaN + ); } // Function to calculate the odds ratio and other relevant info for each permutation // [[Rcpp::export]] -NumericMatrix permuteOR(IntegerMatrix& A, IntegerVector& Gstar, int p = 100) { - +NumericMatrix permuteOR(const IntegerMatrix& A, const IntegerVector& Gstar, int p = 100) { + int n = A.ncol(); NumericMatrix or_table(p, 8); // calculate odds ratio for each permutation for (int i = 0; i < p; ++i) { - IntegerVector Astar = shuffle_and_get_lower_triangle(A); + IntegerVector idx = sample(n, n, false) - 1; + IntegerVector Astar = get_triangle_from_index(A, idx); or_table(i, _) = getOR(Astar, Gstar); } @@ -93,46 +81,45 @@ NumericMatrix permuteOR(IntegerMatrix& A, IntegerVector& Gstar, int p = 100) { } // [[Rcpp::export]] -double getFDR_over(double actual, NumericVector permuted) { - double n = permuted.size(); - double fdr = 0.0; - for (int i = 0; i < n; ++i) { - if (permuted[i] >= actual) { - fdr += 1.0; - } - } - fdr /= n; - return fdr; -} +List getFDR(double actual, const NumericVector& permuted) { + int n = permuted.size(); + int count_over = 0; + int count_under = 0; -// [[Rcpp::export]] -double getFDR_under(double actual, NumericVector permuted) { - double n = permuted.size(); - double fdr = 0.0; + // Count values above and below the actual value for (int i = 0; i < n; ++i) { - if (permuted[i] <= actual) { - fdr += 1.0; - } + double current = permuted[i]; + if (current >= actual) ++count_over; + if (current <= actual) ++count_under; } - fdr /= n; - return fdr; + + // Calculate FDR for both "over" and "under" + double fdr_over = static_cast(count_over) / n; + double fdr_under = static_cast(count_under) / n; + + // Return both FDR values as a named list + return List::create( + Named("over") = fdr_over, + Named("under") = fdr_under + ); } // Function to calculate the odds ratio and FDR, given the adjacency matrix A and the knowledge graph G // [[Rcpp::export]] -NumericMatrix graflex(IntegerMatrix& A, IntegerMatrix& G, int p = 100) { +NumericVector graflex(const IntegerMatrix& A, const IntegerMatrix& G, int p = 100) { // get the actual odds ratio - IntegerVector Astar = get_lower_triangle(A); - IntegerVector Gstar = get_lower_triangle(G); - NumericMatrix actual = getOR(Astar, Gstar); + IntegerVector Astar = get_triangle(A); + IntegerVector Gstar = get_triangle(G); + NumericVector actual = getOR(Astar, Gstar); // get distribution of odds ratios on permuted data NumericMatrix permuted = permuteOR(A, Gstar, p); // calculate the FDR - actual(0, 6) = getFDR_under(actual(0, 4), permuted(_, 4)); - actual(0, 7) = getFDR_over(actual(0, 4), permuted(_, 4)); + List fdr = getFDR(actual(4), permuted(_, 4)); + actual(6) = as(fdr["under"]); + actual(7) = as(fdr["over"]); return actual; } From 51945ea6efb84b4cbfdd089dd0a7b36db28ed449 Mon Sep 17 00:00:00 2001 From: suzannejin Date: Fri, 13 Sep 2024 18:30:18 +0200 Subject: [PATCH 4/5] add more tests graflex --- tests/testthat/test-GRAFLEX.R | 92 ++++++++++++++++++++++++++++++----- 1 file changed, 81 insertions(+), 11 deletions(-) diff --git a/tests/testthat/test-GRAFLEX.R b/tests/testthat/test-GRAFLEX.R index f3126d3..efb85c4 100644 --- a/tests/testthat/test-GRAFLEX.R +++ b/tests/testthat/test-GRAFLEX.R @@ -6,14 +6,21 @@ library(Rcpp) # old code for graflex # # ===================== # +cppFunction(" + IntegerVector sample_idx(int n) { + IntegerVector idx = sample(n, n, false) - 1; + return idx; + } +") + permuteOR_old <- function(A, G, p = 500) { Gstar <- G[lower.tri(G)] res <- lapply(1:p, function(i) { # Shuffle the adjacency matrix # index <- sample(1:ncol(A)) - # A <- A[index, index] - # Astar <- A[lower.tri(A)] - Astar <- propr:::shuffle_and_get_lower_triangle(A) + index <- sample_idx(ncol(A))+1 # use the same random sampling generator + A <- A[index, index] + Astar <- A[lower.tri(A)] getOR_old(Astar, Gstar) }) @@ -83,19 +90,53 @@ runGraflex_old <- function(A, K, p = 500) { # run tests # # ===================== # -test_that("check if get_lower_triangle works correctly", { +test_that("check if get_triangle works correctly", { n <- 50 # Create a square and symmetric matrix mat <- matrix(sample(0:1, n*n, replace = TRUE), nrow = n, byrow = TRUE) # Expected lower triangle expected <- mat[lower.tri(mat)] - # Compute the lower triangle - result <- propr:::get_lower_triangle(mat) + # Compute the triangle + result <- propr:::get_triangle(mat) # Check if the result matches the expected lower triangle expect_equal(result, expected) }) +test_that("check if get_triangle_from_index works correctly", { + n <- 50 + + # Create a square and symmetric matrix + mat <- matrix(sample(0:1, n*n, replace = TRUE), nrow = n, byrow = TRUE) + + # Create a random permutation of the indices + index <- sample(1:n) + + # expected lower triangle + expected <- mat[index, index] + expected <- expected[lower.tri(expected)] + + # Compute the lower triangle + result <- propr:::get_triangle_from_index(mat, index-1) + + # Check if the result matches the expected lower triangle + expect_equal(result, expected) +}) + +test_that("check that getFDR works properly", { + + # create a vector of values between 0 and 1 + x <- runif(100) + actual <- runif(1) + + # compute FDR + res <- propr:::getFDR(actual, x) + + # check + expect_equal(res$under, sum(x <= actual) / length(x)) + expect_equal(res$over, sum(x >= actual) / length(x)) +}) + test_that("check if odds ratio is computed correctly", { # Create two vectors of 1 and 0 @@ -113,11 +154,11 @@ test_that("check if odds ratio is computed correctly", { or_table <- propr:::getOR(Astar, Gstar) # check - expect_equal(a, or_table[1,1]) - expect_equal(b, or_table[1,2]) - expect_equal(c, or_table[1,3]) - expect_equal(d, or_table[1,4]) - expect_equal(expected_odds_ratio, or_table[1,5]) + expect_equal(a, as.integer(or_table[1])) + expect_equal(b, as.integer(or_table[2])) + expect_equal(c, as.integer(or_table[3])) + expect_equal(d, as.integer(or_table[4])) + expect_equal(expected_odds_ratio, or_table[5]) }) test_that("check if runGraflex produce the expected results", { @@ -197,6 +238,35 @@ test_that("check reproducibility seed works", { expect_false(isTRUE(all.equal(res1, res3))) }) +test_that("check that permuteOR works as the old code", { + + # Create matrices of 0 and 1 + A <- matrix(c(1, 1, 0, 1, 0, + 1, 1, 1, 0, 1, + 0, 1, 1, 0, 1, + 1, 0, 0, 1, 1, + 0, 1, 1, 1, 1), + nrow = 5, byrow = TRUE) + K <- matrix(c(1, 0, + 0, 1, + 1, 0, + 0, 1, + 1, 1), + nrow = 5, byrow = TRUE) + colnames(K) <- c("C1", "C2") + Gk <- K[, 1] %*% t(K[, 1]) + Gstar <- Gk[lower.tri(Gk)] + + # compute permuted odds ratios + set.seed(0) + res1 <- propr:::permuteOR(A, Gstar, 10) + set.seed(0) + res2 <- permuteOR_old(A, Gk, 10) + + # check + expect_true(all(res1 == res2)) +}) + test_that("check if same results are obtained compared to old code", { # Create matrices of 0 and 1 From 1a06c772daeb1880d413575504ab35de7c36b244 Mon Sep 17 00:00:00 2001 From: suzannejin Date: Fri, 13 Sep 2024 19:03:07 +0200 Subject: [PATCH 5/5] update news.md --- NEWS.md | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/NEWS.md b/NEWS.md index 95bf39a..b37c800 100755 --- a/NEWS.md +++ b/NEWS.md @@ -1,3 +1,8 @@ +## propr 5.1.1 +--------------------- +* Speed up `graflex` related functions + + ## propr 5.1.0 --------------------- * Allowed `updateCutoffs` function to compute the FDR values for negative cutoffs