Nothing
tests/testthat/test-known-bounds.R
In causaloptim: An Interface to Specify Causal Graphs and Compute Bounds on Causal Effects
test_that("## simple unconfounded X -> Y", {
# b <- readRDS(test_path("test-graphs", "simple-unconfounded.RData"))
b <- igraph::graph_from_literal(X -+ Y, Ur -+ Y) |> initialize_graph()
eff <- "p{Y(X = 1) = 1} - p{Y(X = 0) = 1}"
obj <- analyze_graph(b, constraints = NULL, effectt = eff)
# with new version of optimize_effect:
new_version_bound <- optimize_effect_2(obj)
expect_equal(new_version_bound$bounds, c(lower = "\nMAX {\n p0_0 - p0_1\n}\n",
upper = "\nMIN {\n p0_0 - p0_1\n}\n"))
# visual comparison:
# cat(bound$bounds) # old version output string
# cat(new_version_bound$bounds) # new version output string
})
test_that("## simple confounded X -> Y", {
# b <- readRDS(test_path("test-graphs", "simple-confounded.RData"))
b <- igraph::graph_from_literal(X -+ Y, Ur -+ X, Ur -+ Y) |> initialize_graph()
eff <- "p{Y(X = 1) = 1} - p{Y(X = 0) = 1}"
obj <- analyze_graph(b, constraints = NULL, effectt = eff)
# with new version of optimize_effect:
new_version_bound <- optimize_effect_2(obj)
expect_equal(new_version_bound$bounds, c(lower = "\nMAX {\n -p10_ - p01_\n}\n",
upper = "\nMIN {\n 1 - p10_ - p01_\n}\n"))
# visual comparison:
# cat(bound$bounds) # old version output string
# cat(new_version_bound$bounds) # new version output string
## with categorical outcome (3 levels)
V(b)$nvals <- c(3,2,2)
eff <- "p{Y(X = 2) = 1} - p{Y(X = 0) = 1}"
obj <- analyze_graph(b, constraints = NULL, effectt = eff)
bound <- optimize_effect_2(obj)
expect_equal(bound$bounds, c(lower = "\nMAX {\n -p10_ - p20_ - p01_ - p11_\n}\n",
upper = "\nMIN {\n 1 - p20_ - p01_\n}\n"))
})
test_that("instrumental variable", {
## instrument Z -> X -> Y
#b <- readRDS(test_path("test-graphs", "instrument.RData"))
b <- igraph::graph_from_literal(Z -+ X, X -+ Y, Ul -+ Z, Ur -+ X, Ur -+ Y) |>
initialize_graph()
eff <- "p{Y(X = 1) = 1} - p{Y(X = 0) = 1}"
obj <- analyze_graph(b, constraints = NULL, effectt = eff)
# with new version of optimize_effect:
new_version_bound <- optimize_effect_2(obj)
expect_equal(new_version_bound$bounds, c(lower = "\nMAX {\n p00_0 - p00_1 - p10_1 - p01_1,\n p00_0 - p00_1 - p10_0 - p10_1 - p01_0,\n p00_0 - p00_1 + p10_0 - 2p10_1 - 2p01_1,\n -p10_1 - p01_1,\n -p10_0 - p01_0,\n -p00_0 + p00_1 - 2p10_0 + p10_1 - 2p01_0,\n -p00_0 + p00_1 - p10_0 - p10_1 - p01_1,\n -p00_0 + p00_1 - p10_0 - p01_0\n}\n",
upper = "\nMIN {\n 1 - p10_1 - p01_0,\n 1 + p00_0 + p10_0 - 2p10_1 - p01_1,\n 2 - p00_1 - p10_0 - p10_1 - 2p01_0,\n 1 - p10_1 - p01_1,\n 1 - p10_0 - p01_0,\n 1 + p00_1 - 2p10_0 + p10_1 - p01_0,\n 2 - p00_0 - p10_0 - p10_1 - 2p01_1,\n 1 - p10_0 - p01_1\n}\n"))
# visual comparison:
#cat(bound$bounds) # old version output string
#cat(new_version_bound$bounds) # new version output string
## with monotonocity
mono <- list("X(Z = 1) >= X(Z = 0)")
obj <- analyze_graph(b, constraints = mono, effectt = eff)
# with new version of optimize_effect:
new_version_bound <- optimize_effect_2(obj)
expect_equal(new_version_bound$bounds, c(lower = "\nMAX {\n p00_0 - p00_1 - p10_1 - p01_1\n}\n",
upper = "\nMIN {\n 1 - p10_1 - p01_0\n}\n"))
# visual comparison:
#cat(bound$bounds) # old version output string
#cat(new_version_bound$bounds) # new version output string
## compliers minus defiers effect
eff <- "p{X(Z = 1) = 1; X(Z = 0) = 0} - p{X(Z = 1) = 0; X(Z = 0) = 1}"
obj <- analyze_graph(b, constraints = NULL, effectt = eff)
# with new version of optimize_effect:
new_version_bound <- optimize_effect_2(obj)
expect_equal(new_version_bound$bounds, c(lower = "\nMAX {\n p00_0 - p00_1 + p01_0 - p01_1\n}\n",
upper = "\nMIN {\n p00_0 - p00_1 + p01_0 - p01_1\n}\n"))
# visual comparison:
#cat(bound$bounds) # old version output string
#cat(new_version_bound$bounds) # new version output string
## treatment effect among the treated? This one is not actually liner under the dag
eff <- "p{Y(X = 1) = 1; X = 1} - p{Y(X = 0) = 1; X = 1}"
ivmod <- create_causalmodel(graph = b, respvars = NULL, p.vals = expand.grid(Z = 0:1, X = 0:1, Y = 0:1),
prob.form = list(out = c("X", "Y"), cond = "Z"))
expect_false(check_linear_objective(ivmod,eff))
expect_error(obj <- analyze_graph(b, constraints = NULL, effectt = eff))
# bound <- optimize_effect(obj)
# with new version of optimize_effect:
# new_version_bound <- optimize_effect_2(obj)
#all(new_version_bound$bounds == c("\nMAX {\n p00_0 - p00_1 - p10_0 - p10_1,\n p00_0 - p00_1 - p10_1 - p01_1,\n p00_0 - p00_1 + p10_0 - 2p10_1 - 2p01_1,\n -p10_0 - p01_0,\n -p10_1 - p01_1,\n -p00_0 + p00_1 - 2p10_0 + p10_1 - 2p01_0,\n -p00_0 + p00_1 - p10_0 - p10_1,\n -p00_0 + p00_1 - p10_0 - p01_0\n}\n", "\nMIN {\n 1 - p10_1 - p01_0,\n 2 - p00_0 - p00_1 - p10_0 - p10_1 - 2p01_0,\n 1 + p00_0 + p10_0 - 2p10_1 - p01_1,\n 1 - p10_1 - p01_1,\n 1 - p10_0 - p01_0,\n 1 + p00_1 - 2p10_0 + p10_1 - p01_0,\n 1 - p10_0 - p01_1,\n 2 - p00_0 - p00_1 - p10_0 - p10_1 - 2p01_1\n}\n"))
# visual comparison:
#cat(bound$bounds) # old version output string
#cat(new_version_bound$bounds) # new version output string
})
test_that("Mediator", {
## mediator X -> Z -> Y
## bounds should match https://onlinelibrary.wiley.com/doi/full/10.1111/j.1541-0420.2007.00949.x
#b <- readRDS(test_path("test-graphs", "mediator.RData"))
b <- igraph::graph_from_literal(X -+ Z, Z -+ Y, X -+ Y, Ur -+ Z, Ur -+ Y) |>
initialize_graph()
## total effect: identifiable
eff <- "p{Y(Z(X = 1), X = 1) = 1}"
obj <- analyze_graph(b, constraints = NULL, effectt = eff)
#bound1 <- optimize_effect(obj) # with old version
# with new version of optimize_effect:
new_version_bound1 <- optimize_effect_2(obj)
expect_equal(new_version_bound1$bounds,
c(lower = "\nMAX {\n 1 - p00_1 - p10_1\n}\n", upper= "\nMIN {\n 1 - p00_1 - p10_1\n}\n"))
# visual comparison:
#cat(bound1$bounds) # old version output string
#cat(new_version_bound1$bounds) # new version output string
eff <- "p{Y(X = 1, Z(X = 1)) = 1}"
obj <- analyze_graph(b, constraints = NULL, effectt = eff)
#bound2 <- optimize_effect(obj) # with old version
# with new version of optimize_effect:
new_version_bound2 <- optimize_effect_2(obj)
# visual comparison:
#cat(bound2$bounds) # old version output string
#cat(new_version_bound2$bounds) # new version output string
#all(bound1$bounds == bound2$bounds)
expect_true(all(new_version_bound1$bounds == new_version_bound2$bounds)) # with new version
## controlled direct effect
eff <- "p{Y(X = 1, Z = 0) = 1} - p{Y(X = 0, Z = 0) = 1}"
obj <- analyze_graph(b, constraints = NULL, effectt = eff)
#bound <- optimize_effect(obj) # with old version
## check against equation (3)
#all(bound$bounds == c("\nMAX {\np00_0 + p01_1 - 1\n}\n\n", "\nMIN {\n- p00_1 - p01_0 + 1\n}\n\n"))
# with new version of optimize_effect:
new_version_bound <- optimize_effect_2(obj)
expect_true(all(new_version_bound$bounds == c("\nMAX {\n -1 + p00_0 + p01_1\n}\n", "\nMIN {\n 1 - p00_1 - p01_0\n}\n")))
# visual comparison:
#cat(bound$bounds) # old version output string
#cat("\nMAX {\np00_0 + p01_1 - 1\n}\n\n", "\nMIN {\n- p00_1 - p01_0 + 1\n}\n\n")#
#cat(new_version_bound$bounds) # new version output string
## with monotoncity
mono2 <- list("Z(X = 1) >= Z(X = 0)",
"Y(X = 1, Z = 0) >= Y(X = 0, Z = 0)",
"Y(X = 1, Z = 1) >= Y(X = 0, Z = 1)",
"Y(X = 0, Z = 1) >= Y(X = 0, Z = 0)",
"Y(X = 1, Z = 1) >= Y(X = 1, Z = 0)")
obj <- analyze_graph(b, constraints = mono2, effectt = eff)
#bound <- optimize_effect(obj) # with old version
## check againts equation (5)
#all(bound$bounds == c("\nMAX {\n- p01_0 + p01_1\n0\n}\n\n", "\nMIN {\n- p00_1 - p10_1 - p01_0 + 1\n}\n\n"))
# with new version of optimize_effect:
new_version_bound <- optimize_effect_2(obj)
expect_true(all(new_version_bound$bounds == c("\nMAX {\n -p01_0 + p01_1,\n 0\n}\n", "\nMIN {\n 1 - p00_1 - p10_1 - p01_0\n}\n")))
# visual comparison:
#cat(bound$bounds) # old version output string
#cat("\nMAX {\n- p01_0 + p01_1\n0\n}\n\n", "\nMIN {\n- p00_1 - p10_1 - p01_0 + 1\n}\n\n")#
#cat(new_version_bound$bounds) # new version output string
## natural direct effect
eff <- "p{Y(X = 1, Z(X = 0)) = 1} - p{Y(X = 0, Z(X = 0)) = 1}"
obj <- analyze_graph(b, constraints = NULL, effectt = eff)
#bound <- optimize_effect(obj) # with old version
# with new version of optimize_effect:
new_version_bound <- optimize_effect_2(obj)
#all(new_version_bound$bounds == c("min-bound: MAX {\n -p00_1 + p10_0 - p10_1 - p01_0 - p01_1,\n -2 + 2p00_0 + p10_0 + p01_0 + p01_1,\n -1 + p00_0 + p10_0\n}\n", "max-bound: MIN {\n 2p00_0 + p10_0 - p10_1 + p01_0,\n p00_0 + p10_0,\n 1 - p00_1 + p10_0 - p01_0\n}\n"))
# visual comparison:
#cat(bound$bounds) # old version output string
#cat(new_version_bound$bounds) # new version output string
obj <- analyze_graph(b, constraints = mono2, effectt = eff)
#bound <- optimize_effect(obj) # with old version
# with new version of optimize_effect:
new_version_bound <- optimize_effect_2(obj)
#all(new_version_bound$bounds == c("min-bound: MAX {\n p10_0 - p10_1 - p01_0 + p01_1,\n p10_0 - p10_1,\n -p01_0 + p01_1,\n 0\n}\n", "max-bound: MIN {\n p00_0 - p00_1 + p10_0 - p10_1\n}\n"))
# visual comparison:
#cat(bound$bounds) # old version output string
#cat(new_version_bound$bounds) # new version output string
## compare to Arvids paper 2009
obj <- analyze_graph(b, constraints = list("Z(X=0)<=Z(X=1)"), effectt = eff)
#bound <- optimize_effect(obj) # with old version
# with new version of optimize_effect:
new_version_bound <- optimize_effect_2(obj)
#all(new_version_bound$bounds == c("min-bound: MAX {\n p10_0 - p10_1 - p01_0 + p01_1,\n -1 + p00_0 + p10_0 + p01_1\n}\n", "max-bound: MIN {\n p00_0 - p00_1 + p10_0,\n 2p00_0 - 2p00_1 + p10_0 - p10_1 + p01_0 - p01_1\n}\n"))
# visual comparison:
#cat(bound$bounds) # old version output string
#cat(new_version_bound$bounds) # new version output string
})
test_that("Multiple IV numeric comparison", {
## comparison of multiple IV bounds results
b <- graph_from_literal(Z1 -+ X, Z2 -+ X, Z2 -+ Z1, Ul -+ Z1, Ul -+ Z2, X -+ Y, Ur -+ X, Ur -+ Y) |> initialize_graph()
# V(b)$leftside <- c(1, 0, 1, 1, 0, 0)
# V(b)$latent <- c(0, 0, 0, 1, 0, 1)
# V(b)$nvals <- c(2, 2, 2, 2, 2, 2)
# E(b)$rlconnect <- c(0, 0, 0, 0, 0, 0, 0, 0)
# E(b)$edge.monotone <- c(0, 0, 0, 0, 0, 0, 0, 0)
mivmod <- create_causalmodel(graph = b, p.vals = expand.grid(Z1 = 0:1, Z2 = 0:1,
X = 0:1, Y = 0:1),
prob.form = list(out = c("X", "Y"), cond = c("Z1", "Z2")))
obj <- analyze_graph(b, constraints = NULL, effectt = "p{Y(X = 1) = 1} - p{Y(X = 0) = 1}")
mivold <- readRDS(test_path("test-graphs", "MIV-bounds-result.RData"))
mivnew <- optimize_effect_2(obj)
f2new <- interpret_bounds(mivnew$bounds, obj$parameters)
f2old <- interpret_bounds(mivold$bounds, obj$parameters)
sim.qs <- rbeta(length(obj$variables), .05, 1)
sim.qs <- sim.qs / sum(sim.qs)
names(sim.qs) <- obj$variables
inenv <- new.env()
for(j in 1:length(sim.qs)) {
assign(names(sim.qs)[j], sim.qs[j], inenv)
}
res <- lapply(as.list(obj$constraints[-1]), function(x){
x1 <- strsplit(x, " = ")[[1]]
x0 <- paste(x1[1], " = ", x1[2])
eval(parse(text = x0), envir = inenv)
})
params <- lapply(obj$parameters, function(x) get(x, envir = inenv))
names(params) <- obj$parameters
expect_equal(do.call(f2new, params), do.call(f2old, params))
})
test_that("Missing paths filling in works properly", {
b2 <- graph_from_literal(X -+ Y, Ul -+ X, X -+ M1, X -+ M2, M1 -+ Y, M2 -+ Y,
Ur -+ M1, Ur -+ M2, Ur -+ Y, M1 -+ M2) |>
initialize_graph()
# V(b2)$leftside <- c(1, 0, 1, 0, 0, 0)
# V(b2)$latent <- c(0, 0, 1, 0, 0, 1)
# V(b2)$nvals <- c(2, 2, 2, 2, 2, 2)
# E(b2)$rlconnect <- rep(0, 10)
# E(b2)$edge.monotone <- rep(0, 10)
nofill <- "p{Y(X = 1, M1 = 1, M2(X = 1, M1(X = 1))) = 1}"
nofill <- "p{Y(X = 1, M1 = 1, M2(X = 1, M1 = 1)) = 1}"
eff2 <- parse_effect(nofill)$vars[[1]][[1]]
thisintervene <- unlist(list_to_path(eff2, "Y"))
basevars <- sapply(strsplit(names(thisintervene), " -> "), "[", 1)
## check for missing paths from intervention sets to outcome
outcome <- V(b2)[2]
parents <- adjacent_vertices(b2, v = outcome, mode = "in")[[1]]
expect_true(length(setdiff(names(parents[which(names(parents) != "Ur")]),
names(eff2))) == 0)
if(length(setdiff(names(parents[which(names(parents) != "Ur")]),
names(eff2))) > 0) {
isets <- unique(btm_var(eff2))
missingpaths <- lapply(isets, function(cc) {
allpaths <- igraph::all_simple_paths(b2, from = cc, to = "Y", mode = "out")
paths2 <- unlist(lapply(allpaths, function(x) paste(names(x), collapse = " -> ")))
setdiff(paths2, names(thisintervene))
})
}
b1 <- graph_from_literal(X -+ M1, M1 -+ Y, X -+ Y, Ul -+ X, Ur -+ M1, Ur -+ Y) |>
initialize_graph()
# V(b1)$leftside <- c(1, 0, 0, 1, 0)
# V(b1)$latent <- c(0, 0, 0, 1, 1)
# V(b1)$nvals <- c(2, 2, 2, 2, 2)
# E(b1)$rlconnect <- rep(0, 6)
# E(b1)$edge.monotone <- rep(0, 6)
fill <- "p{Y(X = 1) = 1}"
eff1 <- parse_effect(fill)$vars[[1]][[1]]
outcome <- V(b1)[3]
thisintervene <- unlist(list_to_path(eff1, "Y"))
basevars <- sapply(strsplit(names(thisintervene), " -> "), "[", 1)
## check for missing paths from intervention sets to outcome
## only do this if any of the top level intervention sets doesn't contain all
## parents of the outcome
## the logic being that if the user wrote that as an effect, then
## the intention was to propogate that intervention set forward through
## all paths in the graph to the outcome
parents <- adjacent_vertices(b1, v = outcome, mode = "in")[[1]]
if(length(setdiff(names(parents[which(names(parents) != "Ur")]),
names(eff1))) > 0) {
isets <- unique(btm_var(eff1))
missingpaths <- lapply(isets, function(cc) {
allpaths <- igraph::all_simple_paths(b1, from = cc, to = "Y", mode = "out")
paths2 <- unlist(lapply(allpaths, function(x) paste(names(x), collapse = " -> ")))
setdiff(paths2, names(thisintervene))
})
expect_equal(missingpaths[[1]], "X -> M1 -> Y")
}
})
test_that("IV with outcome dependent sampling", {
## see https://www.tandfonline.com/doi/full/10.1080/01621459.2020.1832502#d1e1777
## result 3
graph_casecont <- initialize_graph(graph_from_literal(X -+ Y, Y -+ S, Ur -+ X, Ur -+ Y, Ur -+ S))
prob.form <- list(out = c("X", "Y", "S"), cond = NULL)
p.vals <- expand.grid(X = 0:1, Y = 0:1, S = 1)
casecont <- create_causalmodel(graph_casecont, respvars= NULL, p.vals = p.vals, prob.form = prob.form)
bnds.cc <- create_linearcausalproblem(casecont, "p{Y(X = 1) = 1} - p{Y(X = 0) = 1}") |> optimize_effect_2()
expect_true(all(bnds.cc$bounds == c("\nMAX {\n -1 + p001_ + p111_\n}\n", "\nMIN {\n 1 - p101_ - p011_\n}\n")))
graph_ivcc <- initialize_graph(graph_from_literal(Z -+ X, X -+ Y, Y -+ S, Ur -+ X, Ur -+ Y, Ur -+ S))
prob.form <- list(out = c("X", "Y", "S"), cond = "Z")
p.vals <- expand.grid(X = 0:1, Y = 0:1, S = 1, Z = 0:1)
ivcc <- create_causalmodel(graph_ivcc, respvars= NULL, p.vals = p.vals, prob.form = prob.form)
bnds.ivcc <- create_linearcausalproblem(ivcc, "p{Y(X = 1) = 1} - p{Y(X = 0) = 1}") |> optimize_effect_2()
expect_true(all(bnds.ivcc$bounds == c("\nMAX {\n -1 - p011_0 - p111_0 + p001_1 + 2p111_1,\n -1 - p001_0 - p101_0 - p011_0 - p111_0 + 2p001_1 + 2p111_1,\n -1 + p001_0 + 2p111_0 - p011_1 - p111_1,\n -1 + 2p001_0 + 2p111_0 - p001_1 - p101_1 - p011_1 - p111_1,\n -1 + 2p001_0 + p111_0 - p001_1 - p101_1,\n -1 + p001_0 + p111_0,\n -1 + p001_1 + p111_1,\n -1 + p001_0 + p111_1,\n -1 + p111_0 + p001_1,\n -1 - p001_0 - p101_0 + 2p001_1 + p111_1\n}\n",
"\nMIN {\n 1 - 2p101_0 - 2p011_0 + p001_1 + p101_1 + p011_1 + p111_1,\n 1 - p101_1 - p011_1,\n 1 - p101_0 - p011_0,\n 1 + p001_0 + p101_0 + p011_0 + p111_0 - 2p101_1 - 2p011_1,\n 1 + p001_0 + p101_0 - 2p101_1 - p011_1,\n 1 - p011_0 - p101_1,\n 1 - p101_0 - 2p011_0 + p011_1 + p111_1,\n 1 - p101_0 - p011_1,\n 1 + p011_0 + p111_0 - p101_1 - 2p011_1,\n 1 - 2p101_0 - p011_0 + p001_1 + p101_1\n}\n" )))
})
Try the causaloptim package in your browser
Any scripts or data that you put into this service are public.
causaloptim documentation built on Oct. 17, 2024, 9:08 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
test_that("## simple unconfounded X -> Y", {
# b <- readRDS(test_path("test-graphs", "simple-unconfounded.RData"))
b <- igraph::graph_from_literal(X -+ Y, Ur -+ Y) |> initialize_graph()
eff <- "p{Y(X = 1) = 1} - p{Y(X = 0) = 1}"
obj <- analyze_graph(b, constraints = NULL, effectt = eff)
# with new version of optimize_effect:
new_version_bound <- optimize_effect_2(obj)
expect_equal(new_version_bound$bounds, c(lower = "\nMAX {\n p0_0 - p0_1\n}\n",
upper = "\nMIN {\n p0_0 - p0_1\n}\n"))
# visual comparison:
# cat(bound$bounds) # old version output string
# cat(new_version_bound$bounds) # new version output string
})
test_that("## simple confounded X -> Y", {
# b <- readRDS(test_path("test-graphs", "simple-confounded.RData"))
b <- igraph::graph_from_literal(X -+ Y, Ur -+ X, Ur -+ Y) |> initialize_graph()
eff <- "p{Y(X = 1) = 1} - p{Y(X = 0) = 1}"
obj <- analyze_graph(b, constraints = NULL, effectt = eff)
# with new version of optimize_effect:
new_version_bound <- optimize_effect_2(obj)
expect_equal(new_version_bound$bounds, c(lower = "\nMAX {\n -p10_ - p01_\n}\n",
upper = "\nMIN {\n 1 - p10_ - p01_\n}\n"))
# visual comparison:
# cat(bound$bounds) # old version output string
# cat(new_version_bound$bounds) # new version output string
## with categorical outcome (3 levels)
V(b)$nvals <- c(3,2,2)
eff <- "p{Y(X = 2) = 1} - p{Y(X = 0) = 1}"
obj <- analyze_graph(b, constraints = NULL, effectt = eff)
bound <- optimize_effect_2(obj)
expect_equal(bound$bounds, c(lower = "\nMAX {\n -p10_ - p20_ - p01_ - p11_\n}\n",
upper = "\nMIN {\n 1 - p20_ - p01_\n}\n"))
})
test_that("instrumental variable", {
## instrument Z -> X -> Y
#b <- readRDS(test_path("test-graphs", "instrument.RData"))
b <- igraph::graph_from_literal(Z -+ X, X -+ Y, Ul -+ Z, Ur -+ X, Ur -+ Y) |>
initialize_graph()
eff <- "p{Y(X = 1) = 1} - p{Y(X = 0) = 1}"
obj <- analyze_graph(b, constraints = NULL, effectt = eff)
# with new version of optimize_effect:
new_version_bound <- optimize_effect_2(obj)
expect_equal(new_version_bound$bounds, c(lower = "\nMAX {\n p00_0 - p00_1 - p10_1 - p01_1,\n p00_0 - p00_1 - p10_0 - p10_1 - p01_0,\n p00_0 - p00_1 + p10_0 - 2p10_1 - 2p01_1,\n -p10_1 - p01_1,\n -p10_0 - p01_0,\n -p00_0 + p00_1 - 2p10_0 + p10_1 - 2p01_0,\n -p00_0 + p00_1 - p10_0 - p10_1 - p01_1,\n -p00_0 + p00_1 - p10_0 - p01_0\n}\n",
upper = "\nMIN {\n 1 - p10_1 - p01_0,\n 1 + p00_0 + p10_0 - 2p10_1 - p01_1,\n 2 - p00_1 - p10_0 - p10_1 - 2p01_0,\n 1 - p10_1 - p01_1,\n 1 - p10_0 - p01_0,\n 1 + p00_1 - 2p10_0 + p10_1 - p01_0,\n 2 - p00_0 - p10_0 - p10_1 - 2p01_1,\n 1 - p10_0 - p01_1\n}\n"))
# visual comparison:
#cat(bound$bounds) # old version output string
#cat(new_version_bound$bounds) # new version output string
## with monotonocity
mono <- list("X(Z = 1) >= X(Z = 0)")
obj <- analyze_graph(b, constraints = mono, effectt = eff)
# with new version of optimize_effect:
new_version_bound <- optimize_effect_2(obj)
expect_equal(new_version_bound$bounds, c(lower = "\nMAX {\n p00_0 - p00_1 - p10_1 - p01_1\n}\n",
upper = "\nMIN {\n 1 - p10_1 - p01_0\n}\n"))
# visual comparison:
#cat(bound$bounds) # old version output string
#cat(new_version_bound$bounds) # new version output string
## compliers minus defiers effect
eff <- "p{X(Z = 1) = 1; X(Z = 0) = 0} - p{X(Z = 1) = 0; X(Z = 0) = 1}"
obj <- analyze_graph(b, constraints = NULL, effectt = eff)
# with new version of optimize_effect:
new_version_bound <- optimize_effect_2(obj)
expect_equal(new_version_bound$bounds, c(lower = "\nMAX {\n p00_0 - p00_1 + p01_0 - p01_1\n}\n",
upper = "\nMIN {\n p00_0 - p00_1 + p01_0 - p01_1\n}\n"))
# visual comparison:
#cat(bound$bounds) # old version output string
#cat(new_version_bound$bounds) # new version output string
## treatment effect among the treated? This one is not actually liner under the dag
eff <- "p{Y(X = 1) = 1; X = 1} - p{Y(X = 0) = 1; X = 1}"
ivmod <- create_causalmodel(graph = b, respvars = NULL, p.vals = expand.grid(Z = 0:1, X = 0:1, Y = 0:1),
prob.form = list(out = c("X", "Y"), cond = "Z"))
expect_false(check_linear_objective(ivmod,eff))
expect_error(obj <- analyze_graph(b, constraints = NULL, effectt = eff))
# bound <- optimize_effect(obj)
# with new version of optimize_effect:
# new_version_bound <- optimize_effect_2(obj)
#all(new_version_bound$bounds == c("\nMAX {\n p00_0 - p00_1 - p10_0 - p10_1,\n p00_0 - p00_1 - p10_1 - p01_1,\n p00_0 - p00_1 + p10_0 - 2p10_1 - 2p01_1,\n -p10_0 - p01_0,\n -p10_1 - p01_1,\n -p00_0 + p00_1 - 2p10_0 + p10_1 - 2p01_0,\n -p00_0 + p00_1 - p10_0 - p10_1,\n -p00_0 + p00_1 - p10_0 - p01_0\n}\n", "\nMIN {\n 1 - p10_1 - p01_0,\n 2 - p00_0 - p00_1 - p10_0 - p10_1 - 2p01_0,\n 1 + p00_0 + p10_0 - 2p10_1 - p01_1,\n 1 - p10_1 - p01_1,\n 1 - p10_0 - p01_0,\n 1 + p00_1 - 2p10_0 + p10_1 - p01_0,\n 1 - p10_0 - p01_1,\n 2 - p00_0 - p00_1 - p10_0 - p10_1 - 2p01_1\n}\n"))
# visual comparison:
#cat(bound$bounds) # old version output string
#cat(new_version_bound$bounds) # new version output string
})
test_that("Mediator", {
## mediator X -> Z -> Y
## bounds should match https://onlinelibrary.wiley.com/doi/full/10.1111/j.1541-0420.2007.00949.x
#b <- readRDS(test_path("test-graphs", "mediator.RData"))
b <- igraph::graph_from_literal(X -+ Z, Z -+ Y, X -+ Y, Ur -+ Z, Ur -+ Y) |>
initialize_graph()
## total effect: identifiable
eff <- "p{Y(Z(X = 1), X = 1) = 1}"
obj <- analyze_graph(b, constraints = NULL, effectt = eff)
#bound1 <- optimize_effect(obj) # with old version
# with new version of optimize_effect:
new_version_bound1 <- optimize_effect_2(obj)
expect_equal(new_version_bound1$bounds,
c(lower = "\nMAX {\n 1 - p00_1 - p10_1\n}\n", upper= "\nMIN {\n 1 - p00_1 - p10_1\n}\n"))
# visual comparison:
#cat(bound1$bounds) # old version output string
#cat(new_version_bound1$bounds) # new version output string
eff <- "p{Y(X = 1, Z(X = 1)) = 1}"
obj <- analyze_graph(b, constraints = NULL, effectt = eff)
#bound2 <- optimize_effect(obj) # with old version
# with new version of optimize_effect:
new_version_bound2 <- optimize_effect_2(obj)
# visual comparison:
#cat(bound2$bounds) # old version output string
#cat(new_version_bound2$bounds) # new version output string
#all(bound1$bounds == bound2$bounds)
expect_true(all(new_version_bound1$bounds == new_version_bound2$bounds)) # with new version
## controlled direct effect
eff <- "p{Y(X = 1, Z = 0) = 1} - p{Y(X = 0, Z = 0) = 1}"
obj <- analyze_graph(b, constraints = NULL, effectt = eff)
#bound <- optimize_effect(obj) # with old version
## check against equation (3)
#all(bound$bounds == c("\nMAX {\np00_0 + p01_1 - 1\n}\n\n", "\nMIN {\n- p00_1 - p01_0 + 1\n}\n\n"))
# with new version of optimize_effect:
new_version_bound <- optimize_effect_2(obj)
expect_true(all(new_version_bound$bounds == c("\nMAX {\n -1 + p00_0 + p01_1\n}\n", "\nMIN {\n 1 - p00_1 - p01_0\n}\n")))
# visual comparison:
#cat(bound$bounds) # old version output string
#cat("\nMAX {\np00_0 + p01_1 - 1\n}\n\n", "\nMIN {\n- p00_1 - p01_0 + 1\n}\n\n")#
#cat(new_version_bound$bounds) # new version output string
## with monotoncity
mono2 <- list("Z(X = 1) >= Z(X = 0)",
"Y(X = 1, Z = 0) >= Y(X = 0, Z = 0)",
"Y(X = 1, Z = 1) >= Y(X = 0, Z = 1)",
"Y(X = 0, Z = 1) >= Y(X = 0, Z = 0)",
"Y(X = 1, Z = 1) >= Y(X = 1, Z = 0)")
obj <- analyze_graph(b, constraints = mono2, effectt = eff)
#bound <- optimize_effect(obj) # with old version
## check againts equation (5)
#all(bound$bounds == c("\nMAX {\n- p01_0 + p01_1\n0\n}\n\n", "\nMIN {\n- p00_1 - p10_1 - p01_0 + 1\n}\n\n"))
# with new version of optimize_effect:
new_version_bound <- optimize_effect_2(obj)
expect_true(all(new_version_bound$bounds == c("\nMAX {\n -p01_0 + p01_1,\n 0\n}\n", "\nMIN {\n 1 - p00_1 - p10_1 - p01_0\n}\n")))
# visual comparison:
#cat(bound$bounds) # old version output string
#cat("\nMAX {\n- p01_0 + p01_1\n0\n}\n\n", "\nMIN {\n- p00_1 - p10_1 - p01_0 + 1\n}\n\n")#
#cat(new_version_bound$bounds) # new version output string
## natural direct effect
eff <- "p{Y(X = 1, Z(X = 0)) = 1} - p{Y(X = 0, Z(X = 0)) = 1}"
obj <- analyze_graph(b, constraints = NULL, effectt = eff)
#bound <- optimize_effect(obj) # with old version
# with new version of optimize_effect:
new_version_bound <- optimize_effect_2(obj)
#all(new_version_bound$bounds == c("min-bound: MAX {\n -p00_1 + p10_0 - p10_1 - p01_0 - p01_1,\n -2 + 2p00_0 + p10_0 + p01_0 + p01_1,\n -1 + p00_0 + p10_0\n}\n", "max-bound: MIN {\n 2p00_0 + p10_0 - p10_1 + p01_0,\n p00_0 + p10_0,\n 1 - p00_1 + p10_0 - p01_0\n}\n"))
# visual comparison:
#cat(bound$bounds) # old version output string
#cat(new_version_bound$bounds) # new version output string
obj <- analyze_graph(b, constraints = mono2, effectt = eff)
#bound <- optimize_effect(obj) # with old version
# with new version of optimize_effect:
new_version_bound <- optimize_effect_2(obj)
#all(new_version_bound$bounds == c("min-bound: MAX {\n p10_0 - p10_1 - p01_0 + p01_1,\n p10_0 - p10_1,\n -p01_0 + p01_1,\n 0\n}\n", "max-bound: MIN {\n p00_0 - p00_1 + p10_0 - p10_1\n}\n"))
# visual comparison:
#cat(bound$bounds) # old version output string
#cat(new_version_bound$bounds) # new version output string
## compare to Arvids paper 2009
obj <- analyze_graph(b, constraints = list("Z(X=0)<=Z(X=1)"), effectt = eff)
#bound <- optimize_effect(obj) # with old version
# with new version of optimize_effect:
new_version_bound <- optimize_effect_2(obj)
#all(new_version_bound$bounds == c("min-bound: MAX {\n p10_0 - p10_1 - p01_0 + p01_1,\n -1 + p00_0 + p10_0 + p01_1\n}\n", "max-bound: MIN {\n p00_0 - p00_1 + p10_0,\n 2p00_0 - 2p00_1 + p10_0 - p10_1 + p01_0 - p01_1\n}\n"))
# visual comparison:
#cat(bound$bounds) # old version output string
#cat(new_version_bound$bounds) # new version output string
})
test_that("Multiple IV numeric comparison", {
## comparison of multiple IV bounds results
b <- graph_from_literal(Z1 -+ X, Z2 -+ X, Z2 -+ Z1, Ul -+ Z1, Ul -+ Z2, X -+ Y, Ur -+ X, Ur -+ Y) |> initialize_graph()
# V(b)$leftside <- c(1, 0, 1, 1, 0, 0)
# V(b)$latent <- c(0, 0, 0, 1, 0, 1)
# V(b)$nvals <- c(2, 2, 2, 2, 2, 2)
# E(b)$rlconnect <- c(0, 0, 0, 0, 0, 0, 0, 0)
# E(b)$edge.monotone <- c(0, 0, 0, 0, 0, 0, 0, 0)
mivmod <- create_causalmodel(graph = b, p.vals = expand.grid(Z1 = 0:1, Z2 = 0:1,
X = 0:1, Y = 0:1),
prob.form = list(out = c("X", "Y"), cond = c("Z1", "Z2")))
obj <- analyze_graph(b, constraints = NULL, effectt = "p{Y(X = 1) = 1} - p{Y(X = 0) = 1}")
mivold <- readRDS(test_path("test-graphs", "MIV-bounds-result.RData"))
mivnew <- optimize_effect_2(obj)
f2new <- interpret_bounds(mivnew$bounds, obj$parameters)
f2old <- interpret_bounds(mivold$bounds, obj$parameters)
sim.qs <- rbeta(length(obj$variables), .05, 1)
sim.qs <- sim.qs / sum(sim.qs)
names(sim.qs) <- obj$variables
inenv <- new.env()
for(j in 1:length(sim.qs)) {
assign(names(sim.qs)[j], sim.qs[j], inenv)
}
res <- lapply(as.list(obj$constraints[-1]), function(x){
x1 <- strsplit(x, " = ")[[1]]
x0 <- paste(x1[1], " = ", x1[2])
eval(parse(text = x0), envir = inenv)
})
params <- lapply(obj$parameters, function(x) get(x, envir = inenv))
names(params) <- obj$parameters
expect_equal(do.call(f2new, params), do.call(f2old, params))
})
test_that("Missing paths filling in works properly", {
b2 <- graph_from_literal(X -+ Y, Ul -+ X, X -+ M1, X -+ M2, M1 -+ Y, M2 -+ Y,
Ur -+ M1, Ur -+ M2, Ur -+ Y, M1 -+ M2) |>
initialize_graph()
# V(b2)$leftside <- c(1, 0, 1, 0, 0, 0)
# V(b2)$latent <- c(0, 0, 1, 0, 0, 1)
# V(b2)$nvals <- c(2, 2, 2, 2, 2, 2)
# E(b2)$rlconnect <- rep(0, 10)
# E(b2)$edge.monotone <- rep(0, 10)
nofill <- "p{Y(X = 1, M1 = 1, M2(X = 1, M1(X = 1))) = 1}"
nofill <- "p{Y(X = 1, M1 = 1, M2(X = 1, M1 = 1)) = 1}"
eff2 <- parse_effect(nofill)$vars[[1]][[1]]
thisintervene <- unlist(list_to_path(eff2, "Y"))
basevars <- sapply(strsplit(names(thisintervene), " -> "), "[", 1)
## check for missing paths from intervention sets to outcome
outcome <- V(b2)[2]
parents <- adjacent_vertices(b2, v = outcome, mode = "in")[[1]]
expect_true(length(setdiff(names(parents[which(names(parents) != "Ur")]),
names(eff2))) == 0)
if(length(setdiff(names(parents[which(names(parents) != "Ur")]),
names(eff2))) > 0) {
isets <- unique(btm_var(eff2))
missingpaths <- lapply(isets, function(cc) {
allpaths <- igraph::all_simple_paths(b2, from = cc, to = "Y", mode = "out")
paths2 <- unlist(lapply(allpaths, function(x) paste(names(x), collapse = " -> ")))
setdiff(paths2, names(thisintervene))
})
}
b1 <- graph_from_literal(X -+ M1, M1 -+ Y, X -+ Y, Ul -+ X, Ur -+ M1, Ur -+ Y) |>
initialize_graph()
# V(b1)$leftside <- c(1, 0, 0, 1, 0)
# V(b1)$latent <- c(0, 0, 0, 1, 1)
# V(b1)$nvals <- c(2, 2, 2, 2, 2)
# E(b1)$rlconnect <- rep(0, 6)
# E(b1)$edge.monotone <- rep(0, 6)
fill <- "p{Y(X = 1) = 1}"
eff1 <- parse_effect(fill)$vars[[1]][[1]]
outcome <- V(b1)[3]
thisintervene <- unlist(list_to_path(eff1, "Y"))
basevars <- sapply(strsplit(names(thisintervene), " -> "), "[", 1)
## check for missing paths from intervention sets to outcome
## only do this if any of the top level intervention sets doesn't contain all
## parents of the outcome
## the logic being that if the user wrote that as an effect, then
## the intention was to propogate that intervention set forward through
## all paths in the graph to the outcome
parents <- adjacent_vertices(b1, v = outcome, mode = "in")[[1]]
if(length(setdiff(names(parents[which(names(parents) != "Ur")]),
names(eff1))) > 0) {
isets <- unique(btm_var(eff1))
missingpaths <- lapply(isets, function(cc) {
allpaths <- igraph::all_simple_paths(b1, from = cc, to = "Y", mode = "out")
paths2 <- unlist(lapply(allpaths, function(x) paste(names(x), collapse = " -> ")))
setdiff(paths2, names(thisintervene))
})
expect_equal(missingpaths[[1]], "X -> M1 -> Y")
}
})
test_that("IV with outcome dependent sampling", {
## see https://www.tandfonline.com/doi/full/10.1080/01621459.2020.1832502#d1e1777
## result 3
graph_casecont <- initialize_graph(graph_from_literal(X -+ Y, Y -+ S, Ur -+ X, Ur -+ Y, Ur -+ S))
prob.form <- list(out = c("X", "Y", "S"), cond = NULL)
p.vals <- expand.grid(X = 0:1, Y = 0:1, S = 1)
casecont <- create_causalmodel(graph_casecont, respvars= NULL, p.vals = p.vals, prob.form = prob.form)
bnds.cc <- create_linearcausalproblem(casecont, "p{Y(X = 1) = 1} - p{Y(X = 0) = 1}") |> optimize_effect_2()
expect_true(all(bnds.cc$bounds == c("\nMAX {\n -1 + p001_ + p111_\n}\n", "\nMIN {\n 1 - p101_ - p011_\n}\n")))
graph_ivcc <- initialize_graph(graph_from_literal(Z -+ X, X -+ Y, Y -+ S, Ur -+ X, Ur -+ Y, Ur -+ S))
prob.form <- list(out = c("X", "Y", "S"), cond = "Z")
p.vals <- expand.grid(X = 0:1, Y = 0:1, S = 1, Z = 0:1)
ivcc <- create_causalmodel(graph_ivcc, respvars= NULL, p.vals = p.vals, prob.form = prob.form)
bnds.ivcc <- create_linearcausalproblem(ivcc, "p{Y(X = 1) = 1} - p{Y(X = 0) = 1}") |> optimize_effect_2()
expect_true(all(bnds.ivcc$bounds == c("\nMAX {\n -1 - p011_0 - p111_0 + p001_1 + 2p111_1,\n -1 - p001_0 - p101_0 - p011_0 - p111_0 + 2p001_1 + 2p111_1,\n -1 + p001_0 + 2p111_0 - p011_1 - p111_1,\n -1 + 2p001_0 + 2p111_0 - p001_1 - p101_1 - p011_1 - p111_1,\n -1 + 2p001_0 + p111_0 - p001_1 - p101_1,\n -1 + p001_0 + p111_0,\n -1 + p001_1 + p111_1,\n -1 + p001_0 + p111_1,\n -1 + p111_0 + p001_1,\n -1 - p001_0 - p101_0 + 2p001_1 + p111_1\n}\n",
"\nMIN {\n 1 - 2p101_0 - 2p011_0 + p001_1 + p101_1 + p011_1 + p111_1,\n 1 - p101_1 - p011_1,\n 1 - p101_0 - p011_0,\n 1 + p001_0 + p101_0 + p011_0 + p111_0 - 2p101_1 - 2p011_1,\n 1 + p001_0 + p101_0 - 2p101_1 - p011_1,\n 1 - p011_0 - p101_1,\n 1 - p101_0 - 2p011_0 + p011_1 + p111_1,\n 1 - p101_0 - p011_1,\n 1 + p011_0 + p111_0 - p101_1 - 2p011_1,\n 1 - 2p101_0 - p011_0 + p001_1 + p101_1\n}\n" )))
})
Try the causaloptim package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.