Overview

The R package R6causal implements an R6 class called SCM. The class aims to simplify working with structural causal models. The missing data mechanism can be defined as a part of the structural model.

The class contains methods for

defining a structural causal model via functions, text or conditional probability tables
printing basic information on the model
plotting the graph for the model using packages igraph or qgraph
simulating data from the model
applying an intervention
checking the identifiability of a query using the R packages causaleffect and dosearch
defining the missing data mechanism
simulating incomplete data from the model according to the specified missing data mechanism
checking the identifiability in a missing data problem using the R package dosearch
checking the identifiability of a counterfactual query using the R package cfid

In addition, there are functions for

running experiments
counterfactual inference using simulation
evaluating fairness of a prediction model

The class ParallelWorld inherits SCM and defines a structural causal model that describes parallel worlds for counterfactual inference.

The class LinearGaussianSCM inherits SCM and defines a structural causal model where all functions are linear and all background variables follow Gaussian distribution.

Setup

library(R6causal)
library(data.table)
#> 
#> Attaching package: 'data.table'
#> The following object is masked from 'package:base':
#> 
#>     %notin%
library(stats)
data.table::setDTthreads(2)

Defining the model

Structural causal model (SCM) for a backdoor situation can be defined as follows

backdoor <- SCM$new("backdoor",
  uflist = list(
    uz = function(n) {return(runif(n))},
    ux = function(n) {return(runif(n))},
    uy = function(n) {return(runif(n))}
  ),
  vflist = list(
    z = function(uz) {
      return(as.numeric(uz < 0.4))},
    x = function(ux, z) {
      return(as.numeric(ux < 0.2 + 0.5*z))},
    y = function(uy, z, x) {
      return(as.numeric(uy < 0.1 + 0.4*z + 0.4*x))}
  )
)

A shortcut notation for this is

backdoor_text <- SCM$new("backdoor",
  uflist = list(
    uz = "n : runif(n)", 
    ux = "n : runif(n)", 
    uy = "n : runif(n)" 
  ),
  vflist = list(
    z = "uz : as.numeric(uz < 0.4)", 
    x = "ux, z : as.numeric(ux < 0.2 + 0.5*z)",
    y = "uy, z, x : as.numeric(uy < 0.1 + 0.4*z + 0.4*x)" 
  )
)

Alternatively the functions of SCM can be specified via conditional probability tables

backdoor_condprob <- SCM$new("backdoor",
  uflist = list(
    uz = function(n) {return(runif(n))},
    ux = function(n) {return(runif(n))},
    uy = function(n) {return(runif(n))}
  ),
  vflist = list(
    z = function(uz) {
      return( generate_condprob( ycondx = data.table(z = c(0,1), 
                                                     prob = c(0.6,0.4)), 
                               x = data.table(uz = uz), 
                               Umerge_expr = "uz"))},
    x = function(ux, z) {
      return( generate_condprob( ycondx = data.table(x = c(0,1,0,1), 
                                                     z = c(0,0,1,1), 
                                                     prob = c(0.8,0.2,0.3,0.7)),
                                             x = data.table(z = z, ux = ux), 
                                             Umerge_expr = "ux"))},
    y = function(uy, z, x) {
      return( generate_condprob( ycondx = data.table(y= rep(c(0,1), 4), 
                                                     z = c(0,0,1,1,0,0,1,1), 
                                                     x = c(0,0,0,0,1,1,1,1), 
                                                     prob = c(0.9,0.1,0.5,0.5,
                                                              0.5,0.5,0.1,0.9)),
                                             x = data.table(z = z, x = x, uy = uy), 
                                             Umerge_expr = "uy"))}
  )
)

It is possible to mix the styles and define some elements of a function list as functions, some as text and some as conditional probability tables.

Defining a linear Gaussian SCM

A linear Gaussian SCM can be defined giving the coefficients for the structural equations:

