Fix generate_corr() for k > 2 analyses

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))

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
+      }
     }
   }
 

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,

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]))
+})