Bayes' Town: A place of data simulation

Generate Bayes’ Town

With this knowledge we can generate data on our citizens for this particular scenario. The core functionality is powered by rJAGS but to make it a bit easier to use, the process is wrapped up into an R6 object called BayesTown. This should be familiar to anyone from object-orientated programming (OOP) languages like Python but perhaps a bit foreign for those more comfortable with multiple dispatch, which is the common approach in R (S3 and S4) and Julia. To create a population of 1,000 people we provide some predefined information to the new (construct) method. The inputParam and simModel values will be discussed later in the methodology section.

newTown <- BayesTown$new(population_size = 1000,
                      jags_data = inputParam, 
                      n.chains = 2,
                      n.iter = 5000,
                      jags_model = simModel)

The output provides some basic information on the model graph and a snapshot of the data sample. This includes all the information from the simulation, including the various data variables, coefficients, and distribution parameters. Due to sampling variation, subsequent runs with the same parameters will create slightly different populations.

Explore simulated population

To explore the simulated data, we extract the variables we observed in the diagram above. This is easily done by method chaining set_variables (to sample only particular variables/parameters) and resample (to rerun the simulation).

# Import library for some exploration work
library(dplyr)
library(magrittr)

# Data sample just to include the set of variables in the list
newTown$set_variables(c('age', 'sex', 'urbanRural', 'disease'), mode = 'include')$resample()

# Output first 5 rows of simulated data
head(newTown$data, 5)

# Summarize some basic information using dplyr
newTown$data %>% 
  group_by(sex, urbanRural) %>%
  summarise(meanAge = mean(age),
            diseasePerc = 100 * (sum(disease)/nrow(.))) 

Here is a sample of the simulated data:

##         age disease sex urbanRural
## 1  24.43352       1   0          1
## 2  22.56561       0   1          1
## 3  28.10493       0   1          1
## 4  33.43865       0   0          1
## 5  37.12840       0   1          1
## 6  56.35582       1   1          1
## 7  37.24761       0   1          1
## 8  47.65436       0   0          0
## 9  30.83504       0   1          1
## 10 33.03936       0   0          1

And here are some summary statistics:

## # A tibble: 4 x 4
## # Groups:   sex [2]
##     sex urbanRural meanAge diseasePerc
##   <dbl>      <dbl>   <dbl>       <dbl>
## 1     0          0    36.4         1.6
## 2     0          1    35.1        11.5
## 3     1          0    40.4         1.7
## 4     1          1    40.3        14.9

Just from this basic summary, we can see that (a) urban areas have a larger burden of disease and (b) females and older age groups on average have a slightly higher proportion of disease. Tables are nice, but pictures are better. Let’s create a few plots to see how the characteristics of Bayes’ Town citizens relate to disease status.

Age among females is slightly higher on average than their male counterparts.

We could plot each variable against the outcome (disease) separately but since the data has just 4 variables it is possible to observe the relationships at the same time. This may not be always preferable, as the visualization can become quite dense; that is to say, it could be hard to digest all of the dimensions presented from the data. How the variables are assigned to plot are important for aiding interpretation. As our interest is the presence/absence of Disease X, we assign this variable as the plot color to ensure it is included in each panel. In this instance plotting all the variables does provide some interesting insights:

• Males appear to have a larger representation in the younger age groups
  
• Disease is present in higher proportions among older age groups

Although I used a density plot in the example above, some may prefer using a box-plot as a suitable alternative.

Age seems to have a strong impact on the disease outcome which may be masking more subtle relationships. When we remove age we observe urban areas have a higher proportion of Disease X, and there is subtle evidence of a higher proportion of disease among females. None of this should come as a surprise, as we set these relationships in advance but it is nice to validate them through these figures.

Simulation toolbox

Agreeing on the usefulness of simulating data, how can it be done? There are many approaches with varying complexity:

• Agent/individual-based models
• Random variables simulated from probability distributions
• Creating Bayesian networks based upon priors and/or learned from data

Most approaches are available in R or its host of packages, albeit with a tradeoff between flexibility and complexity. For instance, NetLogo is a programming environment for agent-based simulation and can connect to R through the RNetLogo package; this however still requires knowledge of the NetLogo syntax. Other simulation approaches in R are available from HydeNet, greta, RStan, and simstudy. Some packages are limited as they only allow certain data types (e.g. discrete or continuous) and the possible relationships that can exist between them (e.g. linear trend). Typically, ‘low-level’ interfaces (i.e. a package that gets you close to the ‘engine’) allow more freedom but increase the learning curve and chances of blunder. Furthermore, a selection of methods listed require data to operate and for some purposes this may be unsuitable.

I am partial to Bayesian networks, which is my method of choice in this post, as I find their associated diagrams a powerful visualization tool to understand the model. They can also be created from our imagination, or trained on any available data we happen to have. Flexible and powerful!

Basic simulation example

Before I continue to explain the method I used to create Bayes’ Town, now is a good time to demonstrate that simulation techniques are not out of anyone’s reach.