lgbackdoor <- LinearGaussianSCM$new("Linear Gaussian Backdoor",
                                    linear_gaussian = list(
                                      uflist = list(ux = function(n) {rnorm(n)},
                                                    uy = function(n) {rnorm(n)},
                                                    uz = function(n) {rnorm(n)}),
                                      vnames = c("x","y","z"),
                                      vcoefmatrix = matrix(c(0,0.4,0,0,0,0,0.6,0.8,0),3,3),
                                      ucoefvector = c(1,1,1),
                                      ccoefvector = c(0,0,0)))
print(lgbackdoor)
#> Name of the model:  Linear Gaussian Backdoor 
#> 
#> Graph: 
#>  z -> x 
#>  x -> y 
#>  z -> y 
#>   
#> Functions of background (exogenous) variables: 
#> 
#> $ux
#> function (n) 
#> {
#>     rnorm(n)
#> }
#> 
#> $uy
#> function (n) 
#> {
#>     rnorm(n)
#> }
#> 
#> $uz
#> function (n) 
#> {
#>     rnorm(n)
#> }
#> 
#> Functions of endogenous variables: 
#> 
#> $x
#> function (z, ux) 
#> {
#>     return(0 + 0.6 * z + 1 * ux)
#> }
#> <environment: 0x59f027fdba48>
#> 
#> $y
#> function (x, z, uy) 
#> {
#>     return(0 + 0.4 * x + 0.8 * z + 1 * uy)
#> }
#> <environment: 0x59f027fe0928>
#> 
#> $z
#> function (uz) 
#> {
#>     return(0 + 1 * uz)
#> }
#> <environment: 0x59f027fe1850>
#> 
#> Topological order of endogenous variables: 
#> [1] "z" "x" "y"
#> 
#> No missing data mechanism

It is also possible to generate the underlying DAG and the coefficients randomly:

randomlg <- LinearGaussianSCM$new("Random Linear Gaussian",
                                  random_linear_gaussian = list(
                                  nv = 6, 
                                  edgeprob=0.5, 
                                  vcoefdistr = function(n) {rnorm(n)}, 
                                  ccoefdistr = function(n) {rnorm(n)}, 
                                  ucoefdistr = function(n) {rnorm(n)}))
print(randomlg)
#> Name of the model:  Random Linear Gaussian 
#> 
#> Graph: 
#>  v3 -> v1 
#>  v3 -> v2 
#>  v1 -> v4 
#>  v2 -> v4 
#>  v6 -> v4 
#>  v3 -> v6 
#>  v5 -> v6 
#>   
#> Functions of background (exogenous) variables: 
#> 
#> $u1
#> function (n) 
#> {
#>     return(rnorm(n))
#> }
#> <environment: 0x59f028e24ee8>
#> 
#> $u2
#> function (n) 
#> {
#>     return(rnorm(n))
#> }
#> <environment: 0x59f028e263c0>
#> 
#> $u3
#> function (n) 
#> {
#>     return(rnorm(n))
#> }
#> <environment: 0x59f028e280f8>
#> 
#> $u4
#> function (n) 
#> {
#>     return(rnorm(n))
#> }
#> <environment: 0x59f028e2d400>
#> 
#> $u5
#> function (n) 
#> {
#>     return(rnorm(n))
#> }
#> <environment: 0x59f026763c10>
#> 
#> $u6
#> function (n) 
#> {
#>     return(rnorm(n))
#> }
#> <environment: 0x59f02675fa60>
#> 
#> Functions of endogenous variables: 
#> 
#> $v1
#> function (v3, u1) 
#> {
#>     return(0.994876274613695 + 1.01902826651745 * v3 + 0.341784654299718 * 
#>         u1)
#> }
#> <environment: 0x59f026757668>
#> 
#> $v2
#> function (v3, u2) 
#> {
#>     return(0.45408326729864 + 1.51327081989653 * v3 + -0.574519867807911 * 
#>         u2)
#> }
#> <environment: 0x59f026756468>
#> 
#> $v3
#> function (u3) 
#> {
#>     return(1.15521996304971 + -0.642014013791146 * u3)
#> }
#> <environment: 0x59f02674d4a0>
#> 
#> $v4
#> function (v1, v2, v6, u4) 
#> {
#>     return(-1.27321583955983 + 0.62242635086953 * v1 + 1.43222002640366 * 
#>         v2 + 0.354128536029517 * v6 + -1.04268922489232 * u4)
#> }
#> <environment: 0x59f02674c7a8>
#> 
#> $v5
#> function (u5) 
#> {
#>     return(-0.11736079042141 + -0.382778548897075 * u5)
#> }
#> <environment: 0x59f026745f20>
#> 
#> $v6
#> function (v3, v5, u6) 
#> {
#>     return(0.0679395989185808 + 0.874975193771324 * v3 + -0.392650895955235 * 
#>         v5 + 0.137286996420426 * u6)
#> }
#> <environment: 0x59f026738218>
#> 
#> Topological order of endogenous variables: 
#> [1] "v3" "v5" "v1" "v2" "v6" "v4"
#> 
#> No missing data mechanism

