Bayes' Town: A place of data simulation

Welcome to Bayes’ Town

Bayes’ Town is a special place, a place where we can be omnipotent and omniscient. We take comfort in having knowledge and control of all things, even though we know the entire place is apart from reality, a muddled reflection at best and a complete fantasy at worst. Bayes’ Town is a simulation, a useful tool to explore scenarios and test hypotheses, based in reality or otherwise.

For public health, creating a synthetic population can be a valuable tool. Population characteristics, and relations between them, can be modeled based upon prior knowledge or to understand ‘what-if’ scenarios. Sometimes, it even is just a good way to better understand how a model is working. If you know the generative process of a synthetic population then you can see how well certain models perform under those circumstances. More on this later…

Simulating Bayes’ Town’s Blight

As we are the masters of our simulated domain, we can decide exactly how the characteristics of our Bayes’ Town citizens are connected. Citizens of Bayes’ Town are very concerned about disease, specifically the rampant Disease X. Fortunately for us, as the creators we can instantly summon the image of how the population is impacted by this pestilence. We can visualize the relationship as a network between variables of sex, age, and whether or not the individual lives in an urban or rural location. Using DiagrammeR and the DOT language, we can visualize this relationship.

Figure 1: The “known” relationship between Disease X and the citizens of Bayes’ Town.

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

Full simulation run

At the beginning a model and parameters were already determined. Below I provide the code to create the simulated Bayes’ Town. A few helper functions built on top of rjags make the process a bit more friendly to read. The Bayes’ Town example is an expert system using just priors but if we decide to include data in the model we simply add a loop to the likelihood function.

# Set input parameters
inputParam <- list(p_sex = 0.52,
age_shape = 6,
age_scale = 6,
age_param = 0.7)

# Define model
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);
}')

# Simulate new town
newTown <- BayesTown$new(population_size = 1000, jags_data = inputParam, n.chains = 2, n.iter = 5000, jags_model = simModel) # Keep only subset of parameters of interest newTown$set_variables(c('age', 'sex', 'urbanRural', 'disease'), mode = 'include')

# Resample data
newTown\$resample()

We can now use newTown as we saw at the start of the post and explore our newly simulated dataset!

Resources

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.
Allen O'Brien
Infectious Disease Epidemiologist

I am an epidemiologist with a passion for teaching and all things data.