From 5de1447ce6a625d7010bd75ea1128539b2cec74e Mon Sep 17 00:00:00 2001 From: Daniel Date: Fri, 13 Sep 2024 14:03:33 +0200 Subject: [PATCH] fix --- R/check_dag.R | 63 ++++++++-- tests/testthat/_snaps/check_dag.md | 184 +++++++++++++++++++++++++++-- tests/testthat/test-check_dag.R | 28 +++++ 3 files changed, 258 insertions(+), 17 deletions(-) diff --git a/R/check_dag.R b/R/check_dag.R index 39dfc2b6c..6ca03a6f2 100644 --- a/R/check_dag.R +++ b/R/check_dag.R @@ -245,13 +245,33 @@ check_dag <- function(..., # data for checking effects checks <- lapply(c("direct", "total"), function(x) { - adjustment_set <- unlist(dagitty::adjustmentSets(dag, effect = x), use.names = FALSE) + if(is.null(adjusted)) { + adjustment_set <- unlist(dagitty::adjustmentSets(dag, effect = x), use.names = FALSE) + } else { + adjustment_set <- adjusted + } adjustment_nodes <- unlist(dagitty::adjustedNodes(dag), use.names = FALSE) + minimal_adjustments <- as.list(dagitty::adjustmentSets(dag, effect = x)) + collider <- adjustment_nodes[vapply(adjustment_nodes, ggdag::is_collider, logical(1), .dag = dag)] + if (!length(collider)) { + # if we don't have colliders, set to NULL + collider <- NULL + } else { + # if we *have* colliders, remove them from minimal adjustments + minimal_adjustments <- lapply(minimal_adjustments, function(ma) { + setdiff(ma, collider) + }) + } list( - adjustment_not_needed = is.null(adjustment_set) && is.null(adjustment_nodes), - incorrectly_adjusted = is.null(adjustment_set) && !is.null(adjustment_nodes), + adjustment_not_needed = (is.null(adjustment_set) && is.null(adjustment_nodes) || + (identical(sort(unique(adjustment_set)), sort(unique(setdiff(c(exposure, adjustment_nodes), collider)))))) && + is.null(collider), + incorrectly_adjusted = (is.null(adjustment_set) && !is.null(adjustment_nodes)) || + (!identical(sort(unique(adjustment_set)), sort(unique(setdiff(c(exposure, adjustment_nodes), collider))))) || + (!is.null(collider) && collider %in% adjustment_nodes), current_adjustments = adjustment_nodes, - minimal_adjustments = as.list(dagitty::adjustmentSets(dag, effect = x)) + minimal_adjustments = minimal_adjustments, + collider = collider ) }) @@ -260,6 +280,10 @@ check_dag <- function(..., attr(dag, "exposure") <- exposure attr(dag, "adjusted") <- adjusted attr(dag, "adjustment_sets") <- checks[[1]]$current_adjustments + attr(dag, "collider") <- checks[[1]]$collider + # remove collider from sub-attributes + checks[[1]]$collider <- NULL + checks[[2]]$collider <- NULL attr(dag, "check_direct") <- insight::compact_list(checks[[1]]) attr(dag, "check_total") <- insight::compact_list(checks[[2]]) @@ -296,6 +320,7 @@ as.dag <- function(x, ...) { #' @export print.check_dag <- function(x, ...) { effect <- attributes(x)$effect + collider <- attributes(x)$collider # header cat(insight::print_color("# Check for correct adjustment sets", "blue")) @@ -317,6 +342,16 @@ print.check_dag <- function(x, ...) { ) } + # add information on colliders + if (!is.null(collider)) { + exposure_outcome_text <- paste0( + exposure_outcome_text, + "\n- Collider", + ifelse(length(collider) > 1, "s", ""), + ": ", insight::color_text(datawizard::text_concatenate(collider), "cyan") + ) + } + cat(exposure_outcome_text) cat("\n\n") @@ -331,12 +366,12 @@ print.check_dag <- function(x, ...) { } else { out <- attributes(x)$check_total } - .print_dag_results(out, x, i, effect) + .print_dag_results(out, x, i, effect, collider) } } } -.print_dag_results <- function(out, x, i, effect) { +.print_dag_results <- function(out, x, i, effect, collider = NULL) { # missing adjustements - minimal_adjustment can be a list of different # options for minimal adjustements, so we check here if any of the minimal # adjustments are currently sufficient @@ -356,8 +391,18 @@ print.check_dag <- function(x, ...) { attributes(x)$outcome, "`." ) + } else if (!is.null(collider)) { + # Scenario 2: adjusted for (downstream) collider + msg <- paste0( + insight::color_text("Incorrectly adjusted!", "red"), + "\nYour model adjusts for a (downstream) collider, ", + insight::color_text(datawizard::text_concatenate(collider, enclose = "`"), "cyan"), + ". To estimate the ", i, " effect, do ", + insight::color_text("not", "italic"), + " adjust for it, to avoid collider-bias." + ) } else if (isTRUE(out$incorrectly_adjusted)) { - # Scenario 2: incorrectly adjusted, adjustments where none is allowed + # Scenario 3: incorrectly adjusted, adjustments where none is allowed msg <- paste0( insight::color_text("Incorrectly adjusted!", "red"), "\nTo estimate the ", i, " effect, do ", @@ -367,13 +412,13 @@ print.check_dag <- function(x, ...) { "." ) } else if (any(sufficient_adjustments)) { - # Scenario 3: correct adjustment + # Scenario 4: correct adjustment msg <- paste0( insight::color_text("Model is correctly specified.", "green"), "\nAll minimal sufficient adjustments to estimate the ", i, " effect were done." ) } else { - # Scenario 4: missing adjustments + # Scenario 5: missing adjustments msg <- paste0( insight::color_text("Incorrectly adjusted!", "red"), "\nTo estimate the ", i, " effect, ", diff --git a/tests/testthat/_snaps/check_dag.md b/tests/testthat/_snaps/check_dag.md index aa3ebc45e..a63da9ea0 100644 --- a/tests/testthat/_snaps/check_dag.md +++ b/tests/testthat/_snaps/check_dag.md @@ -1,3 +1,128 @@ +# check_dag + + Code + print(dag) + Output + # Check for correct adjustment sets + - Outcome: y + - Exposure: x + + Identification of direct and total effects + + Model is correctly specified. + No adjustment needed to estimate the direct and total effect of `x` on `y`. + + +--- + + Code + print(dag) + Output + # Check for correct adjustment sets + - Outcome: y + - Exposure: x + - Adjustment: b + + Identification of direct and total effects + + Incorrectly adjusted! + To estimate the direct and total effect, do not adjust for `b`. + + +--- + + Code + print(dag) + Output + # Check for correct adjustment sets + - Outcome: y + - Exposure: x + + Identification of direct and total effects + + Incorrectly adjusted! + To estimate the direct and total effect, do not adjust for . + + +--- + + Code + print(dag) + Output + # Check for correct adjustment sets + - Outcome: y + - Exposure: x + - Adjustment: c + + Identification of direct and total effects + + Incorrectly adjusted! + To estimate the direct and total effect, do not adjust for `c`. + + +--- + + Code + print(dag) + Output + # Check for correct adjustment sets + - Outcome: y + - Exposure: x + - Adjustment: c + + Identification of direct and total effects + + Incorrectly adjusted! + To estimate the direct and total effect, do not adjust for `c`. + + +--- + + Code + print(dag) + Output + # Check for correct adjustment sets + - Outcome: mpg + - Exposure: wt + - Adjustments: cyl, disp and gear + + Identification of direct and total effects + + Model is correctly specified. + No adjustment needed to estimate the direct and total effect of `wt` on `mpg`. + + +# check_dag, multiple adjustment sets + + Code + print(dag) + Output + # Check for correct adjustment sets + - Outcome: exam + - Exposure: podcast + + Identification of direct and total effects + + Incorrectly adjusted! + To estimate the direct and total effect, do not adjust for . + + +--- + + Code + print(dag) + Output + # Check for correct adjustment sets + - Outcome: exam + - Exposure: podcast + - Adjustments: alertness and prepared + + Identification of direct and total effects + + Incorrectly adjusted! + To estimate the direct and total effect, do not adjust for `alertness` and `prepared`. + + # check_dag, different adjustements for total and direct Code @@ -10,12 +135,12 @@ Identification of direct effects Incorrectly adjusted! - To estimate the direct effect, at least adjust for `x1` and `x2`. Currently, the model does not adjust for any variables. + To estimate the direct effect, do not adjust for . Identification of total effects Incorrectly adjusted! - To estimate the total effect, at least adjust for `x1`. Currently, the model does not adjust for any variables. + To estimate the total effect, do not adjust for . --- @@ -31,12 +156,12 @@ Identification of direct effects Incorrectly adjusted! - To estimate the direct effect, at least adjust for `x1` and `x2`. Currently, the model only adjusts for `x1`. You possibly also need to adjust for `x2` to block biasing paths. + To estimate the direct effect, do not adjust for `x1`. Identification of total effects - Model is correctly specified. - All minimal sufficient adjustments to estimate the total effect were done. + Incorrectly adjusted! + To estimate the total effect, do not adjust for `x1`. --- @@ -52,7 +177,7 @@ Identification of direct effects Incorrectly adjusted! - To estimate the direct effect, at least adjust for `x1` and `x2`. Currently, the model only adjusts for `x2`. You possibly also need to adjust for `x1` to block biasing paths. + To estimate the direct effect, do not adjust for `x2`. Identification of total effects @@ -72,8 +197,8 @@ Identification of direct effects - Model is correctly specified. - All minimal sufficient adjustments to estimate the direct effect were done. + Incorrectly adjusted! + To estimate the direct effect, do not adjust for `x1` and `x2`. Identification of total effects @@ -81,3 +206,46 @@ To estimate the total effect, do not adjust for `x1` and `x2`. +# check_dag, collider bias + + Code + print(dag) + Output + # Check for correct adjustment sets + - Outcome: SMD_ICD11 + - Exposure: agegroup + - Adjustments: edgroup3, gender_kid, pss4_kid_sum_2sd and residence + + Identification of direct effects + + Model is correctly specified. + No adjustment needed to estimate the direct effect of `agegroup` on `SMD_ICD11`. + + Identification of total effects + + Model is correctly specified. + No adjustment needed to estimate the total effect of `agegroup` on `SMD_ICD11`. + + +--- + + Code + print(dag) + Output + # Check for correct adjustment sets + - Outcome: SMD_ICD11 + - Exposure: agegroup + - Adjustments: edgroup3, gender_kid, pss4_kid_sum_2sd, residence and sm_h_total_kid + - Collider: sm_h_total_kid + + Identification of direct effects + + Incorrectly adjusted! + Your model adjusts for a (downstream) collider, `sm_h_total_kid`. To estimate the direct effect, do not adjust for it, to avoid collider-bias. + + Identification of total effects + + Incorrectly adjusted! + Your model adjusts for a (downstream) collider, `sm_h_total_kid`. To estimate the total effect, do not adjust for it, to avoid collider-bias. + + diff --git a/tests/testthat/test-check_dag.R b/tests/testthat/test-check_dag.R index 1ac94d520..7bde379f1 100644 --- a/tests/testthat/test-check_dag.R +++ b/tests/testthat/test-check_dag.R @@ -133,3 +133,31 @@ test_that("check_dag, different adjustements for total and direct", { ) expect_snapshot(print(dag)) }) + +test_that("check_dag, collider bias", { + dag <- check_dag( + SMD_ICD11 ~ agegroup + gender_kid + edgroup3 + residence + pss4_kid_sum_2sd + sm_h_total_kid, + pss4_kid_sum_2sd ~ gender_kid, + sm_h_total_kid ~ gender_kid + agegroup, + adjusted = c( + "agegroup", "gender_kid", "edgroup3", "residence", + "pss4_kid_sum_2sd" + ), + outcome = "SMD_ICD11", + exposure = "agegroup" + ) + expect_snapshot(print(dag)) + + dag <- check_dag( + SMD_ICD11 ~ agegroup + gender_kid + edgroup3 + residence + pss4_kid_sum_2sd + sm_h_total_kid, + pss4_kid_sum_2sd ~ gender_kid, + sm_h_total_kid ~ gender_kid + agegroup, + adjusted = c( + "agegroup", "gender_kid", "edgroup3", "residence", + "pss4_kid_sum_2sd", "sm_h_total_kid" + ), + outcome = "SMD_ICD11", + exposure = "agegroup" + ) + expect_snapshot(print(dag)) +})