Printing the model

The print method presents the basic information on the model

backdoor
#> Name of the model:  backdoor 
#> 
#> Graph: 
#>  z -> x 
#>  z -> y 
#>  x -> y 
#>   
#> Functions of background (exogenous) variables: 
#> 
#> $uz
#> function (n) 
#> {
#>     return(runif(n))
#> }
#> 
#> $ux
#> function (n) 
#> {
#>     return(runif(n))
#> }
#> 
#> $uy
#> function (n) 
#> {
#>     return(runif(n))
#> }
#> 
#> Functions of endogenous variables: 
#> 
#> $z
#> function (uz) 
#> {
#>     return(as.numeric(uz < 0.4))
#> }
#> 
#> $x
#> function (ux, z) 
#> {
#>     return(as.numeric(ux < 0.2 + 0.5 * z))
#> }
#> 
#> $y
#> function (uy, z, x) 
#> {
#>     return(as.numeric(uy < 0.1 + 0.4 * z + 0.4 * x))
#> }
#> 
#> Topological order of endogenous variables: 
#> [1] "z" "x" "y"
#> 
#> No missing data mechanism

Plotting the graph

The plotting method of the package igraph is used by default. If qgraph is available, its plotting method can be used as well. The argument subset controls which variables are plotted. Plotting parameters are passed to the plotting method.

backdoor$plot(vertex.size = 25) # with package 'igraph'

backdoor$plot(subset = "v") # only observed variables

if (requireNamespace("qgraph", quietly = TRUE)) backdoor$plot(method = "qgraph")

# alternative look with package 'qgraph'

Simulating data

Calling method simulate() creates or updates data table simdata.

backdoor$simulate(10)
backdoor$simdata
#>             uz         ux         uy     z     x     y
#>          <num>      <num>      <num> <num> <num> <num>
#>  1: 0.60541397 0.44072661 0.56770626     0     0     0
#>  2: 0.25877978 0.24649332 0.28064333     1     1     1
#>  3: 0.30488493 0.98923513 0.27242901     1     0     1
#>  4: 0.06904353 0.48486931 0.92127339     1     1     0
#>  5: 0.41610374 0.35687211 0.30543364     0     0     0
#>  6: 0.01150970 0.34392137 0.93502956     1     1     0
#>  7: 0.72304303 0.54422719 0.50026632     0     0     0
#>  8: 0.46653116 0.01560096 0.08414311     0     1     1
#>  9: 0.82308997 0.98123438 0.13427670     0     0     0
#> 10: 0.36184499 0.90960819 0.04239095     1     0     1
backdoor$simulate(8)
backdoor$simdata
#>           uz        ux        uy     z     x     y
#>        <num>     <num>     <num> <num> <num> <num>
#> 1: 0.5303116 0.9748824 0.3782929     0     0     0
#> 2: 0.2400277 0.4213903 0.3874816     1     1     1
#> 3: 0.7991782 0.1931011 0.5702527     0     1     0
#> 4: 0.1221450 0.7062660 0.8546533     1     0     0
#> 5: 0.6070940 0.5582551 0.6224968     0     0     0
#> 6: 0.5806234 0.5181475 0.1343228     0     0     0
#> 7: 0.9369502 0.6242546 0.7375925     0     0     0
#> 8: 0.1676231 0.5420063 0.5282841     1     1     1
backdoor_text$simulate(20)
backdoor_condprob$simulate(30)

