The hardware and bandwidth for this mirror is donated by dogado GmbH, the Webhosting and Full Service-Cloud Provider. Check out our Wordpress Tutorial.
If you wish to report a bug, or if you are interested in having us mirror your free-software or open-source project, please feel free to contact us at mirror[@]dogado.de.
In this small vignette, we introduce the sim_from_dag()
function, which can be used to simulate complex data from arbitrary
causal directed acyclic graphs (DAGs). The simulated data may include
continuous, binary, categorical, count or time-to-event variables. This
function is most useful if the DAG is static, meaning that there are no
time-varying variables. It is theoretically possible to use this
function to simulate data from DAGs with a time structure as well, but
there are some difficulties associated with it that will be discussed
later.
A causal DAG is a DAG in which all nodes correspond to variables and the directed edges correspond to direct causal relationships between these variables. A direct edge from node \(A\) to node \(B\) implies that there is direct causal effect of \(A\) on \(B\). On the other hand, if there is no edge from node \(A\) to node \(B\), there is no direct causal relationship between these variables. Using a DAG in this way makes it easy to encode the causal structure of a given system, which is very useful for causal inference. This general idea is a centerpiece of the structural approach to causality developed by Pearl (2002) and Spirtes et al. (1993). We strongly encourage the reader to make themselves familiar with some of this literature before moving on.
It is very simple to generate data from a defined causal DAG. To see why we first need to introduce the concept of root nodes and child nodes. A root node is a node in a DAG that does not have any edges pointing to it (no incoming arrows). A child node on the other hand is a node that has at least one incoming edge. In other words, root nodes have no direct causes but child nodes do. Every node pointing into another node is considered a parent of that child node. For example, consider the DAG in figure 1.
A small DAG with four nodes
Nodes \(A\) and \(B\) are root nodes because they do not have any directed edges pointing into them. Nodes \(C\) and \(D\) on the other hand are child nodes. The parents of node \(C\) are both \(A\) and \(B\) because both of these nodes have directed edge towards \(C\). Note that node \(B\) is not a parent of node \(D\) because there is no edge from \(B\) to \(D\).
As the name implies, DAGs do not have cycles. Therefore every DAG has at least one root node. Generating data for these nodes is the first step to simulate data for the whole DAG. Since root nodes have no parents, we can simply generate random data from them using an appropriate distribution. Once we have data for all root nodes, we can generate their directly connected child nodes next as a function of the root nodes (and perhaps additional random error). These direct child nodes are then used as input for the next child nodes in line and so on. This continues until every node has been generated. Since every DAG can be topologically sorted (Chickering 1995), this will always work. All we need is to specify the DAG and the functional relationship between each node and its parents.
Because the sim_from_dag()
function uses the method
described above, it requires information about the causal structure and
the exact form of the relationship between child nodes and their
parents. All of this information has to be included in the
dag
argument, which should be a DAG
object
created using the empty_dag()
function and grown using
node()
calls as described below. This can be done
completely manually (which is the usual strategy when conducting
simulation studies) or (partially) using existing data (which may be
useful when the interest is in getting a toy data set resembling real
data as closely as possible).
Regardless of which strategy you want to use, first you have to
initialize an empty DAG
object like this:
Afterwards you can add an unlimited amount of root nodes and child nodes to it. Multiple different types are implemented.
The values for the root_nodes
are simply sampled from
some defined distributions. Therefore, any function that generates
random draws from some distribution may be used here. Popular
alternatives for continuous data are the normal-, beta-,
gamma-distributions which are implemented in base R inside the
rnorm()
, rbeta()
and rgamma()
functions. For binary or categorical data we could use the custom
functions rbernoulli()
or rcategorical()
instead.
The simDAG
package implements the following types of
child_nodes
directly:
node_gaussian
: A node based on linear regression
(continuous data).node_binomial
: A node based on logistic regression
(binary data).node_multinomial
: A node based on multinomial logistic
regression (categorical data).node_poisson
: A node based on poisson regression (count
data).node_negative_binomial
: A node based on negative
binomial regression (count data).node_cox
: A node based on cox regression (time-to-event
data).node_conditional_prob
: A node based on conditional
probabilities (binary / categorical data).node_conditional_distr
: A node based on conditional
distributions (any data type).All of these nodes have their own documentation page containing a
detailed description on how data is generated from them. Although this
collection of nodes covers a lot of data types, it is still a somewhat
limited collection. If, for example, we wanted to add a child node that
is normally distributed but also truncated at specific values, we could
not do this using just the offered node functions. For this reason, the
sim_from_dag()
function also allows the user to use custom
functions as nodes, which makes it possible to model any kind of data
and any kind of relationship.
Suppose that node \(A\) in the
figure above stands for age
, \(B\) stands for sex
, \(C\) stands for the Body-Mass-Index
(BMI
) and \(D\) stands for
death
. We have to start by defining what the root nodes
should look like. We use the following code to define age
and sex
:
All nodes are defined by calling the node()
function and
adding the output to the dag
object using a simple
+
. This syntax is heavily inspired by the
simCausal
R-package (Sofrygin et al. 2017). Here, we assume
that age
is a continuous normally distributed variable with
a mean of 50 and a standard deviation of 4 (If this was a real
simulation study we would probably use a truncated normal distribution
to ensure that age is not negative). This can be done by setting the
dist
parameter to "rnorm"
, which is the
standard R function for generating random values from a normal
distribution. All arguments listed in the params
parameter
will be passed to this function. Similarly, we define sex
to be a Bernoulli distributed variable (taking only the values 0/1). We
assume that there is an even gender distribution by setting
p = 0.5
.
Next, we have to define what the relationship between the child nodes and their parents should look like. We may use the following code:
dag <- dag +
node("bmi", type="gaussian", parents=c("sex", "age"), betas=c(1.1, 0.4),
intercept=12, error=2) +
node("death", type="binomial", parents=c("age", "bmi"), betas=c(0.1, 0.3),
intercept=-15)
Since the bmi
node is dependent on both sex
and age
, we have to list both of these nodes as the parents
of bmi
. We then specify that the bmi
should be
a continuous variable modeled using a linear regression by setting
type="gaussian"
. The concrete regression equation is
defined through the use of the intercept
,
betas
and error
arguments. Our specification
for the bmi
node corresponds to the following equation:
\[ bmi = 12 + sex \cdot 1.1 + age \cdot 0.4 + N(0, 2), \]
where \(N(0, 2)\) indicates that the error term is modelled as a normally distributed variable with mean 0 and a standard deviation of 2.
Since death
has only two states (alive vs. dead), we use
a logistic regression model here instead. We can do this easily by
setting type="binomial"
. The rest of the syntax essentially
stays the same. The regression equation for death
as
described by the code above is then:
\[ logit(death) = -15 + age \cdot 0.1 + bmi \cdot 0.3. \]
To check whether we got the causal relationships right, we can call
the plot()
function on the DAG object. The output should
look very similar to the hand-drawn DAG above.
We can also directly print the underlying structural equations using
the summary()
function:
summary(dag)
#> A DAG object using the following structural equations:
#>
#> age ~ N(50, 4)
#> sex ~ Bernoulli(0.5)
#> bmi ~ N(12 + 1.1*sex + 0.4*age, 2)
#> death ~ Bernoulli(logit(-15 + 0.1*age + 0.3*bmi))
This is all correct. We can now use this DAG
object to
generate random data using the sim_from_dag()
function:
Setting a seed for the random number generator is necessary to obtain replicable results. The data generated using this code looks like this:
head(sim_dat, 5)
#> age sex bmi death
#> <num> <lgcl> <num> <lgcl>
#> 1: 55.48383 TRUE 33.25311 FALSE
#> 2: 47.74121 FALSE 29.58815 FALSE
#> 3: 51.45251 FALSE 30.12967 FALSE
#> 4: 52.53145 TRUE 32.07877 FALSE
#> 5: 51.61707 FALSE 36.09082 TRUE
Binary variables such as sex
and death
are
by default treated as logical variables, because this is the most memory
efficient way to store them. We can now check the distributions and
relationships in this dataset to confirm that it indeed corresponds to
our specified causal DAG. Starting with the root nodes:
This seems to be correct. Note that this is a finite dataset, which
means that the results will never exactly match the theoretical
distributions. But it’s definitely close enough here. To check if the
child nodes were modeled correctly, we simply fit the corresponding
models using the glm()
function:
mod_bmi <- glm(bmi ~ age + sex, data=sim_dat, family="gaussian")
summary(mod_bmi)
#>
#> Call:
#> glm(formula = bmi ~ age + sex, family = "gaussian", data = sim_dat)
#>
#> Deviance Residuals:
#> Min 1Q Median 3Q Max
#> -7.4923 -1.3951 0.0158 1.3831 8.0821
#>
#> Coefficients:
#> Estimate Std. Error t value Pr(>|t|)
#> (Intercept) 12.184254 0.253188 48.12 <2e-16 ***
#> age 0.396020 0.005039 78.58 <2e-16 ***
#> sexTRUE 1.177159 0.040563 29.02 <2e-16 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> (Dispersion parameter for gaussian family taken to be 4.112615)
#>
#> Null deviance: 70134 on 9999 degrees of freedom
#> Residual deviance: 41114 on 9997 degrees of freedom
#> AIC: 42524
#>
#> Number of Fisher Scoring iterations: 2
mod_death <- glm(death ~ age + bmi, data=sim_dat, family="binomial")
summary(mod_death)
#>
#> Call:
#> glm(formula = death ~ age + bmi, family = "binomial", data = sim_dat)
#>
#> Deviance Residuals:
#> Min 1Q Median 3Q Max
#> -2.4957 -0.9638 -0.5215 1.0239 2.5186
#>
#> Coefficients:
#> Estimate Std. Error z value Pr(>|z|)
#> (Intercept) -14.319144 0.371511 -38.54 <2e-16 ***
#> age 0.093229 0.007057 13.21 <2e-16 ***
#> bmi 0.289714 0.011331 25.57 <2e-16 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> (Dispersion parameter for binomial family taken to be 1)
#>
#> Null deviance: 13773 on 9999 degrees of freedom
#> Residual deviance: 11805 on 9997 degrees of freedom
#> AIC: 11811
#>
#> Number of Fisher Scoring iterations: 3
Evidently, the coefficients do match the causal coefficients we specified earlier.
If the data should resemble a specific real data set, it makes sense
to base the values for the causal coefficients on that specific data
set. This can be done by fitting a single model for each child node,
extracting the estimated coefficients from the fitted models and putting
those into an appropriate DAG
object. If the assumed DAG is
big, this can be a time-extensive task. The dag_from_data()
function automates this process. This function takes a node list
containing only minimal information about the causal structure and node
type and outputs a fully specified DAG
object.
For example, lets assume that the data we just generated
(sim_dat
) was our data set of interest. Let us also assume
that we know the true underlying causal diagram and have a rough idea
about the nature of the relationship between the nodes (e.g. we know or
can reasonably guess the node type). Now all we have to do is create a
partially specified DAG
in accordance to these assumptions
first:
dag <- empty_dag() +
node("age", type="rnorm") +
node("sex", type="rbernoulli") +
node("bmi", type="gaussian", parents=c("sex", "age")) +
node("death", type="binomial", parents=c("age", "bmi"))
This looks a lot like the code used above, except that we are not
explicitly defining the actual beta coefficients. We only define the
causal structure and the node types. Now we can call the
dag_from_data()
function:
It returns an object that includes a fully specified
DAG
, which can be used directly in the
sim_from_dag()
function:
The dag_from_data()
function essentially just fits the
corresponding models one by one for each node and extracts the relevant
data from the models to fill in the gaps in the empty nodes. If we set
return_models
to TRUE
in the
dag_from_data()
function call above, we can actually see
that it used the exact same models we fit earlier to check if the
simulation was valid.
Most real data sets include time-varying covariates, e.g. variables
that are measured at multiple points in time that are subject to
changes. It is possible to generate this type of data using the
sim_from_dag()
function as well. All we need to do is to
define an appropriate DAG
that directly specifies how the
variables change over time. For example, we can extend the simple
DAG
from above to include a dimension of time:
A small DAG with four nodes
Here, nodes \(A\) and \(B\) are time-constant variables that only
have a causal effect on the initial state of \(C\) and \(D\), while nodes \(C\) and \(D\) change over time interdependently. If
we want to simulate data from a DAG that looks like this using the
sim_from_dag()
function, we have to add a node to the for
every point in time that we want to consider.
We will quickly go through a somewhat simpler example, considering only 2 points in time. We define our nodes in the following way:
dag <- empty_dag() +
node("age", type="rnorm", mean=50, sd=4) +
node("sex", type="rbernoulli", p=0.5) +
node("bmi_t1", type="gaussian", betas=c(1.1, 0.4), parents=c("sex", "age"),
intercept=12, error=2) +
node("death_t1", type="binomial", parents=c("age", "sex", "bmi_t1"),
betas=c(0.1, 0.3, 0.1), intercept=-15) +
node("bmi_t2", type="gaussian", parents="bmi_t1", betas=c(1.1), intercept=0,
error=2) +
node("death_t2", type="binomial", betas=c(0.1, 0.3),
parents=c("age", "bmi_t2"), intercept=-15)
sim_dat <- sim_from_dag(dag=dag, n_sim=10000)
In this example, the bmi
at \(t = 1\) is a function of both
sex
and age
, but the bmi
at \(t = 2\) is only a function of the previous
bmi
. The death
node is determined by the
initial age
and by the time-varying bmi
. This
surely is not the most realistic example. It is only meant to show how
the sim_from_dag()
function may be used to incorporate
time-dependent covariates. If many points in time should be considered
or there are very complex time-dependent structures that may not be
easily described using a DAG like the one above, the
sim_discrete_time()
function also included in this package
may be used instead.
Judea Pearl (2009). Causality: Models, Reasoning and Inference. 2nd ed. Cambridge: Cambridge University Press
Peter Spirtes, Clark Glymour, and Richard Scheines (2000) Causation, Prediction, and Search. 2nd ed. The MIT Press, Cambridge
Chickering, D.M. (1995). A transformational characterization of equivalent Bayesian network structures. Proceedings of the 11th Conference on Uncertainty in Artificial Intelligence, Montreal, Canada, 87-98.
Oleg Sofrygin, Mark J. van der Laan, and Romain Neugebauer (2017). simcausal R Package: Conducting Transparent and Reproducible Simulation Studies of Causal Effect Estimation with Complex Longitudinal Data. In: Journal of Statistical Software. 81.2, pp. 1-47
These binaries (installable software) and packages are in development.
They may not be fully stable and should be used with caution. We make no claims about them.
Health stats visible at Monitor.