Fix index mapping bug in bpbounds_calc_tri_z3 and update examples

bpbounds_calc_tri_z3 was assigning p[2]→p100 and p[3]→p010 (and the corresponding pairs for z=1 and z=2), swapping the x=0,y=1 and x=1,y=0 cells relative to the (x,y,z) input ordering used by xtabs and bpbounds_calc_tri_z2. Corrected to match bpbounds_calc_tri_z2.

Replace the hardcoded approximate conditional probabilities for the Meleady trivariate example with count-data-derived probabilities via xtabs, which satisfy the IV constraints under the corrected mapping. Update test expectations accordingly.

Closes #3

Author: remlapmot

NEWS.md

@@ -4,6 +4,8 @@
 
 * **dplyr** removed as a soft dependency and instances of the magrittr/dplyr pipe, `%>%`, replaced with the native pipe, `|>`.
 
+* Fixed incorrect index mapping in `bpbounds_calc_tri_z3()`: the x=0,y=1 and x=1,y=0 conditional probability cells were swapped for each category of Z, giving wrong bounds for the trivariate 3-category instrument case (thanks @sachsmc).
+
 # bpbounds 0.1.6
 
 * Tweak formatting of code in helpfile examples and vignette

R/bpbounds_calc_tri_z3.R

@@ -1,17 +1,18 @@
 bpbounds_calc_tri_z3 <- function(p) {
   # assuming that p is in the following order
   # notation for conditional probabilities is p(y,x|z)
+  # input vector uses (x,y,z) ordering, consistent with bpbounds_calc_tri_z2
   p000 = p[1]
-  p100 = p[2]
-  p010 = p[3]
+  p100 = p[3]
+  p010 = p[2]
   p110 = p[4]
   p001 = p[5]
-  p101 = p[6]
-  p011 = p[7]
+  p101 = p[7]
+  p011 = p[6]
   p111 = p[8]
   p002 = p[9]
-  p102 = p[10]
-  p012 = p[11]
+  p102 = p[11]
+  p012 = p[10]
   p112 = p[12]
 
   # bounds on probabilities

tests/testthat/test-bpbounds.R

@@ -231,44 +231,37 @@ test_that("Balke and Pearl Table 2 example: trivariate data with 2 category inst
           })
 
 
-# Meleady AJCN 2003; Trivariate data with 3 category instrument - Table 3 of paper ----
+# Meleady AJCN 2003; 3 category instrument - Table 3 of paper ----
+dat <- data.frame(
+  count = c(341, 47, 297, 17, 63, 18, 272, 41, 269, 38, 56, 35),
+  z = c(0, 0, 1, 1, 2, 2, 0, 0, 1, 1, 2, 2),
+  y = c(0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1),
+  x = c(0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1)
+)
+longdat <- tidyr::uncount(dat, weights = count)
+
 ## Trivariate data
-mt3 <- c(.83, .05, .11, .01, .88, .06, .05, .01, .72, .05, .20, .03)
-p3 <- array(mt3,
-           dim = c(2, 2, 3),
-           dimnames = list(
-             x = c(0, 1),
-             y = c(0, 1),
-             z = c(0, 1, 2)
-           ))
-p3 <- as.table(p3)
+xt3 <- xtabs(count ~ x + y + z, data = dat)
+p3  <- prop.table(xt3, margin = 3)
 
 test_that("Mendelian randomization with 3 category instrument, trivariate data",
           {
             bpres <- bpbounds(p3)
             expect_true(bpres$inequality)
-            expect_equal(bpres$bplb, -0.090, tol = 1e-4)
-            expect_equal(bpres$bpub, 0.74, tol = 1e-4)
-            expect_equal(bpres$p10low, 0.06, tol = 1e-4)
-            expect_equal(bpres$p10upp, 0.12, tol = 1e-4)
-            expect_equal(bpres$p11low, 0.03, tol = 1e-4)
-            expect_equal(bpres$p11upp, 0.800, tol = 1e-4)
-            expect_equal(bpres$crrlb, 0.25, tol = 1e-4)
-            expect_equal(bpres$crrub, 13.3333, tol = 1e-4)
+            expect_equal(bpres$bplb, -0.3101, tol = 1e-4)
+            expect_equal(bpres$bpub, 0.4622, tol = 1e-4)
+            expect_equal(bpres$p10low, 0.4332, tol = 1e-4)
+            expect_equal(bpres$p10upp, 0.5136, tol = 1e-4)
+            expect_equal(bpres$p11low, 0.2035, tol = 1e-4)
+            expect_equal(bpres$p11upp, 0.8953, tol = 1e-4)
+            expect_equal(bpres$crrlb, 0.3962, tol = 1e-4)
+            expect_equal(bpres$crrub, 2.0670, tol = 1e-4)
             expect_false(bpres$monoinequality)
             sbp <- summary(bpres)
             print(sbp)
           })
 
 ## Bivariate data
-dat <- data.frame(
-  count = c(341, 47, 297, 17, 63, 18, 272, 41, 269, 38, 56, 35),
-  z = c(0, 0, 1, 1, 2, 2, 0, 0, 1, 1, 2, 2),
-  y = c(0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1),
-  x = c(0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1)
-)
-longdat <- tidyr::uncount(dat, weights = count)
-
 gtab <- xtabs(~ y + z, data = longdat)
 gp <- prop.table(gtab, margin = 2)
 gp

vignettes/bpbounds.Rmd

@@ -138,18 +138,17 @@ Mendelian randomization is an approach in epidemiology due to @daveysmith-ije-20
 This example uses data from @meleady-ajcn-2003. It is trivariate data with a 3-category instrument and binary phenotype and outcomes. The instrument is the 677CT polymorphism (rs1801133) in the Methylenetetrahydrofolate Reductase gene, involved in folate metabolism, as an instrumental variable to investigate the effect of homocysteine on cardiovascular
 disease. 
 
-The data are presented to us as conditional probabilities, so we take care to enter them in the correct position in the vectors.
+The data from Table 3 of @meleady-ajcn-2003 are given as counts below.
 
 ```{r}
-mt3 <- c(.83, .05, .11, .01, .88, .06, .05, .01, .72, .05, .20, .03)
-p3 <- array(mt3,
-            dim = c(2, 2, 3),
-            dimnames = list(
-              x = c(0, 1),
-              y = c(0, 1),
-              z = c(0, 1, 2)
-            ))
-p3     <- as.table(p3)
+dat <- data.frame(
+  count = c(341, 47, 297, 17, 63, 18, 272, 41, 269, 38, 56, 35),
+  z = c(0, 0, 1, 1, 2, 2, 0, 0, 1, 1, 2, 2),
+  y = c(0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1),
+  x = c(0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1)
+)
+xt3    <- xtabs(count ~ x + y + z, data = dat)
+p3     <- prop.table(xt3, margin = 3)
 bpres3 <- bpbounds(p3)
 sbp3   <- summary(bpres3)
 print(sbp3)