Applying an intervention

In an intervention, the structural equation of the target variable is changed.

backdoor_x1 <- backdoor$clone()  # making a copy
backdoor_x1$intervene("x",1) # applying the intervention
backdoor_x1$plot(method = "qgraph") # to see that arrows incoming to x are cut

backdoor_x1$simulate(10) # simulating from the intervened model
backdoor_x1$simdata
#>              uz          ux        uy     z     x     y
#>           <num>       <num>     <num> <num> <num> <num>
#>  1: 0.549629843 0.981114544 0.5802219     0     1     0
#>  2: 0.009270503 0.581465549 0.4700830     1     1     1
#>  3: 0.537568542 0.815849268 0.5933757     0     1     0
#>  4: 0.172508582 0.083774623 0.1140479     1     1     1
#>  5: 0.434797424 0.001022082 0.3589887     0     1     1
#>  6: 0.874668694 0.534036309 0.4094776     0     1     1
#>  7: 0.061445728 0.823698662 0.1278462     1     1     1
#>  8: 0.691796059 0.201262525 0.8082143     0     1     0
#>  9: 0.880458870 0.317214495 0.5830123     0     1     0
#> 10: 0.720024795 0.573694530 0.8709501     0     1     0

An intervention can redefine a structural equation

backdoor_yz <- backdoor$clone()  # making a copy
backdoor_yz$intervene("y", 
  function(uy, z) {return(as.numeric(uy < 0.1 + 0.8*z ))}) # making y a function of z only
backdoor_yz$plot(method = "qgraph") # to see that arrow x -> y is cut

Running an experiment (set of interventions)

The function run_experiment applies a set of interventions, simulates data and collects the results.

backdoor_experiment <- run_experiment(backdoor, 
                                      intervene = list(x = c(0,1)), 
                                      response = "y", 
                                      n = 10000)
str(backdoor_experiment)
#> List of 2
#>  $ interventions:Classes 'data.table' and 'data.frame':  2 obs. of  1 variable:
#>   ..$ x: num [1:2] 0 1
#>   ..- attr(*, ".internal.selfref")=<pointer: 0x59f024256bd0> 
#>   ..- attr(*, "sorted")= chr "x"
#>  $ response_list:List of 1
#>   ..$ y:Classes 'data.table' and 'data.frame':   10000 obs. of  2 variables:
#>   .. ..$ V1: num [1:10000] 0 0 0 0 0 1 0 0 0 0 ...
#>   .. ..$ V2: num [1:10000] 1 0 1 1 1 0 1 1 1 0 ...
#>   .. ..- attr(*, ".internal.selfref")=<pointer: 0x59f024256bd0>
colMeans(backdoor_experiment$response_list$y)
#>     V1     V2 
#> 0.2490 0.6552

Applying the ID algorithm, Do-search and cfid

There are direct plugins to R packages causaleffect, dosearch and cfid that can be used to solve identifiability problems.

backdoor$causal.effect(y = "y", x = "x")
#> [1] "\\sum_{z}P(y|z,x)P(z)"
backdoor$dosearch(data = "p(x,y,z)", query = "p(y|do(x))")
#> \sum_{z}\left(p(z)p(y|x,z)\right)
backdoor$cfid(gamma = cfid::conj(cfid::cf("Y",0), cfid::cf("X",0, c(Z=1))) ) 
#> The query P(y /\ x_{z'}) is not identifiable from P_*.

Counterfactual inference (a simple case)

Let us assume that intervention do(X=0) was applied and the response Y = 0 was recorded. What is the probability that in this situation the intervention do(X=1) would have led to the response Y = 1? We estimate this probability by means of simulation.