Keeping with the theme of health, let’s simulate the probability of death by boredom in Bayes’ Town. Hopefully you are not suffering from this affliction as it is quite serious, as we will soon discover.

As a binary outcome, we can use the binomial distribution for modeling death. The probability of death will be a function of an intercept term and the binary explanatory variable for ‘boredom’ (0: not bored, 1: bored). The logit link will map the linear combination of the intercept and ‘boredom’ term between 0% and 100% probability. In math jargon this looks like:

\[\begin{align*} & y_i \sim Binomial(n, p_i) \\ & logit(p_i) = \alpha + \beta x_i\\ \end{align*}\]

If you are like me, you may get nervous when a page becomes overtaken by formulae. I promise, it is not as bad as it looks. I usually feel more at ease once I see the formulae translated to code, which we will do now.

# Set 100,000 people in our population
n <- 1e5

# Half of our population is bored
bored <- rep(c(0,1), each = n/2)

# Bored has a small increase to the log-odds of death
beta <- 0.3

# The baseline chance of death is low
alpha <- 0.05

# The log-odds of the outcome (arithmetic from the formula above!)
logodds <- alpha + beta * bored

# Probability of death (size = 1, only one death possible per person!)
death_probability <- rbinom(n, prob = plogis(logodds), size = 1)

# Final dataset
data <- data.frame(boredom = bored,
                   death = death_probability)

We can see a relationship between boredom and death in a simple 2x2 table. The values on the bottom of the table compare the odds of not being bored among the living and dead; these results show a 1.4x higher odds of death among those experiencing boredom.

library(knitr)
library(kableExtra)

# Function to calculate odds
odds <- function(x) round(x[[1]]/x[[2]], 3)

# Create table with odds, adjust rounding, and format
summary_tbl <- table(Boredom = data$boredom,
                     Death = data$death) %>% 
  addmargins(margin = c(1,2),
             FUN = list(Odds = odds),
             quiet = T)

summary_tbl[c(1:3), c(1,3)] <- as.character(summary_tbl[c(1:3), c(1,3)])

summary_tbl %>%
  kable(caption = 'Relationship between boredom and death.') %>%
  kable_styling(full_width = T) %>%
  pack_rows('Boredom', 1,2, label_row_css = 'border-bottom:0px solid;') %>%
  add_header_above(c(" ", "Death" = 3), align = 'l') %>%
  scroll_box(box_css = "", extra_css = 'overflow-x: auto; width: 100%')

If that isn’t enough to convince you that the simulated data created a positive correlation between boredom and death, let’s pass the data through a logisitic regression model to see how well we can retrieve the original parameters for the intercept (0.05) and beta term (0.3).

coef(glm(death ~ boredom, data = data, family = 'binomial')) %>%
  round(3) %>%
  kable(col.names = 'Estimate', caption = 'Parameter estimates for model on boredom and death.') %>%
  kable_styling() %>%
  scroll_box(box_css = "",
             extra_css = 'overflow-x: auto; width: 50%')

That is pretty close! We could take this simulation further and create a dependency of boredom upon additional terms. I encourage you to explore other situations, such as collinear variables. As the simulation becomes increasingly complex other approaches may be easier to maintain. This is what we will discuss next.

Model definition

Similar to our simple example, we will start by assigning our parameter priors, essentially a set of distributions for the intercept and coefficient terms in the model.

\[ \begin{align*} & \textbf{Prior parameters}\\ & a_{disease} \sim Normal(-10, 10) \\ & b_{age} \sim Normal(0.015, 10) \\ & b_{sex} \sim Normal(0.005, 10) \\ & b_{urbanity} \sim Normal(5, 10) \\ \\ \end{align*} \] \(a_{disease}\) is the baseline odds of disease and the other parameters define the magnitude by which they increase the odds of disease. Next, we create distributions for our variables of interest (alternatively, these could be set as fixed values).

\[ \begin{align*} & \textbf{Random variable priors}\\ & sex \sim Binom(propMale, 1) \\ & urbanRural \sim Binom(0.8, 1) \\ & \begin{aligned} age \sim Gamma((age_{shape} + &age\_param * sex),\\ 1/age_{scale}) \end{aligned}\\ \\ \end{align*} \]

We also will have several input parameters for our variables. We can adjust these or the prior parameters to create strong/weaker, or entirely different, relationships.

\[ \begin{align*} & \textbf{Input parameters}\\ & propMale = 0.52 \\ & age_{shape} = 6 \\ & age_{scale} = 6 \\ & age\_param = 0.7 \\ \\ \end{align*} \]

And finally, we need to assign a likelihood function for the prior predictive simulation. This is the distribution for Disease X based upon the relationship of age, sex, and location.

\[ \begin{align*} & \textbf{Likelihood function}\\ & disease \sim Binom(p, 1) \\ & \begin{aligned} logit(p) &= a_{disease} \\ &+ b_{age}(age) \\ &+ b_{sex}(sex) \\ &+ b_{urbanity}(urbanRural)\\\\ \end{aligned} \end{align*} \]

