Calculated distribution of a query from a prior or posterior distribution of parameters
Usage
query_distribution(
model,
queries,
given = NULL,
using = "parameters",
parameters = NULL,
n_draws = 4000,
join_by = "|",
case_level = FALSE,
query = NULL
)Arguments
- model
A
causal_model. A model object generated bymake_model.- queries
A character vector or list of character vectors specifying queries on potential outcomes such as "Y[X=1] - Y[X=0]"
- given
A character vector specifying givens for each query. A given is a quoted expression that evaluates to logical statement.
givenallows the query to be conditioned on *observational* distribution. A value of TRUE is interpreted as no conditioning.- using
A character. Whether to use priors, posteriors or parameters
- parameters
A vector or list of vectors of real numbers in [0,1]. A true parameter vector to be used instead of parameters attached to the model in case
usingspecifiesparameters- n_draws
An integer. Number of draws.rm
- join_by
A character. The logical operator joining expanded types when
querycontains wildcard (.). Can take values"&"(logical AND) or"|"(logical OR). When restriction contains wildcard (.) andjoin_byis not specified, it defaults to"|", otherwise it defaults toNULL.- case_level
Logical. If TRUE estimates the probability of the query for a case.
- query
alias for queries
Value
A DataFrame where columns contain draws from the distribution
of the potential outcomes specified in query
Examples
model <- make_model("X -> Y") %>%
set_parameters(c(.5, .5, .1, .2, .3, .4))
# \donttest{
# simple queries
query_distribution(model, query = "(Y[X=1] > Y[X=0])",
using = "priors") |>
head()
#> (Y[X=1] > Y[X=0])
#> 1 0.06111379
#> 2 0.44154416
#> 3 0.26766680
#> 4 0.36289837
#> 5 0.23667545
#> 6 0.04890897
# multiple queries
query_distribution(model,
query = list("(Y[X=1] > Y[X=0])",
"(Y[X=1] < Y[X=0])"),
using = "priors")|>
head()
#> (Y[X=1] > Y[X=0]) (Y[X=1] < Y[X=0])
#> 1 0.61847568 0.1268278
#> 2 0.16350915 0.1386330
#> 3 0.09975490 0.0578272
#> 4 0.57829491 0.0235060
#> 5 0.01991716 0.7515266
#> 6 0.17745259 0.5982866
# multiple queries and givens
query_distribution(model,
query = list("(Y[X=1] > Y[X=0])", "(Y[X=1] < Y[X=0])"),
given = list("Y==1", "(Y[X=1] <= Y[X=0])"),
using = "priors")|>
head()
#> (Y[X=1] > Y[X=0]) | Y==1 (Y[X=1] < Y[X=0]) | (Y[X=1] <= Y[X=0])
#> 1 0.13542224 0.57771672
#> 2 0.44302955 0.43379670
#> 3 0.32289961 0.06444665
#> 4 0.03733217 0.38763635
#> 5 0.11432872 0.51898372
#> 6 0.60695022 0.33701013
# linear queries
query_distribution(model, query = "(Y[X=1] - Y[X=0])")
#> (Y[X=1] - Y[X=0])
#> 1 0.1
# queries conditional on observables
query_distribution(model, query = "(Y[X=1] > Y[X=0])",
given = "X==1 & Y ==1")
#> (Y[X=1] > Y[X=0])
#> 1 0.4285714
# Linear query conditional on potential outcomes
query_distribution(model, query = "(Y[X=1] - Y[X=0])",
given = "Y[X=1]==0")
#> (Y[X=1] - Y[X=0])
#> 1 -0.6666667
# Use join_by to amend query interpretation
query_distribution(model, query = "(Y[X=.] == 1)", join_by = "&")
#> Generated expanded expression:
#> (Y[X=0] == 1 | Y[X=1] == 1)
#> (Y[X=.] == 1)
#> 1 0.9
# Probability of causation query
query_distribution(model,
query = "(Y[X=1] > Y[X=0])",
given = "X==1 & Y==1",
using = "priors") |> head()
#> (Y[X=1] > Y[X=0])
#> 1 0.1405846
#> 2 0.5870869
#> 3 0.9578517
#> 4 0.5280612
#> 5 0.8793493
#> 6 0.7937734
# Case level probability of causation query
query_distribution(model,
query = "(Y[X=1] > Y[X=0])",
given = "X==1 & Y==1",
case_level = TRUE,
using = "priors")
#> (Y[X=1] > Y[X=0])
#> 1 0.4950269
# Query posterior
update_model(model, make_data(model, n = 3)) |>
query_distribution(query = "(Y[X=1] - Y[X=0])", using = "posteriors") |>
head()
#>
#> SAMPLING FOR MODEL 'simplexes' NOW (CHAIN 1).
#> Chain 1:
#> Chain 1: Gradient evaluation took 0.000201 seconds
#> Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 2.01 seconds.
#> Chain 1: Adjust your expectations accordingly!
#> Chain 1:
#> Chain 1:
#> Chain 1: Iteration: 1 / 2000 [ 0%] (Warmup)
#> Chain 1: Iteration: 200 / 2000 [ 10%] (Warmup)
#> Chain 1: Iteration: 400 / 2000 [ 20%] (Warmup)
#> Chain 1: Iteration: 600 / 2000 [ 30%] (Warmup)
#> Chain 1: Iteration: 800 / 2000 [ 40%] (Warmup)
#> Chain 1: Iteration: 1000 / 2000 [ 50%] (Warmup)
#> Chain 1: Iteration: 1001 / 2000 [ 50%] (Sampling)
#> Chain 1: Iteration: 1200 / 2000 [ 60%] (Sampling)
#> Chain 1: Iteration: 1400 / 2000 [ 70%] (Sampling)
#> Chain 1: Iteration: 1600 / 2000 [ 80%] (Sampling)
#> Chain 1: Iteration: 1800 / 2000 [ 90%] (Sampling)
#> Chain 1: Iteration: 2000 / 2000 [100%] (Sampling)
#> Chain 1:
#> Chain 1: Elapsed Time: 1.558 seconds (Warm-up)
#> Chain 1: 1.477 seconds (Sampling)
#> Chain 1: 3.035 seconds (Total)
#> Chain 1:
#>
#> SAMPLING FOR MODEL 'simplexes' NOW (CHAIN 2).
#> Chain 2:
#> Chain 2: Gradient evaluation took 0.000128 seconds
#> Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 1.28 seconds.
#> Chain 2: Adjust your expectations accordingly!
#> Chain 2:
#> Chain 2:
#> Chain 2: Iteration: 1 / 2000 [ 0%] (Warmup)
#> Chain 2: Iteration: 200 / 2000 [ 10%] (Warmup)
#> Chain 2: Iteration: 400 / 2000 [ 20%] (Warmup)
#> Chain 2: Iteration: 600 / 2000 [ 30%] (Warmup)
#> Chain 2: Iteration: 800 / 2000 [ 40%] (Warmup)
#> Chain 2: Iteration: 1000 / 2000 [ 50%] (Warmup)
#> Chain 2: Iteration: 1001 / 2000 [ 50%] (Sampling)
#> Chain 2: Iteration: 1200 / 2000 [ 60%] (Sampling)
#> Chain 2: Iteration: 1400 / 2000 [ 70%] (Sampling)
#> Chain 2: Iteration: 1600 / 2000 [ 80%] (Sampling)
#> Chain 2: Iteration: 1800 / 2000 [ 90%] (Sampling)
#> Chain 2: Iteration: 2000 / 2000 [100%] (Sampling)
#> Chain 2:
#> Chain 2: Elapsed Time: 1.416 seconds (Warm-up)
#> Chain 2: 1.243 seconds (Sampling)
#> Chain 2: 2.659 seconds (Total)
#> Chain 2:
#>
#> SAMPLING FOR MODEL 'simplexes' NOW (CHAIN 3).
#> Chain 3:
#> Chain 3: Gradient evaluation took 0.00012 seconds
#> Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 1.2 seconds.
#> Chain 3: Adjust your expectations accordingly!
#> Chain 3:
#> Chain 3:
#> Chain 3: Iteration: 1 / 2000 [ 0%] (Warmup)
#> Chain 3: Iteration: 200 / 2000 [ 10%] (Warmup)
#> Chain 3: Iteration: 400 / 2000 [ 20%] (Warmup)
#> Chain 3: Iteration: 600 / 2000 [ 30%] (Warmup)
#> Chain 3: Iteration: 800 / 2000 [ 40%] (Warmup)
#> Chain 3: Iteration: 1000 / 2000 [ 50%] (Warmup)
#> Chain 3: Iteration: 1001 / 2000 [ 50%] (Sampling)
#> Chain 3: Iteration: 1200 / 2000 [ 60%] (Sampling)
#> Chain 3: Iteration: 1400 / 2000 [ 70%] (Sampling)
#> Chain 3: Iteration: 1600 / 2000 [ 80%] (Sampling)
#> Chain 3: Iteration: 1800 / 2000 [ 90%] (Sampling)
#> Chain 3: Iteration: 2000 / 2000 [100%] (Sampling)
#> Chain 3:
#> Chain 3: Elapsed Time: 1.356 seconds (Warm-up)
#> Chain 3: 1.502 seconds (Sampling)
#> Chain 3: 2.858 seconds (Total)
#> Chain 3:
#>
#> SAMPLING FOR MODEL 'simplexes' NOW (CHAIN 4).
#> Chain 4:
#> Chain 4: Gradient evaluation took 0.000127 seconds
#> Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 1.27 seconds.
#> Chain 4: Adjust your expectations accordingly!
#> Chain 4:
#> Chain 4:
#> Chain 4: Iteration: 1 / 2000 [ 0%] (Warmup)
#> Chain 4: Iteration: 200 / 2000 [ 10%] (Warmup)
#> Chain 4: Iteration: 400 / 2000 [ 20%] (Warmup)
#> Chain 4: Iteration: 600 / 2000 [ 30%] (Warmup)
#> Chain 4: Iteration: 800 / 2000 [ 40%] (Warmup)
#> Chain 4: Iteration: 1000 / 2000 [ 50%] (Warmup)
#> Chain 4: Iteration: 1001 / 2000 [ 50%] (Sampling)
#> Chain 4: Iteration: 1200 / 2000 [ 60%] (Sampling)
#> Chain 4: Iteration: 1400 / 2000 [ 70%] (Sampling)
#> Chain 4: Iteration: 1600 / 2000 [ 80%] (Sampling)
#> Chain 4: Iteration: 1800 / 2000 [ 90%] (Sampling)
#> Chain 4: Iteration: 2000 / 2000 [100%] (Sampling)
#> Chain 4:
#> Chain 4: Elapsed Time: 1.404 seconds (Warm-up)
#> Chain 4: 1.229 seconds (Sampling)
#> Chain 4: 2.633 seconds (Total)
#> Chain 4:
#> (Y[X=1] - Y[X=0])
#> 1 -0.5747837
#> 2 -0.3171723
#> 3 -0.1635851
#> 4 -0.7593423
#> 5 -0.2373443
#> 6 0.2677196
# Case level queries provide the inference for a case, which is a scalar
# The case level query *updates* on the given information
# For instance, here we have a model for which we are quite sure that X
# causes Y but we do not know whether it works through two positive effects
# or two negative effects. Thus we do not know if M=0 would suggest an
# effect or no effect
set.seed(1)
model <-
make_model("X -> M -> Y") |>
update_model(data.frame(X = rep(0:1, 8), Y = rep(0:1, 8)), iter = 10000)
#>
#> SAMPLING FOR MODEL 'simplexes' NOW (CHAIN 1).
#> Chain 1:
#> Chain 1: Gradient evaluation took 0.000245 seconds
#> Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 2.45 seconds.
#> Chain 1: Adjust your expectations accordingly!
#> Chain 1:
#> Chain 1:
#> Chain 1: Iteration: 1 / 10000 [ 0%] (Warmup)
#> Chain 1: Iteration: 1000 / 10000 [ 10%] (Warmup)
#> Chain 1: Iteration: 2000 / 10000 [ 20%] (Warmup)
#> Chain 1: Iteration: 3000 / 10000 [ 30%] (Warmup)
#> Chain 1: Iteration: 4000 / 10000 [ 40%] (Warmup)
#> Chain 1: Iteration: 5000 / 10000 [ 50%] (Warmup)
#> Chain 1: Iteration: 5001 / 10000 [ 50%] (Sampling)
#> Chain 1: Iteration: 6000 / 10000 [ 60%] (Sampling)
#> Chain 1: Iteration: 7000 / 10000 [ 70%] (Sampling)
#> Chain 1: Iteration: 8000 / 10000 [ 80%] (Sampling)
#> Chain 1: Iteration: 9000 / 10000 [ 90%] (Sampling)
#> Chain 1: Iteration: 10000 / 10000 [100%] (Sampling)
#> Chain 1:
#> Chain 1: Elapsed Time: 16.707 seconds (Warm-up)
#> Chain 1: 17.177 seconds (Sampling)
#> Chain 1: 33.884 seconds (Total)
#> Chain 1:
#>
#> SAMPLING FOR MODEL 'simplexes' NOW (CHAIN 2).
#> Chain 2:
#> Chain 2: Gradient evaluation took 0.000226 seconds
#> Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 2.26 seconds.
#> Chain 2: Adjust your expectations accordingly!
#> Chain 2:
#> Chain 2:
#> Chain 2: Iteration: 1 / 10000 [ 0%] (Warmup)
#> Chain 2: Iteration: 1000 / 10000 [ 10%] (Warmup)
#> Chain 2: Iteration: 2000 / 10000 [ 20%] (Warmup)
#> Chain 2: Iteration: 3000 / 10000 [ 30%] (Warmup)
#> Chain 2: Iteration: 4000 / 10000 [ 40%] (Warmup)
#> Chain 2: Iteration: 5000 / 10000 [ 50%] (Warmup)
#> Chain 2: Iteration: 5001 / 10000 [ 50%] (Sampling)
#> Chain 2: Iteration: 6000 / 10000 [ 60%] (Sampling)
#> Chain 2: Iteration: 7000 / 10000 [ 70%] (Sampling)
#> Chain 2: Iteration: 8000 / 10000 [ 80%] (Sampling)
#> Chain 2: Iteration: 9000 / 10000 [ 90%] (Sampling)
#> Chain 2: Iteration: 10000 / 10000 [100%] (Sampling)
#> Chain 2:
#> Chain 2: Elapsed Time: 16.866 seconds (Warm-up)
#> Chain 2: 17.604 seconds (Sampling)
#> Chain 2: 34.47 seconds (Total)
#> Chain 2:
#>
#> SAMPLING FOR MODEL 'simplexes' NOW (CHAIN 3).
#> Chain 3:
#> Chain 3: Gradient evaluation took 0.000186 seconds
#> Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 1.86 seconds.
#> Chain 3: Adjust your expectations accordingly!
#> Chain 3:
#> Chain 3:
#> Chain 3: Iteration: 1 / 10000 [ 0%] (Warmup)
#> Chain 3: Iteration: 1000 / 10000 [ 10%] (Warmup)
#> Chain 3: Iteration: 2000 / 10000 [ 20%] (Warmup)
#> Chain 3: Iteration: 3000 / 10000 [ 30%] (Warmup)
#> Chain 3: Iteration: 4000 / 10000 [ 40%] (Warmup)
#> Chain 3: Iteration: 5000 / 10000 [ 50%] (Warmup)
#> Chain 3: Iteration: 5001 / 10000 [ 50%] (Sampling)
#> Chain 3: Iteration: 6000 / 10000 [ 60%] (Sampling)
#> Chain 3: Iteration: 7000 / 10000 [ 70%] (Sampling)
#> Chain 3: Iteration: 8000 / 10000 [ 80%] (Sampling)
#> Chain 3: Iteration: 9000 / 10000 [ 90%] (Sampling)
#> Chain 3: Iteration: 10000 / 10000 [100%] (Sampling)
#> Chain 3:
#> Chain 3: Elapsed Time: 17.212 seconds (Warm-up)
#> Chain 3: 19.106 seconds (Sampling)
#> Chain 3: 36.318 seconds (Total)
#> Chain 3:
#>
#> SAMPLING FOR MODEL 'simplexes' NOW (CHAIN 4).
#> Chain 4:
#> Chain 4: Gradient evaluation took 0.00018 seconds
#> Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 1.8 seconds.
#> Chain 4: Adjust your expectations accordingly!
#> Chain 4:
#> Chain 4:
#> Chain 4: Iteration: 1 / 10000 [ 0%] (Warmup)
#> Chain 4: Iteration: 1000 / 10000 [ 10%] (Warmup)
#> Chain 4: Iteration: 2000 / 10000 [ 20%] (Warmup)
#> Chain 4: Iteration: 3000 / 10000 [ 30%] (Warmup)
#> Chain 4: Iteration: 4000 / 10000 [ 40%] (Warmup)
#> Chain 4: Iteration: 5000 / 10000 [ 50%] (Warmup)
#> Chain 4: Iteration: 5001 / 10000 [ 50%] (Sampling)
#> Chain 4: Iteration: 6000 / 10000 [ 60%] (Sampling)
#> Chain 4: Iteration: 7000 / 10000 [ 70%] (Sampling)
#> Chain 4: Iteration: 8000 / 10000 [ 80%] (Sampling)
#> Chain 4: Iteration: 9000 / 10000 [ 90%] (Sampling)
#> Chain 4: Iteration: 10000 / 10000 [100%] (Sampling)
#> Chain 4:
#> Chain 4: Elapsed Time: 17.337 seconds (Warm-up)
#> Chain 4: 15.37 seconds (Sampling)
#> Chain 4: 32.707 seconds (Total)
#> Chain 4:
Q <- "Y[X=1] > Y[X=0]"
G <- "X==1 & Y==1 & M==1"
QG <- "(Y[X=1] > Y[X=0]) & (X==1 & Y==1 & M==1)"
# In this case these are very different:
query_distribution(model, Q, given = G, using = "posteriors")[[1]] |> mean()
#> [1] 0.4320531
query_distribution(model, Q, given = G, using = "posteriors",
case_level = TRUE)
#> Y[X=1] > Y[X=0]
#> 1 0.6760235
# These are equivalent:
# 1. Case level query via function
query_distribution(model, Q, given = G,
using = "posteriors", case_level = TRUE)
#> Y[X=1] > Y[X=0]
#> 1 0.6760235
# 2. Case level query by hand using Bayes
distribution <- query_distribution(
model, list(QG = QG, G = G), using = "posteriors")
mean(distribution$QG)/mean(distribution$G)
#> [1] 0.6760235
# }