cfdata <- counterfactual(backdoor, situation = list(do = list(target = "x", ifunction = 0), 
                                                    condition = data.table( x = 0, y = 0)), 
                         target = "x", ifunction = 1, n = 100000, 
                         method = "rejection")
mean(cfdata$y)
#> [1] 0.53698

The result differs from P(Y = 1 | do(X = 1))

backdoor_x1$simulate(100000)
mean(backdoor_x1$simdata$y)
#> [1] 0.65883

Counterfactual inference (parallel worlds)

Parallel world graphs (a generalization of a twin graph) are used for counterfactual inference with several counterfactual interventions . The package implements class ParallelWorld which heritates class SCM. A ParallelWorld object is created from an SCM object by specifying the interventions for each world. By default the variables of the parallel worlds are named with suffixes “_1”, “_2”, …

In the example below, we have the original world (variables x, z, y) and its two variants. In the variant 1 (variables x_1, z_1, y_1), the value of x (variable x_1 in the object) is set to be 0. In the variant 2 (variables x_2, z_2, y_2), the value of x (variable x_2 in the object) is set to be 0 and the value of z (variable z_2 in the object) is set to be 1.

backdoor_parallel <- ParallelWorld$new(
                         backdoor,
                        dolist=list(
                             list(target = "x", 
                                 ifunction = 0),
                             list(target = list("z","x"), 
                                  ifunction = list(1,0))
                         )
 )
backdoor_parallel
#> Name of the model:  backdoor 
#> 
#> Graph: 
#>  uz -> z 
#>  z -> x 
#>  uy -> y 
#>  z -> y 
#>  x -> y 
#>  uz -> z_1 
#>  uy -> y_1 
#>  z_1 -> y_1 
#>  x_1 -> y_1 
#>  uy -> y_2 
#>  z_2 -> y_2 
#>  x_2 -> y_2 
#>   
#> Functions of background (exogenous) variables: 
#> 
#> $uz
#> function (n) 
#> {
#>     return(runif(n))
#> }
#> <bytecode: 0x59f02fe8a170>
#> 
#> $ux
#> function (n) 
#> {
#>     return(runif(n))
#> }
#> <bytecode: 0x59f02fe8a170>
#> 
#> $uy
#> function (n) 
#> {
#>     return(runif(n))
#> }
#> <bytecode: 0x59f02fe8a170>
#> 
#> Functions of endogenous variables: 
#> 
#> $z
#> function (uz) 
#> {
#>     return(as.numeric(uz < 0.4))
#> }
#> <bytecode: 0x59f02ffc21e0>
#> 
#> $x
#> function (ux, z) 
#> {
#>     return(as.numeric(ux < 0.2 + 0.5 * z))
#> }
#> <bytecode: 0x59f03033e110>
#> 
#> $y
#> function (uy, z, x) 
#> {
#>     return(as.numeric(uy < 0.1 + 0.4 * z + 0.4 * x))
#> }
#> <bytecode: 0x59f030641900>
#> 
#> $z_1
#> function (uz) 
#> {
#>     return(as.numeric(uz < 0.4))
#> }
#> 
#> $x_1
#> function (...) 
#> {
#>     return(constant)
#> }
#> <environment: 0x59f034f7b6c0>
#> 
#> $y_1
#> function (uy, z_1, x_1) 
#> {
#>     return(as.numeric(uy < 0.1 + 0.4 * z_1 + 0.4 * x_1))
#> }
#> 
#> $z_2
#> function (...) 
#> {
#>     return(constant)
#> }
#> <environment: 0x59f034b96938>
#> 
#> $x_2
#> function (...) 
#> {
#>     return(constant)
#> }
#> <environment: 0x59f034b93d68>
#> 
#> $y_2
#> function (uy, z_2, x_2) 
#> {
#>     return(as.numeric(uy < 0.1 + 0.4 * z_2 + 0.4 * x_2))
#> }
#> 
#> Topological order of endogenous variables: 
#> [1] "x_1" "z_2" "x_2" "z"   "z_1" "y_2" "x"   "y_1" "y"  
#> 
#> No missing data mechanism
if (requireNamespace("qgraph", quietly = TRUE)) backdoor_parallel$plot(method = "qgraph")

