causact: From narrative to computational insight

An Initial Step

Let’s get started capturing the BankPass problem as a graphical model. The first thought should be to make a simple model of how the real-world ends up generating your variable(s) of interest; and don’t worry about how good the model is - just get something written.

library(causact)
dag_create() %>%
  dag_node("Get Card", "Y") %>%
  dag_node("Adventurous Personality", "X",
           child = "Y") %>%
  dag_render(shortLabel = TRUE)

The above output is a visual depiction of what mathematicians and computer scientists call a graph. To them, a graph is a set of related objects where the objects are called nodes and any pair of related nodes are known as an edge. For our purposes, a node is a visual depiction of a random variable (i.e. an oval) and the presence of a visible edge indicates a relationship between the connected random variables.

Ensuring Nodes Are RV’s

Vetting the above model reveals some shortcomings and our first vetting criteria is to ensure the two nodes can be modelled as random variables with each node representing a mapping of real-world outcomes to a numerical values.

For the “Get Card” node, let’s label it X, the mapping of outcomes to values is the conventional mapping for a Bernoulli-distributed random variable:

$$ X \equiv \begin{cases} 0, & \textrm{if customer does not get the card }\\ 1, & \textrm{if customer gets the TRC card } \end{cases} $$

where successes are mapped to 1 and failures mapped to 0.

However, how does one map “Adventurous Personality” to values? In fact, what measurement or specific criteria am I using for “Adventurous Personality”? This is an ill-defined variable. Since the story actually suggests that Car Model is the data we have, let’s create a graphical model where Car Model helps determine the likelihood of Get Card. Using Y and X as respective labels for Car Model and Get Card, we show a revised candidate graphical model below.

library(causact)
dag_create() %>%
  dag_node("Get Card", "X") %>%
  dag_node("Car Model", "Y",
           child = "X") %>%
  dag_render(shortLabel = FALSE)

As shown above, the difficult to measure personality node gets replaced by the car model node. To verify we can measure Car Model as a random variable, let’s be explicit about the mapping of “Car Model” outcomes to real numbers. For example:

$$ Y \equiv \begin{cases} 1, & \textrm{if customer owns a Toyota Corolla}\\ 2, & \textrm{if customer owns a Jeep Wrangler}\\ \textrm{ } \vdots & \hspace{3cm} \vdots\\ K, & \textrm{if customer owns the } K^{th} \textrm{car model.} \end{cases} $$

Adding data, likelihood, and prior

A dataset to accompany this story is built-in to the causact package:

# ------------------------------------
# Revisiting the Get Credit Card Model
# ------------------------------------

# data included in the causact package
# 1,000 customers and the car they drive
carModelDF   #show data
#> # A tibble: 1,000 x 3
#> # Rowwise: 
#>    customerID carModel       getCard
#>         <int> <chr>            <int>
#>  1          1 Toyota Corolla       0
#>  2          2 Subaru Outback       1
#>  3          3 Toyota Corolla       0
#>  4          4 Subaru Outback       0
#>  5          5 Kia Forte            1
#>  6          6 Toyota Corolla       0
#>  7          7 Toyota Corolla       0
#>  8          8 Subaru Outback       1
#>  9          9 Toyota Corolla       0
#> 10         10 Toyota Corolla       0
#> # ... with 990 more rows

We take a bottom-up approach starting with our target measurement Get Card to add in mathematical details about how the data is generated. Specifically, we declare any observation x to be a Bernoulli-distributed random variable with yet-to-be-defined parameter theta. (For a list of supported distributions type ?causact::distributions)

To model additional nodes, consider any parameters on the right-hand side of the statistical model that remain undefined; in our case theta is undefined. Hence, if x ∼ Bernoulli(θ), then θ should be a parent node of x in the graphical model. And at this point, we say ay caramba!. There is a mismatch between Car Model, the parent of Get Card in the previous graphical model, and θ, the implied parent of Get Card using the statistical model.

Realizing that Signup Probability, θ, is truly the variable of interest as we are interested in modelling the probability of someone signing up for the card. The below figure graphically models the probability of sign-up (θ) as the parent. Notice that nodes with data, i.e. Get Card, have darker fill; observed nodes get a darker fill and latent or unobserved nodes remain lighter. The right-hand side (rhs) of darker nodes represent the likelihood part of Bayes rule and the right-hand side of lighter nodes reprent the prior part of Bayes rule.

graph = dag_create() %>%
  dag_node("Get Card","x",
           rhs = bernoulli(theta),
           data = carModelDF$getCard) %>%
  dag_node("Signup Probability","theta",
           rhs = uniform(0,1),
           child = "x")
graph %>% dag_render()

The above generative recipe implies there is only one true θ value - which is unknown so we place a flat uniform prior over all possible values for θ. Given θ, 1,000 observations are outcomes of Bernoulli trials using that one θ value. The posterior distribution reallocates plausibility among θ values to those most consistent with our prior and data. A representative sample from the posterior distribution is easily obtained using dag_numpyro():