These formulae may look confusing, but they are outlining very similar concepts to the simple example between boredom and death. If you are a little lost, the DAG from the beginning should help. In plain English we are saying that:

• Sex is slightly more likely to be male
• Location has a much higher chance to be urban
• Age is slightly dependent upon sex; females reach slightly older age groups on average
• Disease X is a linear relationship between age, sex, and location
• The strength of age, sex, and location are pulled from a distribution of values that we can tweak to increase/decrease their strength.

In R code it looks quite similar:

# Set input parameters (could be included directly in the model below if we arent changing often)
inputParam <- list(p_sex = 0.52, 
                   age_shape = 6,
                   age_scale = 6,
                   age_param = 0.7)

# Define model (character string is code for JAGS)
simModel <-  as.character('
                      model {
                        a_disease ~ dnorm(-10, 10);
                        b_age ~ dnorm(.015, 10);
                        b_sex ~ dnorm(.005, 10);
                        b_urbanity ~ dnorm(5, 10);
                      
                        sex ~ dbinom(p_sex, 1);
                        urbanRural ~ dbinom(0.8, 1);
                        age ~ dgamma((age_shape + age_param * sex), 1/age_scale);
                      
                        logit(p) <- a_disease + b_age * age + b_sex * sex + b_urbanity * urbanRural;
                        disease ~ dbinom(p, 1);
                      }')

This is the set of rules that governs Bayes’ Town citizens. Right now it provides a single snapshot but it is possible to extend this to a dynamic BN to create multiple time slices.

Bayes’ Theorem is used to solve this type of question and the equation may look simple enough… \[P(A|B) = \frac{P(B|A)P(A)}{P(B)}\] I do not provide details here as there are a plethora of resources online that go through Bayes’ Theorem; however, I will point out that a special ‘engine’ is often required to solve these models. When the distributions are complex, the solutions become non-trivial. Whilst some examples such as the beta-binomial have closed-form solutions, others require special sampling methods to solve, which we will discuss next.

MCMC and rjags

With the rules defined to create Disease X in Bayes’ Town, we now need to perform the simulation. The procedure used is called markov chain monte carlo (MCMC), a stochastic method to estimate posterior probability distributions. In short this means we collect samples of possible parameter values instead of attempting to directly compute them. The results are then understood in terms of their sample distribution. However, generating these samples is not always easy and there are various MCMC methods to tackle this challenge, albeit with their own limitations and trade-offs.

1. Metropolis algorithm: original and ‘simple’ MCMC implementation.
2. Metropolis-Hastings algorithm: an extension to the metropolis algorithm.
3. Gibbs sampling: more efficient sampling by using adaptive proposals (BUGS and JAGS software use this method).
4. Hamiltonian Monte Carlo (HMC): runs a (resource expensive) physics simulation to generate informative samples. HMC provides useful diagnostics but currently cannot sample discrete parameters (not easily at least). STAN uses this method and is considered harder to learn than BUGS.

Numerous R packages implement the assortment of MCMC methods but require, or are built on top of, an R interface to STAN (e.g. rstan) or JAGS (e.g. rjags). For my current purposes, I decided to use rjags (authored by Martyn Plummer who is also the creator of JAGS) since it can generate discrete values and follows a syntax similar to R. Using rjags directly is a bit more complex (as we saw in the model definitions) but provides flexibility and customization.

There are three steps to run a model with JAGS: (1) Define, (2) Compile, and (3) Simulate. We have already done the first step, for the rest I created helper functions to pass the model to JAGS and retrieve the samples! This may be unnecessary but I found it to be a helpful abstraction. We can always dip into rjags if we need something more specific.

The R6 package was used to encapsulate the BayesTown object created from the simulation. This object is created by running BayesTown$new() with the provided model (simModel) and input parameters (inputParam). JAGS will then compile and perform the simulation, returning the results to the assigned object class. A subset of variables can be retained via the set_variables() method and a new draw can be done via resample(). Various diagnostics are also accessible. If this is still a bit unclear, the complete example in the next section may be of help.

Useful data sources

1. Kaggle
2. WHO
3. US Gov
4. CDC
5. NHS data catalogue
6. Alberta Open Datasets
7. Canada Open Data
   
8. UCI Machine Learning Repository

Software and details on simulation methods

1. NetLogo and RNetLogo
2. Paper that uses NetLogo for ABM of infectious disease.
3. simstudy
4. epimodel, mathematical epidemiology models in R
5. HydeNet
6. bnlearn
7. bayesvl
8. Blog that introduces data simulation with probability distributions

Information on Bayesian statistics

While there are many blogs and online articles about Bayesian statistics and Bayesian networks, I found books were my best resource; enjoyable and worth the price tag.

1. Statistical Rethinking and related YouTube channel
2. Bayesian Networks with Examples in R
   
3. Bayesian Networks without Tears, a great introduction to the topic.