## ------------------------------------------------------------------------------------------------------------------------------------------------ #| include: true #' --- #' title: "Replication code for 'Making, Updating, and Querying Causal Models using CausalQueries'" #' date: `r Sys.Date()` #' author: Till Tietz, Lily Medina, Georgiy Syunyaev, and Macartan Humphreys #' format: #' html: #' toc: true #' toc-depth: 2 #' --- #' #' Generated by `knitr::spin()` from `code.R`. #' #' Note: final section of code for microbenchmarking is slow #' ## Set up ## ------------------------------------------------------------------------------------------------------------------------------------------------ #| label: preamble #| include: true library("tidyverse") library("CausalQueries") library("microbenchmark") library("parallel") library("future") library("future.apply") library("knitr") library("rstan") options(kableExtra.latex.load_packages = FALSE) library("kableExtra") options(mc.cores = parallel::detectCores()) set.seed(1, "L'Ecuyer-CMRG") theme_set(theme_bw()) kabble <- function( x, align = "c", format = ifelse(knitr::is_html_output(), "html", "latex"), ...) { x |> knitr::kable( format = format, digits = 2, align = align, booktabs = TRUE, linesep = "", longtable = TRUE, ...) |> kableExtra::kable_classic_2(latex_options = c("scale_down")) } ## ------------------------------------------------------------------------------------------------------------------------------------------------ #| include: true #' # Motivating example ## ------------------------------------------------------------------------------------------------------------------------------------------------ #| echo: true data("lipids_data") lipids_data ## ------------------------------------------------------------------------------------------------------------------------------------------------ #| include: true #' # Motivating example ## ------------------------------------------------------------------------------------------------------------------------------------------------ #| echo: true #| eval: true #| purl: true #| message: false lipids_model <- make_model("Z -> X -> Y; X <-> Y") ## ------------------------------------------------------------------------------------------------------------------------------------------------ #| echo: true #| eval: true #| purl: true #| message: false lipids_model <- update_model(lipids_model, lipids_data) ## ------------------------------------------------------------------------------------------------------------------------------------------------ #| echo: true #| eval: true #| purl: true lipids_queries <- query_model(lipids_model, queries = list( ATE = "Y[X = 1] - Y[X = 0]", PoC = "Y[X = 1] - Y[X = 0] :|: X == 0 & Y == 0", LATE = "Y[X = 1] - Y[X = 0] :|: X[Z = 1] > X[Z = 0]"), using = "posteriors") ## ------------------------------------------------------------------------------------------------------------------------------------------------ #| label: tbl-lipids #| tbl-cap: "Replication of \\citet{chickering_clinicians_1996}." #| echo: true lipids_queries |> dplyr::select(label, query, given, mean, sd, starts_with("cred")) |> kabble(align = "c") ## ------------------------------------------------------------------------------------------------------------------------------------------------ #| fig-cap: "Illustration of queries plotted." #| label: queryplot #| fig-width: 8 #| fig-height: 3 lipids_queries |> plot() ## ------------------------------------------------------------------------------------------------------------------------------------------------ #| echo: true #| include: false #| message: false # make_model("Z -> X -> Y; X <-> Y") |> update_model(lipids_data) |> # query_model(queries = list(ATE = "Y[X = 1] - Y[X = 0]", # PoC = "Y[X = 1] - Y[X = 0] :|: X == 0 & Y == 0", # LATE = "Y[X = 1] - Y[X = 0] :|: X[Z = 1] > X[Z = 0]"), # using = "posteriors") |> plot() # ## ------------------------------------------------------------------------------------------------------------------------------------------------ #| include: true #' # Statistical model ## ------------------------------------------------------------------------------------------------------------------------------------------------ #| label: tbl-lipidspar #| tbl-cap: "Nodal types and parameters for Lipids model." #| echo: true with_pars <- lipids_model |> set_parameters(param_type = "prior_draw") with_pars |> grab("parameters_df") |> dplyr::select(node, nodal_type, param_set, param_names, param_value, priors) |> kabble() ## ------------------------------------------------------------------------------------------------------------------------------------------------ #| include: true #' # Making models ## ------------------------------------------------------------------------------------------------------------------------------------------------ #| echo: true #| results: markup model <- make_model("X -> M -> Y <- X") ## ------------------------------------------------------------------------------------------------------------------------------------------------ #| include: true #' ## Graphing ## ------------------------------------------------------------------------------------------------------------------------------------------------ #| label: fig-plots #| echo: true #| out-width: 50% #| fig-height: 6 #| fig-width: 4 #| fig-cap: "Examples of model graphs." #| fig-subcap: #| - "Without options" #| - "With options" #| layout-ncol: 2 #| purl: true lipids_model |> plot() lipids_model |> plot(x_coord = 1:3, y_coord = 3:1, textcol = "black", textsize = 3, shape = c(15, 16, 16), nodecol = "lightgrey", nodesize = 10) ## ------------------------------------------------------------------------------------------------------------------------------------------------ #| include: true #' ## Tailoring models ## ------------------------------------------------------------------------------------------------------------------------------------------------ #| label: tbl-dof #| tbl-cap: "Number of different independent parameters (degrees of freedom) for different three-node models." #| tbl-pos: "b" #| echo: true #| eval: true statements <- list("X -> Y <- W", "X -> Y <- W; X <-> W", "X -> Y <- W; X <-> Y; W <-> Y", "X -> Y <- W; X <-> Y; W <-> Y; X <-> W", "X -> W -> Y <- X", "X -> W -> Y <- X; W <-> Y", "X -> W -> Y <- X; X <-> W; W <-> Y", "X -> W -> Y <- X; X <-> W; W <-> Y; X <-> Y") dof <- function(statement) { make_model(statement, add_causal_types = FALSE) |> grab("parameters_df") |> group_by(param_set) |> summarize(n = n() - 1) |> pull(n) |> sum() } data.frame( Model = c( "X -> Y <- W", "X -> Y <- W; X <-> W", "X -> Y <- W; X <-> Y; W <-> Y", "X -> Y <- W; X <-> Y; W <-> Y; X <-> W", "X -> W -> Y <- X", "X -> W -> Y <- X; W <-> Y", "X -> W -> Y <- X; X <-> W; W <-> Y", "X -> W -> Y <- X; X <-> W; W <-> Y; X <-> Y" ), dof = unlist(lapply(statements, dof)) ) |> kabble( col.names = c("Model", "Degrees of freedom"), align = c("l", "c") ) ## ------------------------------------------------------------------------------------------------------------------------------------------------ #| echo: true #| eval: true model_restricted <- lipids_model |> set_restrictions("X[Z = 1] < X[Z = 0]") ## ------------------------------------------------------------------------------------------------------------------------------------------------ #| label: set-restrictions1 #| echo: true #| purl: true model <- lipids_model |> set_restrictions( labels = list(X = "01", Y = c("00", "01", "11")), keep = TRUE) ## ------------------------------------------------------------------------------------------------------------------------------------------------ #| label: set-restrictions2 #| echo: true #| purl: true model <- lipids_model |> set_restrictions(labels = list(Y = "?0")) ## ------------------------------------------------------------------------------------------------------------------------------------------------ #| label: set-restrictions3 #| echo: true #| purl: true model <- lipids_model |> set_restrictions( labels = list(Y = c('00', '11')), given = 'X.00') ## ------------------------------------------------------------------------------------------------------------------------------------------------ #| label: get-priors #| echo: true #| eval: true lipids_model |> inspect("prior_hyperparameters", nodes = "X") ## ------------------------------------------------------------------------------------------------------------------------------------------------ #| label: set-priors #| echo: true #| purl: true model <- lipids_model |> set_priors(distribution = "jeffreys") ## ------------------------------------------------------------------------------------------------------------------------------------------------ #| label: set-priors-custom #| echo: true #| eval: true #| message: false lipids_model |> set_priors( param_names = c("X.10", "X.01"), alphas = 3:4) |> inspect("prior_hyperparameters", nodes = "X") ## ------------------------------------------------------------------------------------------------------------------------------------------------ #| label: set-priors-statement #| echo: true #| eval: true lipids_model |> set_priors( statement = "X[Z = 1] > X[Z = 0]", alphas = 3) |> inspect("prior_hyperparameters", nodes = "X") ## ------------------------------------------------------------------------------------------------------------------------------------------------ #| label: set-priors-flat query <- make_model("X -> Y") |> set_restrictions(decreasing("X", "Y")) |> query_model("Y[X = 1] - Y[X = 0]", using = "priors") ## ------------------------------------------------------------------------------------------------------------------------------------------------ #| label: get-parameters #| echo: true #| eval: true make_model("X -> Y") |> inspect("parameters") ## ------------------------------------------------------------------------------------------------------------------------------------------------ #| label: set-parameters #| echo: true #| eval: true make_model("X -> Y") |> set_parameters( statement = "Y[X = 1] > Y[X = 0]", parameters = .7) |> inspect("parameters") ## ------------------------------------------------------------------------------------------------------------------------------------------------ #| include: true #' ## Drawing and manipulating data ## ------------------------------------------------------------------------------------------------------------------------------------------------ #| label: make-data #| echo: true #| eval: true #| include: true lipids_model |> make_data(n = 4) ## ------------------------------------------------------------------------------------------------------------------------------------------------ #| label: make-data-incomplete #| echo: true #| eval: true #| message: false sample_data <- lipids_model |> make_data( n = 8, nodes = list(c("Z", "Y"), "X"), probs = list(1, .5), subsets = list(TRUE, "Z == 1 & Y == 0")) ## ------------------------------------------------------------------------------------------------------------------------------------------------ #| label: collapse-data #| echo: true #| eval: true sample_data |> collapse_data(lipids_model) ## ------------------------------------------------------------------------------------------------------------------------------------------------ #| include: true #' # Updating models ## ------------------------------------------------------------------------------------------------------------------------------------------------ #| include: true #' ## How the Stan program works ## ------------------------------------------------------------------------------------------------------------------------------------------------ #| echo: true #| eval: true #| purl: true #| include: true make_model("X -> Y") |> inspect("parameter_mapping") ## ------------------------------------------------------------------------------------------------------------------------------------------------ #| include: true #' ## Incomplete and censored data ## ------------------------------------------------------------------------------------------------------------------------------------------------ #| echo: true #| eval: true #| message: false data <- data.frame(X = rep(0:1, 5), Y = rep(0:1, 5)) list(uncensored = update_model(make_model("X -> Y"), data = data, refresh = 0), censored = update_model(make_model("X -> Y"), data = data, censored_types = c("X1Y0", "X0Y1"), refresh = 0)) |> query_model("Y[X = 1] - Y[X = 0]", using = "posteriors") |> dplyr::select(-using) ## ------------------------------------------------------------------------------------------------------------------------------------------------ #| include: true #' ## Output ## ------------------------------------------------------------------------------------------------------------------------------------------------ #| echo: true #| purl: true #| message: false model <- make_model("X -> Y") |> update_model() posterior <- inspect(model, "posterior_distribution") ## ------------------------------------------------------------------------------------------------------------------------------------------------ #| echo: true #| purl: true #| message: false lipids_model <- lipids_model |> update_model( keep_fit = TRUE, keep_event_probabilities = TRUE) ## ------------------------------------------------------------------------------------------------------------------------------------------------ #| eval: true #| echo: true #| purl: true #| message: false make_model("X -> Y") |> update_model(keep_type_distribution = FALSE, refresh = 0) |> inspect("stan_summary") ## ------------------------------------------------------------------------------------------------------------------------------------------------ #| include: true #' ## Convergence problems and diagnostics ## ------------------------------------------------------------------------------------------------------------------------------------------------ #| label: preservewwarning #| warning: false model <- make_model("X -> M -> Y") |> update_model( data = data.frame(X = rep(0:1, 10000), Y = rep(0:1, 10000)), iter = 5000, refresh = 0) ## ------------------------------------------------------------------------------------------------------------------------------------------------ #| eval: true #| echo: true #| purl: true model ## ------------------------------------------------------------------------------------------------------------------------------------------------ #| eval: true #| echo: true #| purl: true #| message: false model <- make_model("X -> Y") |> update_model(refresh = 0, keep_fit = TRUE) ## ------------------------------------------------------------------------------------------------------------------------------------------------ #| eval: true #| echo: true #| purl: true model |> inspect("stanfit") ## ------------------------------------------------------------------------------------------------------------------------------------------------ #| include: true #' # Querying models ## ------------------------------------------------------------------------------------------------------------------------------------------------ #| include: true #' ## Calculating factual and counterfactual quantities ## ------------------------------------------------------------------------------------------------------------------------------------------------ #| label: realise-outcomes #| echo: true #| include: true make_model("X -> Y") |> realise_outcomes() ## ------------------------------------------------------------------------------------------------------------------------------------------------ #| label: realise-outcomes-do #| echo: true #| include: true make_model("X -> Y") |> realise_outcomes(dos = list(X = 1)) ## ------------------------------------------------------------------------------------------------------------------------------------------------ #| include: true #' ## Causal syntax ## ------------------------------------------------------------------------------------------------------------------------------------------------ make_model("X -> Y") |> get_query_types("Y == 1") ## ------------------------------------------------------------------------------------------------------------------------------------------------ make_model("X -> Y") |> get_query_types("Y[X = 0] == 1") ## ------------------------------------------------------------------------------------------------------------------------------------------------ #| eval: true #| echo: true make_model("X1 -> Y <- X2") |> get_query_types( "X1 == 1 & X2 == 1 & (Y[X1 = 1, X2 = 1] > Y[X1 = 0, X2 = 0])") ## ------------------------------------------------------------------------------------------------------------------------------------------------ make_model("X -> Y") |> get_query_types("Y[X = 1] - Y[X = 0]") ## ------------------------------------------------------------------------------------------------------------------------------------------------ #| include: true #' ## Quantifying queries ## ------------------------------------------------------------------------------------------------------------------------------------------------ #| label: fig-posterior-dist #| fig-cap: 'Posterior on "Probability $Y$ is increasing in $X$".' #| fig-pos: "t" #| fig-align: center #| out-width: "60%" #| purl: true #| message: false data <- data.frame(X = rep(0:1, 50), Y = rep(0:1, 50)) model <- make_model("X -> Y") |> update_model( data, iter = 4000, refresh = 0) model |> grab("posterior_distribution") |> ggplot(aes(Y.01 - Y.10)) + geom_histogram() ## ------------------------------------------------------------------------------------------------------------------------------------------------ #| echo: true queries <- make_model("X -> Y") |> query_distribution( query = list(increasing = "(Y[X = 1] > Y[X = 0])", ATE = "(Y[X = 1] - Y[X = 0])"), using = "priors") ## ------------------------------------------------------------------------------------------------------------------------------------------------ #| label: fig-case-level-query #| fig-cap: "Illustration of case-level queries plotted." #| echo: true #| eval: true #| purl: true #| fig-width: 8 #| fig-height: 3 lipids_model |> query_model( query = "Y[X = 1] - Y[X = 0] :|: X == 1 & Y == 1 & Z == 1", using = "posteriors") |> plot() ## ------------------------------------------------------------------------------------------------------------------------------------------------ #| label: fig-new-case-level-query #| fig-cap: "Illustration of new case-level queries plotted." #| echo: true #| eval: true #| purl: true #| fig-width: 8 #| fig-height: 3 #| message: false make_model("X -> M -> Y") |> update_model( data.frame(X = rep(0:1, 8), Y = rep(0:1, 8)), iter = 4000, refresh = 0) |> query_model("Y[X = 1] > Y[X = 0] :|: X == 1 & Y == 1 & M == 1", using = "posteriors", case_level = c(TRUE, FALSE)) |> plot() ## ------------------------------------------------------------------------------------------------------------------------------------------------ #| label: batch-query #| echo: true #| eval: true #| message: false models <- list( Unrestricted = lipids_model |> update_model(data = lipids_data, refresh = 0), Restricted = lipids_model |> set_restrictions("X[Z = 1] < X[Z = 0]") |> update_model(data = lipids_data, refresh = 0) ) ## ------------------------------------------------------------------------------------------------------------------------------------------------ #| label: batch-query-fig #| echo: true #| eval: true queries <- query_model( models, query = list( ATE = "Y[X=1] - Y[X=0]", POS = "Y[X=1] > Y[X=0] :|: Y==1 & X==1"), case_level = c(FALSE, TRUE), using = c("priors", "posteriors"), expand_grid = TRUE) ## ------------------------------------------------------------------------------------------------------------------------------------------------ #| label: tbl-batch-query #| tbl-cap: "Results for two queries on two models." #| echo: true #| eval: true #| message: false queries |> dplyr::select(-starts_with("cred"), -sd) |> kabble(align = "c") ## ------------------------------------------------------------------------------------------------------------------------------------------------ #| include: true #' default plot associated with this query: ## ------------------------------------------------------------------------------------------------------------------------------------------------ #| label: fig-batch #| fig-cap: "Default plotting for a a set of queries over multiple models." #| fig-width: 8 #| fig-height: 4 #| echo: true #| eval: true #| message: false queries |> plot() ## ------------------------------------------------------------------------------------------------------------------------------------------------ #| include: true #' # Appendix A: Parallelization ## ------------------------------------------------------------------------------------------------------------------------------------------------ #| include: true library("parallel") options(mc.cores = parallel::detectCores()) ## ------------------------------------------------------------------------------------------------------------------------------------------------ #| echo: true #| include: true #| library("future") library("future.apply") chains <- 3 cores <- 8 future::plan(list( future::tweak(future::multisession, workers = floor(cores/(chains + 1))), future::tweak(future::multisession, workers = chains))) model <- make_model("X -> Y") data <- list(data_1 = data.frame(X = 0:1, Y = 0:1), data_2 = data.frame(X = 0:1, Y = 1:0)) results <- future.apply::future_lapply(data, function(d) { update_model(model = model, data = d, chains = chains, refresh = 0)}, future.seed = TRUE) ## ------------------------------------------------------------------------------------------------------------------------------------------------ #| include: true #' # Appendix B: Stan code ## ------------------------------------------------------------------------------------------------------------------------------------------------ #| echo: true #| include: true CausalQueries:::stanmodels$simplexes ## ------------------------------------------------------------------------------------------------------------------------------------------------ #| include: true #' # Appendix C: Benchmarks ## ------------------------------------------------------------------------------------------------------------------------------------------------ #| label: tbl-bench1 #| tbl-cap: "Benchmarking 1." #| echo: true #| eval: true #| warning: false #| message: false iter <- 4000 times <- 8 model <- list( CausalQueries::make_model("X -> Y"), CausalQueries::make_model("X1 -> Y <- X2"), CausalQueries::make_model("X1 -> Y; X2 -> Y; X3 -> Y") ) data <- expand_grid(X1 = 0:1, X2 = 0:1, X3 = 0:1) |> uncount(20) |> mutate(Y = 1 * (X1 - X2 * X3 + rnorm(n()) > 0)) benchmark_model <- microbenchmark::microbenchmark( m1 = CausalQueries::update_model(model[[1]], data, iter = iter, refresh = 0), m2 = CausalQueries::update_model(model[[2]], data, iter = iter, refresh = 0), m3 = CausalQueries::update_model(model[[3]], data, iter = iter, refresh = 0), times = times ) benchmark_1 <- data.frame( Model = c( "X1 -> Y", "X1 -> Y; X2 -> Y", "X1 -> Y; X2 -> Y; X3 -> Y" ), x = c(6, 20, 262), t = (benchmark_model |> summary() |> pull(mean))) benchmark_1 |> kabble(col.names = c("Model", "Number of parameters", "Runtime (seconds)"), align = "c") ## ------------------------------------------------------------------------------------------------------------------------------------------------ #| label: tbl-bench2 #| tbl-cap: "Benchmarking 2." #| echo: true #| message: false #| warning: false #| eval: true model <- CausalQueries::make_model("X -> Y") data <- lapply(10^c(1:4), function(n) { CausalQueries::make_data(model, n) }) benchmark_data <- microbenchmark::microbenchmark( d0 = CausalQueries::update_model(model, data[[1]], iter = iter, refresh = 0), d1 = CausalQueries::update_model(model, data[[2]], iter = iter, refresh = 0), d2 = CausalQueries::update_model(model, data[[3]], iter = iter, refresh = 0), d3 = CausalQueries::update_model(model, data[[4]], iter = iter, refresh = 0), times = times ) benchmark2 <- data.frame( Model = c("X1 -> Y", "X1 -> Y", "X1 -> Y", "X1 -> Y"), x = c(10, 100, 1000, 10000), t = (benchmark_data |> summary() |> pull(mean)) ) benchmark2 |> kabble(col.names = c("Model", "Number of observations", "Runtime (seconds)"), align = "c")