Counterfactual data can be simulated with function counterfactual. In the example below, we know that variable y obtained value 0 in the original world as well as variants 1 and 2. We are interested in the counterfactual distribution of y if x had been set to 1.

cfdata <- counterfactual(backdoor_parallel,
                         situation = list(
                            do = NULL,
                            condition = data.table::data.table( y = 0, y_1 = 0, y_2 = 0)),
                         target = "x",
                         ifunction = 1,
                         n = 100000,
                         method = "rejection")
mean(cfdata$y)
#> [1] 0.12403

The printed value is a simulation based estimate for the counterfactual probability $P(Y=1)$.

An alternative way for answering the same question defines the case of interest as one of the parallel worlds (here variant 3).

backdoor_parallel2 <- ParallelWorld$new(
                         backdoor,
                        dolist=list(
                             list(target = "x", 
                                 ifunction = 0),
                             list(target = list("z","x"), 
                                  ifunction = list(1,0)),
                              list(target = "x", 
                                 ifunction = 1)
                         )
 )
cfdata <- counterfactual(backdoor_parallel2,
                         situation = list(
                            do = NULL,
                            condition = data.table::data.table( y = 0, y_1 = 0, y_2 = 0)),
                         n = 100000, 
                         method = "rejection")
mean(cfdata$y_3)
#> [1] 0.12382

The printed value is a simulation based estimate for the counterfactual probability $P(Y=1)$.

A model with a missing data mechanism

The missing data mechanism is defined in similar manner as the other variables.

backdoor_md <- SCM$new("backdoor_md",
                       uflist = list(
                         uz = "n : runif(n)",
                         ux = "n : runif(n)",
                         uy = "n : runif(n)",
                         urz = "n : runif(n)",
                         urx = "n : runif(n)",
                         ury = "n : runif(n)"
                       ),
                       vflist = list(
                         z = "uz : as.numeric(uz < 0.4)",
                         x = "ux, z : as.numeric(ux < 0.2 + 0.5*z)",
                         y = "uy, z, x : as.numeric(uy < 0.1 + 0.4*z + 0.4*x)"
                       ),
                       rflist = list(
                         z = "urz : as.numeric( urz < 0.9)",
                         x = "urx, z : as.numeric( (urx + z)/2 < 0.9)",
                         y = "ury, z : as.numeric( (ury + z)/2 < 0.9)"
                       ),
                       rprefix = "r_"
)

Plotting the graph for a model with missing data mechanism

backdoor_md$plot(vertex.size = 25, edge.arrow.size=0.5) # with package 'igraph'

backdoor_md$plot(subset = "v") # only observed variables a

if (!requireNamespace("qgraph", quietly = TRUE)) backdoor_md$plot(method = "qgraph") 
# alternative look with package 'qgraph'

Simulating incomplete data

By default both complete data and incomplete data are simulated. The incomplete dataset is named as $simdata_obs.