drawsDF = graph %>% dag_numpyro()
drawsDF %>% dagp_plot()

From this, a reasonable conclusion might be that plausible values for Signup Probability lie (roughly) between 36% and 42%.

Modelling the effect of car model with plates

Our business narrative calls for the probability to change based on car model. The dag_plate() function helps create a generative DAG where theta will vary based on car model. In comparison to mode probabilistic programming paradigms, the causact package shifts the burden of index creation to the computer instead of the analyst. A plate signifies that the random variables inside of it should be repeated; thus, any oval inside a plate is indexed behind the scenes and effectively, a repeat instance is created for each member of the plate. The system automatically adds the index to the model using the plate label for indexing notation. When using dag_plate() the following arguments prove most useful:

• descr - a character string representing a longer more descriptive label for the cluster/plate.
• label - a short character string to use as an index label when generating code.
• data - a vector representing the categorical data whose unique values become the plate index. When using with addDataNode = TRUE, ensure this vector represents observations of a variable whose unique values can be coerced to a factor. Omit this argument and causact assumes that node replication is based on observations of any observed nodes within the plate.
• nodeLabels - a character vector of node labels or descriptions to include in the list of nodes that get repeated one time per index value. Labels must be defined already within the specification of the causact_graph_.
• addDataNode - a logical value. When addDataNode = TRUE, the code adds a node of observed data, in addition to the plate, based on the supplied data argument value. This node is used to extract the proper index when any of the nodes in the plate are fed to children nodes. Verify the graphical model using dag_render() to ensure correct behavior.

Below, we can change the generative DAG to accommodate the story of Signup Probability being different for drivers of different cars.

### Now vary your modelling of theta by car model
### we use the plate notation to say
### "create a copy of theta for each unique car model
graph = dag_create() %>%
  dag_node("Get Card","x",
           rhs = bernoulli(theta),
           data = carModelDF$getCard) %>%
  dag_node(descr = "Signup Probability",
           label = "theta",
           rhs = uniform(0,1),
           child = "x") %>%
  dag_plate(descr = "Car Model",
            label = "y",
            data = carModelDF$carModel, #plate labels
            nodeLabels = "theta",  # node on plate
            addDataNode = TRUE)  # create obs. node
graph %>% dag_render()

This generative DAG used one additional function call, dag_plate(), that resulted in both the creation of the plate, and also the additional node representing Car Model observations. Each of the 1,000 observed customers of the Car Model node were driving one of 4 unique car models of the Car Model plate. The indexing of theta by y, i.e. the theta[y] on the rhs of the Get Card node extracts one of the 4 θ-values to govern each observation of x.

Expert tip: Observed-discrete random variables, like car model, are best modelled via plates. Do not rush to add nodes or indexes for these random variables. Let dag_plate() do the work for you.

Now, when we get the posterior distribution from our code, we get a representative sample of four θ values instead of one.

## get representative sample
drawsDF = graph %>% dag_numpyro()

Visualizing the posterior yields:

drawsDF %>% dagp_plot()

Notice, the index names for the four theta values are automatically generated as abbreviated versions of the unique labels found in carModelDF$carModel.

Digesting the posterior plot, notice that credible interval widths are different for the different theta values. For example, the probability interval for the Toyota Corolla is more narrow than the interval for the other cars. This is simply because there is more data regarding Toyota Corolla drivers than the other vehicles. With more data, we get more confident in our estimate. Here is a quick way to see the distribution of car models across our observations:

library(tidyverse)
carModelDF %>%
  group_by(carModel) %>%
  summarize(numOfObservations = n())
#> # A tibble: 4 x 2
#>   carModel       numOfObservations
#>   <chr>                      <int>
#> 1 Jeep Wrangler                193
#> 2 Kia Forte                     86
#> 3 Subaru Outback               145
#> 4 Toyota Corolla               576

Making A Probabilistic Statement With An Indicator Function

In the previous posterior plot, notice the overlapping beliefs regarding the true sign-up probability for drivers of the Toyota Corolla and Kia Forte. Based on this, a reasonable query might be what is the probability that theta_Kia Forte > theta_ToyotaCorolla? As a typical workflow, use an indicator function on the posterior draws to answer the question. So, to answer “what is the probability that theta_Kia.Forte > theta_Toyota.Corolla?”, we compute:

drawsDF %>%
  mutate(indicatorFunction =
           theta_Kia.Forte > theta_Toyota.Corolla) %>%
  summarize(pctOfDraws = mean(indicatorFunction))
#> # A tibble: 1 x 1
#>   pctOfDraws
#>        <dbl>
#> 1      0.769

And using this ability to make probabalistic statements, one might say “there is an approximately 77% chance that Kia Forte drivers are more likely than Toyota Corolla drivers to sign up for the credit card.