diff --git a/R/calc_seq_p.R b/R/calc_seq_p.R index e123680..b8deed4 100644 --- a/R/calc_seq_p.R +++ b/R/calc_seq_p.R @@ -80,37 +80,38 @@ #' ) #' } calc_seq_p <- function( - test_analysis = 2, - test_hypothesis = "H1, H2, H3", - p_obs = tibble::tibble( - analysis = 1:2, - H1 = c(0.02, 0.0015), - H2 = c(0.01, 0.01), - H3 = c(0.01, 0.004) + test_analysis = 2, + test_hypothesis = "H1, H2, H3", + p_obs = tibble::tibble( + analysis = 1:2, + H1 = c(0.02, 0.0015), + H2 = c(0.01, 0.01), + H3 = c(0.01, 0.004) + ), + alpha_spending_type = 2, + n_analysis = 2, + initial_weight = c(0.3, 0.3, 0.4), + transition_mat = matrix(c( + 0.0000000, 0.4285714, 0.5714286, + 0.4285714, 0.0000000, 0.5714286, + 0.5000000, 0.5000000, 0.0000000 + ), nrow = 3, byrow = TRUE), + z_corr = matrix( + c( + 1.0000000, 0.7627701, 0.6666667, 0.7071068, 0.5393599, 0.4714045, + 0.7627701, 1.0000000, 0.6992059, 0.5393599, 0.7071068, 0.4944132, + 0.6666667, 0.6992059, 1.0000000, 0.4714045, 0.4944132, 0.7071068, + 0.7071068, 0.5393599, 0.4714045, 1.0000000, 0.7627701, 0.6666667, + 0.5393599, 0.7071068, 0.4944132, 0.7627701, 1.0000000, 0.6992059, + 0.4714045, 0.4944132, 0.7071068, 0.6666667, 0.6992059, 1.0000000 ), - alpha_spending_type = 2, - n_analysis = 2, - initial_weight = c(0.3, 0.3, 0.4), - transition_mat = matrix(c( - 0.0000000, 0.4285714, 0.5714286, - 0.4285714, 0.0000000, 0.5714286, - 0.5000000, 0.5000000, 0.0000000 - ), nrow = 3, byrow = TRUE), - z_corr = matrix( - c( - 1.0000000, 0.7627701, 0.6666667, 0.7071068, 0.5393599, 0.4714045, - 0.7627701, 1.0000000, 0.6992059, 0.5393599, 0.7071068, 0.4944132, - 0.6666667, 0.6992059, 1.0000000, 0.4714045, 0.4944132, 0.7071068, - 0.7071068, 0.5393599, 0.4714045, 1.0000000, 0.7627701, 0.6666667, - 0.5393599, 0.7071068, 0.4944132, 0.7627701, 1.0000000, 0.6992059, - 0.4714045, 0.4944132, 0.7071068, 0.6666667, 0.6992059, 1.0000000 - ), - nrow = 6, byrow = TRUE - ), - spending_fun = gsDesign::sfHSD, - spending_fun_par = -4, - info_frac = c(0.5, 1), - interval = c(1e-4, 0.2)) { + nrow = 6, byrow = TRUE + ), + spending_fun = gsDesign::sfHSD, + spending_fun_par = -4, + info_frac = c(0.5, 1), + interval = c(1e-4, 0.2) +) { foo <- function(x) { all_hypothesis <- strsplit(test_hypothesis, split = ", ") %>% unlist() all_hypothesis_idx <- as.numeric(gsub(".*?([0-9]+).*", "\\1", all_hypothesis)) diff --git a/R/generate_corr.R b/R/generate_corr.R index 947784a..0aa7036 100644 --- a/R/generate_corr.R +++ b/R/generate_corr.R @@ -68,13 +68,18 @@ generate_corr <- function(event) { D <- diag(elem$Event) # Within hypothesis across analyses + # D[Hi_As, Hi_At] = events for hypothesis i at min(s, t) for (i in 1:n_hypotheses) { - for (j in 2:n_analyses) { - count <- as.numeric(event %>% - filter(H1 == i & H2 == i & Analysis == j - 1) %>% - select(Event)) - D[i, n_hypotheses * (j - 1) + i] <- count - D[n_hypotheses * (j - 1) + i, i] <- count + for (s in 1:n_analyses) { + for (t in 1:n_analyses) { + if (s == t) next + count <- as.numeric(event %>% + filter(H1 == i & H2 == i & Analysis == min(s, t)) %>% + select(Event)) + row_idx <- n_hypotheses * (s - 1) + i + col_idx <- n_hypotheses * (t - 1) + i + D[row_idx, col_idx] <- count + } } } diff --git a/tests/testthat/test-independent-generate_bounds.R b/tests/testthat/test-independent-generate_bounds.R index 67d474b..97b269a 100644 --- a/tests/testthat/test-independent-generate_bounds.R +++ b/tests/testthat/test-independent-generate_bounds.R @@ -391,7 +391,6 @@ test_that("BH bounds replicate tables A6 and A7", { expect_equal(wA6_result1_test, round(wA6_result1, 4)) - wA6_result2 <- c( 0.00019939, 0.000422943, @@ -484,7 +483,6 @@ test_that("BH bounds replicate tables A6 and A7", { expect_equal(A7_result2_test, round(A7_result2, 4)) - A7_result3 <- c( 3.582064348, NA, @@ -507,9 +505,6 @@ test_that("BH bounds replicate tables A6 and A7", { expect_equal(A7_result3_test, round(A7_result3, 4)) - - - # Table A7 wA7_result1 <- c( 3.509232997, @@ -556,7 +551,6 @@ test_that("BH bounds replicate tables A6 and A7", { expect_equal(wA7_result2_test, round(wA7_result2, 4)) - wA7_result3 <- c( 3.570376445, NA, diff --git a/tests/testthat/test-independent-generate_corr.R b/tests/testthat/test-independent-generate_corr.R index a1699f5..d3a6a7f 100644 --- a/tests/testthat/test-independent-generate_corr.R +++ b/tests/testthat/test-independent-generate_corr.R @@ -29,3 +29,47 @@ test_that("2 endpoints 2 analysis correlation as expected", { expect_equal(matrix(corr %>% as.numeric(), nrow = 4, byrow = TRUE), corr_test) }) + +test_that("2 endpoints 3 analyses correlation is valid and correct", { + # Event counts: H1 (subgroup), H2 (overall), intersection = H1 (nested) + a <- c(50, 80, 100) # H1 events at analyses 1, 2, 3 + b <- c(80, 130, 170) # H2 events at analyses 1, 2, 3 + ab <- c(50, 80, 100) # Intersection events (= H1 for nested populations) + + event <- tibble::tribble( + ~H1, ~H2, ~Analysis, ~Event, + 1, 1, 1, a[1], 2, 2, 1, b[1], 1, 2, 1, ab[1], + 1, 1, 2, a[2], 2, 2, 2, b[2], 1, 2, 2, ab[2], + 1, 1, 3, a[3], 2, 2, 3, b[3], 1, 2, 3, ab[3] + ) + + corr <- generate_corr(event) + + # Matrix should be 6x6 + expect_equal(nrow(corr), 6) + expect_equal(ncol(corr), 6) + + # All diagonal entries should be 1 + expect_equal(diag(corr), rep(1, 6)) + + # All entries should be in [-1, 1] (with floating point tolerance) + expect_true(all(corr >= -1 - 1e-10 & corr <= 1 + 1e-10)) + + # Matrix should be positive definite + expect_true(all(eigen(corr)$values > 0)) + + # Check specific entries (use tolerance for named vector comparison): + # corr(H1_A1, H1_A2) = a[1] / sqrt(a[1] * a[2]) + expect_equal(as.numeric(corr[1, 3]), a[1] / sqrt(a[1] * a[2])) + # corr(H1_A1, H1_A3) = a[1] / sqrt(a[1] * a[3]) + expect_equal(as.numeric(corr[1, 5]), a[1] / sqrt(a[1] * a[3])) + # corr(H1_A2, H1_A3) = a[2] / sqrt(a[2] * a[3]) + expect_equal(as.numeric(corr[3, 5]), a[2] / sqrt(a[2] * a[3])) + # corr(H2_A1, H2_A3) = b[1] / sqrt(b[1] * b[3]) + expect_equal(as.numeric(corr[2, 6]), b[1] / sqrt(b[1] * b[3])) + + # Cross-hypothesis: corr(H1_A1, H2_A2) = ab[1] / sqrt(a[1] * b[2]) + expect_equal(as.numeric(corr[1, 4]), ab[1] / sqrt(a[1] * b[2])) + # corr(H1_A2, H2_A3) = ab[2] / sqrt(a[2] * b[3]) + expect_equal(as.numeric(corr[3, 6]), ab[2] / sqrt(a[2] * b[3])) +})