backdoor_md$simulate(100)
summary(backdoor_md$simdata)
#>        uz                 ux                 uy               urz          
#>  Min.   :0.002313   Min.   :0.003122   Min.   :0.01879   Min.   :0.003848  
#>  1st Qu.:0.285032   1st Qu.:0.176956   1st Qu.:0.26043   1st Qu.:0.269699  
#>  Median :0.483755   Median :0.474820   Median :0.47428   Median :0.510355  
#>  Mean   :0.498774   Mean   :0.470626   Mean   :0.47280   Mean   :0.479431  
#>  3rd Qu.:0.714295   3rd Qu.:0.738180   3rd Qu.:0.65695   3rd Qu.:0.693502  
#>  Max.   :0.998902   Max.   :0.989408   Max.   :0.99489   Max.   :0.995606  
#>                                                                            
#>       urx               ury                z              x       
#>  Min.   :0.02073   Min.   :0.02312   Min.   :0.00   Min.   :0.00  
#>  1st Qu.:0.25458   1st Qu.:0.21659   1st Qu.:0.00   1st Qu.:0.00  
#>  Median :0.50155   Median :0.53516   Median :0.00   Median :0.00  
#>  Mean   :0.48424   Mean   :0.51389   Mean   :0.38   Mean   :0.43  
#>  3rd Qu.:0.64228   3rd Qu.:0.76593   3rd Qu.:1.00   3rd Qu.:1.00  
#>  Max.   :0.99214   Max.   :0.99116   Max.   :1.00   Max.   :1.00  
#>                                                                   
#>        y             z_md             x_md             y_md       
#>  Min.   :0.00   Min.   :0.0000   Min.   :0.0000   Min.   :0.0000  
#>  1st Qu.:0.00   1st Qu.:0.0000   1st Qu.:0.0000   1st Qu.:0.0000  
#>  Median :0.00   Median :0.0000   Median :0.0000   Median :0.0000  
#>  Mean   :0.47   Mean   :0.4043   Mean   :0.4149   Mean   :0.4157  
#>  3rd Qu.:1.00   3rd Qu.:1.0000   3rd Qu.:1.0000   3rd Qu.:1.0000  
#>  Max.   :1.00   Max.   :1.0000   Max.   :1.0000   Max.   :1.0000  
#>                 NAs    :6        NAs    :6        NAs    :11      
#>       r_z            r_x            r_y      
#>  Min.   :0.00   Min.   :0.00   Min.   :0.00  
#>  1st Qu.:1.00   1st Qu.:1.00   1st Qu.:1.00  
#>  Median :1.00   Median :1.00   Median :1.00  
#>  Mean   :0.94   Mean   :0.94   Mean   :0.89  
#>  3rd Qu.:1.00   3rd Qu.:1.00   3rd Qu.:1.00  
#>  Max.   :1.00   Max.   :1.00   Max.   :1.00  
#> 
summary(backdoor_md$simdata_obs)
#>       z_md             x_md             y_md             r_z      
#>  Min.   :0.0000   Min.   :0.0000   Min.   :0.0000   Min.   :0.00  
#>  1st Qu.:0.0000   1st Qu.:0.0000   1st Qu.:0.0000   1st Qu.:1.00  
#>  Median :0.0000   Median :0.0000   Median :0.0000   Median :1.00  
#>  Mean   :0.4043   Mean   :0.4149   Mean   :0.4157   Mean   :0.94  
#>  3rd Qu.:1.0000   3rd Qu.:1.0000   3rd Qu.:1.0000   3rd Qu.:1.00  
#>  Max.   :1.0000   Max.   :1.0000   Max.   :1.0000   Max.   :1.00  
#>  NAs    :6        NAs    :6        NAs    :11                     
#>       r_x            r_y      
#>  Min.   :0.00   Min.   :0.00  
#>  1st Qu.:1.00   1st Qu.:1.00  
#>  Median :1.00   Median :1.00  
#>  Mean   :0.94   Mean   :0.89  
#>  3rd Qu.:1.00   3rd Qu.:1.00  
#>  Max.   :1.00   Max.   :1.00  
#>

By using the argument fixedvars one can keep the complete data unchanged and re-simulate the missing data mechanism.

backdoor_md$simulate(100, fixedvars = c("x","y","z","ux","uy","uz"))
summary(backdoor_md$simdata)
#>        uz                 ux                 uy               urz         
#>  Min.   :0.002313   Min.   :0.003122   Min.   :0.01879   Min.   :0.01426  
#>  1st Qu.:0.285032   1st Qu.:0.176956   1st Qu.:0.26043   1st Qu.:0.29394  
#>  Median :0.483755   Median :0.474820   Median :0.47428   Median :0.48859  
#>  Mean   :0.498774   Mean   :0.470626   Mean   :0.47280   Mean   :0.51458  
#>  3rd Qu.:0.714295   3rd Qu.:0.738180   3rd Qu.:0.65695   3rd Qu.:0.78412  
#>  Max.   :0.998902   Max.   :0.989408   Max.   :0.99489   Max.   :0.99737  
#>                                                                           
#>       urx               ury                 z              x       
#>  Min.   :0.01785   Min.   :0.001106   Min.   :0.00   Min.   :0.00  
#>  1st Qu.:0.26595   1st Qu.:0.227070   1st Qu.:0.00   1st Qu.:0.00  
#>  Median :0.46274   Median :0.530725   Median :0.00   Median :0.00  
#>  Mean   :0.49163   Mean   :0.505790   Mean   :0.38   Mean   :0.43  
#>  3rd Qu.:0.71624   3rd Qu.:0.780873   3rd Qu.:1.00   3rd Qu.:1.00  
#>  Max.   :0.99958   Max.   :0.997806   Max.   :1.00   Max.   :1.00  
#>                                                                    
#>        y             z_md             x_md             y_md       
#>  Min.   :0.00   Min.   :0.0000   Min.   :0.0000   Min.   :0.0000  
#>  1st Qu.:0.00   1st Qu.:0.0000   1st Qu.:0.0000   1st Qu.:0.0000  
#>  Median :0.00   Median :0.0000   Median :0.0000   Median :0.0000  
#>  Mean   :0.47   Mean   :0.3678   Mean   :0.4066   Mean   :0.4362  
#>  3rd Qu.:1.00   3rd Qu.:1.0000   3rd Qu.:1.0000   3rd Qu.:1.0000  
#>  Max.   :1.00   Max.   :1.0000   Max.   :1.0000   Max.   :1.0000  
#>                 NAs    :13       NAs    :9        NAs    :6       
#>       r_z            r_x            r_y      
#>  Min.   :0.00   Min.   :0.00   Min.   :0.00  
#>  1st Qu.:1.00   1st Qu.:1.00   1st Qu.:1.00  
#>  Median :1.00   Median :1.00   Median :1.00  
#>  Mean   :0.87   Mean   :0.91   Mean   :0.94  
#>  3rd Qu.:1.00   3rd Qu.:1.00   3rd Qu.:1.00  
#>  Max.   :1.00   Max.   :1.00   Max.   :1.00  
#> 
summary(backdoor_md$simdata_obs)
#>       z_md             x_md             y_md             r_z      
#>  Min.   :0.0000   Min.   :0.0000   Min.   :0.0000   Min.   :0.00  
#>  1st Qu.:0.0000   1st Qu.:0.0000   1st Qu.:0.0000   1st Qu.:1.00  
#>  Median :0.0000   Median :0.0000   Median :0.0000   Median :1.00  
#>  Mean   :0.3678   Mean   :0.4066   Mean   :0.4362   Mean   :0.87  
#>  3rd Qu.:1.0000   3rd Qu.:1.0000   3rd Qu.:1.0000   3rd Qu.:1.00  
#>  Max.   :1.0000   Max.   :1.0000   Max.   :1.0000   Max.   :1.00  
#>  NAs    :13       NAs    :9        NAs    :6                      
#>       r_x            r_y      
#>  Min.   :0.00   Min.   :0.00  
#>  1st Qu.:1.00   1st Qu.:1.00  
#>  Median :1.00   Median :1.00  
#>  Mean   :0.91   Mean   :0.94  
#>  3rd Qu.:1.00   3rd Qu.:1.00  
#>  Max.   :1.00   Max.   :1.00  
#>

Applying Do-search to a missing data problem

backdoor_md$dosearch(data = "p(x*,y*,z*,r_x,r_y,r_z)", query = "p(y|do(x))")
#> \sum_{z}\left(\frac{p(z,r_z = 1)}{p(r_z = 1)}p(y|z,r_z = 1,x,r_x = 1,r_y = 1)\right)

It is automatically recognized that the problem is a missing data problem when rflist != NULL.

Using R6causal