Directed acyclic graphs (DAGs), and causal graphs in general, provide a framework for making assumptions explicit and identifying confounders or mediators of the relationship between the exposure of interest and outcome that need to be adjusted for in analysis. Recently, I ran into the need to generate data from a DAG for a paper I am writing with my peers Kevin McIntyre and Joshua Wiener. After a quick Google search, I was pleasantly surprised to see there were several options to do so. In particular, the dagR library provides “functions to draw, manipulate, [and] evaluate directed acyclic graphs and simulate corresponding data”.

Besides `dagR`

‘s reference manual, a short letter published in Epidemiology, and a limited collection of examples, I couldn’t find too many resources regarding how to use the functionality provided by `dagR`

. The goal of this blog post is to provide an expository example of how to create a DAG and generate data from it using the `dagR`

library.

To simulate data from a DAG with `dagR`

, we need to:

- Create the DAG of interest using the
`dag.init`

function by specifying its nodes (exposure, outcome, and covariates) and their directed arcs (directed arrows to/from nodes). - Pass the DAG from (1) to the
`dag.sim`

function and specify the number of observations to be generated, arc coefficients, node types (binary or continuous), and parameters of the node distributions (Normal or Bernoulli).

For this tutorial, we are going to try to replicate the simple confounding/common cause DAG presented in Figure 1b as well as the more complex DAG in Figure 2a of Shier and Platt’s (2008) paper, *Reducing bias through directed acyclic graphs*.

```
library(dagR)
set.seed(12345)
```

## Simple confounding / common cause DAG

**Figure 1b of Shier and Platt (2008).**

### Specify

In dagR, the node corresponding to the primary exposure is assigned a ID of `0`

. The node corresponding to the outcome is assigned an ID of `-1`

. Any covariates are assigned an integer ID between `1`

and the total number of covariates `length(covs)`

. Covariates can be specified as either known or measured (`covs = 1`

), or unknown or unmeasured (`covs = 2`

).

Arcs are defined by provided the ID of the first node followed by the ID of the terminal node. Arcs are directed by default (`assocs = 0`

) but can be specified as undirected (`assocs = 1`

).

```
confound_dag <- dag.init(
# REQUIRED: SPECIFY COVARIATES AND THEIR TYPE
covs = c(1), # specify that single covariate (common cause) is "known"
# OPTIONAL: PROVIDE CUSTOM LABELS FOR DAG NODES AND SIM DATA
symbols = c('X', # symbol for primary exposure within DAG
'C', # " " common cause " "
'O'), # " " outcome " "
# REQUIRED: SPECIFY ARCS BETWEEN X, C, AND O
arcs = c(0, -1, # arc from primary exposure to outcome
1, 0, # arc from common cause to exposure
1, -1), # arc from common cause to outcome
# OPTIONAL: PROVIDE NAMES FOR X, C, AND O (USED IN LEGEND)
x.name = 'Primary exposure', # name of primary exposure
y.name = 'Outcome', # name of outcome
cov.names = c('Common cause'), # (vector of) names of covariates
# REQUIRED: TURN OFF WEIRD DEFAULTS
noxy = T # tell dagR we want to specify our *own* arc between X and O
)
```

That’s it! We’ve specified the simple confounding DAG. To make sure we have done so correctly, we can take a look at the details of our DAG object `confound`

as well as its summary `summary(confound)`

.

```
confound_dag
```

## $cov.types ## [1] 0 1 -1 ## ## $x ## [1] 0.0 0.5 1.0 ## ## $y ## [1] 0.0 0.7 0.0 ## ## $arc ## [,1] [,2] ## [1,] 1 3 ## [2,] 2 1 ## [3,] 2 3 ## ## $arc.type ## [1] 0 0 0 ## ## $curve.x ## [1] NA NA NA ## ## $curve.y ## [1] NA NA NA ## ## $xgap ## [1] 0.04 ## ## $ygap ## [1] 0.05 ## ## $len ## [1] 0.1 ## ## $names ## [1] "Primary exposure" "Common cause" "Outcome" ## ## $symbols ## [1] "X" "C" "O" ## ## $version ## [1] "1.2.0" ## ## Directed acyclic graph (DAG) object created by dagR version 1.2.0. ## For a convenient interpreted presentation, use summary().

`$cov.types`

summarizes tell variables in our model in their type, confirming that we have 1 primary exposure with ID`0`

, 1 covariate with ID`1`

, and 1 outcome with ID`-1`

.`$x`

and`$y`

provide the`(x, y)`

co-ordinates of nodes for plotting.`$arc`

confirms the arcs that we’ve drawn. Here, nodes are numbered such that node`1`

is the primary exposure, node`3`

is the outcome, and node`2`

is the common cause.`$arc.type`

confirms that all of our arcs are directed (`0`

), as desired.`$xgap`

,`$ygap`

, and`$len`

are graphical parameters corresponding to the arcs.- Finally, our
`$names`

and`$symbols`

appear to correctly correspond with our`$cov.types`

.

Perfect!

```
summary(confound_dag)
```

## DAG created by dagR version 1.2.0. ## - from X = Primary exposure to Y = Outcome ## - acyclic: TRUE ## - 3 nodes: X, Y, 1 cov(s) (0 unknown/unmeasured) ## - node 2 = Common cause ## - alternative node symbols: X C O ## - 3 arc(s) (0 undirected) ## - currently unadjusted ## - potentially biasing paths not sought ## - possible adjustment sets not examined

Our summary appears to confirm our interpretation of our DAG object.

### Plot

The `dag.draw`

function can be used to plot a DAG object.

```
confound_fig <- dag.draw(
confound_dag, # DAG object
numbering = F, # Number arcs: T/F
legend = T, # Legend: T/F
noxy = 2 # do not draw ? above X -> O arc
)
```

### Simulate

The algorithm used to simulate data is described in the `dagR`

documentation as follows:

*

Simulation steps: *

Simulate data for nodes

`i`

without ancestors, drawing from Normal distribution with mean`mu[i]`

and`stdev[i]`

(continuous node), or drawing from Bernoulli events with probability`mu[i]`

(binary node).Simulate data for nodes

`i`

for which all ancestors already have been simulated by multiplying the ancestor values with the corresponding arc coefficients and summing them up, shifting the resulting values to the mean`mu[i]`

specified for the currently simulated node (logit-transformed if binary), then adding noise drawn from a Normal distribution with mean 0 and standard deviation`stdev[i]`

, finally using the inverse logit of the resulting values as success probabilities for simulating binary data if node is binary.

Suppose we want to generate observations based on our DAG. In our generated sample, we want the primary exposure to be a binary variable with an observed prevalence of , the common cause to be a continuous variable distributed according to , and the outcome to be a continuous variable centered at .

The common cause node is our only node without ancestors. According to the algorithm above, data for will be generated first by drawing observations from a distribution.

The primary exposure only has one ancestor while the outcome has two. Thus, is the next candidate for simulation. is binary and is continuous so data is first generated on the logit scale such that,

where is the direct effect of and is the prevalence of . Then, the desired prevalence can be generated from using an operation similar to,

where the are random noise generated from . Finally, . In this simulation, and .

Finally, to generate centered at ,

where is the direct effect of , is the direct effect of , and the are random noise generated from . In this simulation, , , and .

```
confound_sim <- dag.sim(
# REQUIRED: PROVIDE DAG OBJECT AND DESIRED SAMPLE SIZE
confound_dag,
n = 100,
# REQUIRED: SPECIFY NODE TYPES (1 = BINARY, 0 = CTS)
binary = c(1, # binary primary exposure X
0, # continuous common cause C
0), # continuous outcome O
# REQUIRED: SPECIFY AND MEAN AND SD OF CTS NODES
# (*) MU = PREVALENCE FOR BINARY NODES
mu = c(0.5, # prevalence p of primary exposure X
120, # mean of common cause C
10), # " " outcome O
stdev = c(3, # sd of random noise added when generating X
20, # " " common cause C
1), # " " random noise added when generate O
# REQUIRED: SPECIFY DIRECT EFFECTS OF ARCS
b = c(5, # direct effect b of X -> O
log(2), # direct effect *exp(b)* of C -> X (since X binary)
1.5), # direct effect b of C -> O
# OPTIONAL: SET SEED FOR REPRODUCIBILITY
seed = 12345, # set seed for reproducibility
naming = 2, # use custom labels for output variable names
verbose = F # suppress alg. output
)
```

## Note: undirected arcs are discarded, an xy arc is appended. ## simulating node 2 from Normal, with mean SD 120 20 ## node 1 : calculated from nodes 2 via arcs 2 with exp(b) 2 , shifted to probability 0.5 (added noise SD: 3 ); ## node 3 : calculated from nodes 1 2 1 via arcs 1 3 4 with b 5 1.5 0 , shifted to mean 10 (added noise SD: 1 );

Let’s check out the data that gets spit out by `data.sim`

.

```
head(confound_sim)
```

## X C O ## 1 1 131.71058 21.701275 ## 2 1 134.18932 27.229184 ## 3 0 117.81393 -2.878963 ## 4 0 110.93006 -15.095081 ## 5 1 132.11775 22.855403 ## 6 0 83.64088 -53.952826

```
apply(confound_sim, MARGIN = 2, FUN = mean)
```

## X C O ## 0.54000 124.90394 10.07882

The data has our desired variable types and appears to have produced the correct means and prevalence.

Let’s see if we can use the generated values to derive and as outlined above.

```
# Logit and inverse logit functions for generating X
logit <- function(p){log(p/(1-p))}
ilogit <- function(x){1/(1+exp(-x))}
s_X = 3
s_O = 1
b_CX = log(10)
b_CO = 1.5
b_XO = 5
p_i = ilogit((b_CX * confound_sim$C - mean(b_CX * confound_sim$C) + logit(0.5)) + rnorm(n = 100, mean = 0, sd = s_X))
X_i = sapply(p_i, FUN = function(x) rbinom(n = 1, size = 1, prob = x))
mean(X_i)
```

## [1] 0.56

```
O_i = (b_CO * confound_sim$C + b_XO * X_i - mean(b_CO * confound_sim$C + b_XO * X_i) + 10) + rnorm(n = 100, mean = 0, sd = s_O)
mean(O_i)
```

## [1] 9.946751

Sweet, this is pretty close to what we got from our simulation.

Finally, we can construct some regression models using our simulate data to determine if they behave how we would expect. In a model of on , without adjustment for , we would expect the direct effect of to be distorted due to confounding by . After adjustment for the common cause , the direct effect of should be unbiased.

```
lm(O ~ X, data = confound_sim)
```

## ## Call: ## lm(formula = O ~ X, data = confound_sim) ## ## Coefficients: ## (Intercept) X ## -22.20 59.78

In the model adjusting for alone, it appears the direct effect of has been overestimated as a result of failing to adjust for .

```
lm(O ~ X + C, data = confound_sim)
```

## ## Call: ## lm(formula = O ~ X + C, data = confound_sim) ## ## Coefficients: ## (Intercept) X C ## -180.689 4.860 1.506

After adjustment for , we have approximately recovered our direct effect of which we set as .

## More complex DAG

**Figure 2a of Shier and Platt (2008).**

### Specify

```
complex_dag <- dag.init(
# REQUIRED: SPECIFY COVARIATES AND THEIR TYPE
covs = c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1),
# OPTIONAL: PROVIDE CUSTOM LABELS FOR DAG NODES AND SIM DATA
symbols = c('X', # 0: Warm-up exercises (primary exposure)
'IGP', # 1: Intra-game proprioception
'TMA', # 2: Team motivation, aggression
'PI', # 3: Previous injury
'CS', # 4: Contact sport
'C', # 5: Coach
'PGP', # 6: Pre-game proprioception
'FL', # 7: Fitness level
'G', # 8: Genetics
'NMF', # 9: Neuromuscular fatigue
'CTD', # 10: Connective Tissue Disorder
'TW', # 11: Tissue weakness
'O'), # -1: Injury (outcome)
# REQUIRED: SPECIFY ARCS BETWEEN X, C, AND O
arcs = c(0, 1, # X to IGP
1, -1, # IGP to O
2, 0, # TMA to X
2, 3, # TMA to PI
4, 3, # CS to PI
4, 1, # CS to IGP
6, 0, # PGP to X
5, 2, # C to TMA
5, 7, # C to FL
7, 6, # FL to PGP
7, 9, # FL to NMF
8, 7, # G to FL
8, 9, # G to NMF
8, 10, # G to CTD
10, 9, # CTD to NMF
10, 11, # CTD to TW
11, -1, # TW to O
9, 1, # NMF to IGP
9, -1, # NMF to O
1, -1), # arc from common cause to outcome
# OPTIONAL: PROVIDE NAMES FOR X, C, AND O (USED IN LEGEND)
x.name = 'Warm-up Exercises', # name of primary exposure
y.name = 'Injury', # name of outcome
cov.names = c('Intra-game proprioception',
'Team motivation, aggression',
'Previous injury',
'Contact sport',
'Coach',
'Pre-game proprioception',
'Fitness level',
'Genetics',
'Neuromuscular fatigue',
'Connective Tissue Disorder',
'Tissue weakness'),
# REQUIRED: TURN OFF WEIRD DEFAULTS
noxy = T # tell dagR we want to specify our *own* arc between X and O
)
```

```
complex_dag
```

## $cov.types ## [1] 0 1 1 1 1 1 1 1 1 1 1 1 -1 ## ## $x ## [1] 0.000000000 -0.176148078 -0.106217783 0.005025253 0.150000000 0.318826668 0.500000000 0.681173332 0.850000000 ## [10] 0.994974747 1.106217783 1.176148078 1.000000000 ## ## $y ## [1] 0.0000000 0.1811733 0.3500000 0.4949747 0.6062178 0.6761481 0.7000000 0.6761481 0.6062178 0.4949747 0.3500000 ## [12] 0.1811733 0.0000000 ## ## $arc ## [,1] [,2] ## [1,] 1 2 ## [2,] 2 13 ## [3,] 3 1 ## [4,] 3 4 ## [5,] 5 4 ## [6,] 5 2 ## [7,] 7 1 ## [8,] 6 3 ## [9,] 6 8 ## [10,] 8 7 ## [11,] 8 10 ## [12,] 9 8 ## [13,] 9 10 ## [14,] 9 11 ## [15,] 11 10 ## [16,] 11 12 ## [17,] 12 13 ## [18,] 10 2 ## [19,] 10 13 ## [20,] 2 13 ## ## $arc.type ## [1] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ## ## $curve.x ## [1] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA ## ## $curve.y ## [1] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA ## ## $xgap ## [1] 0.04 ## ## $ygap ## [1] 0.05 ## ## $len ## [1] 0.1 ## ## $names ## [1] "Warm-up Exercises" "Intra-game proprioception" "Team motivation, aggression" ## [4] "Previous injury" "Contact sport" "Coach" ## [7] "Pre-game proprioception" "Fitness level" "Genetics" ## [10] "Neuromuscular fatigue" "Connective Tissue Disorder" "Tissue weakness" ## [13] "Injury" ## ## $symbols ## [1] "X" "IGP" "TMA" "PI" "CS" "C" "PGP" "FL" "G" "NMF" "CTD" "TW" "O" ## ## $version ## [1] "1.2.0" ## ## Directed acyclic graph (DAG) object created by dagR version 1.2.0. ## For a convenient interpreted presentation, use summary().

```
summary(complex_dag)
```

## DAG created by dagR version 1.2.0. ## - from X = Warm-up Exercises to Y = Injury ## - acyclic: TRUE ## - 13 nodes: X, Y, 11 cov(s) (0 unknown/unmeasured) ## - node 2 = Intra-game proprioception ## - node 3 = Team motivation, aggression ## - node 4 = Previous injury ## - node 5 = Contact sport ## - node 6 = Coach ## - node 7 = Pre-game proprioception ## - node 8 = Fitness level ## - node 9 = Genetics ## - node 10 = Neuromuscular fatigue ## - node 11 = Connective Tissue Disorder ## - node 12 = Tissue weakness ## - alternative node symbols: X IGP TMA PI CS C PGP FL G NMF CTD TW O ## - 20 arc(s) (0 undirected) ## - currently unadjusted ## - potentially biasing paths not sought ## - possible adjustment sets not examined

### Plot

By default, `dagR`

arranges covariates in a “semi-circle” pattern. It doesn’t work too well for a large DAG such as this one.

```
complex_fig0 <- dag.draw(complex_dag)
```

We can provide custom co-ordinates to get our DAG to closer resemble the form of Figure 2a.

```
# Co-ordinates specified in the order of symbols
# "X" "IGP" "TMA" "PI" "CS" "C" "PGP" "FL" "G" "NMF" "CTD" "TW" "O"
complex_dag$x <- c(0, 1, 0, 0.5, 1.0, 0.5, 0.5, 1.0, 1.5, 1.5, 2.0, 2.0, 2)
complex_dag$y <- c(0.25, 0.25, 1, 0.5, 1.0, 1.5, 1.0, 1.5, 1.5, 1.0, 1.5, 0.5, 0.25)
```

```
complex_fig1 <- dag.draw(complex_dag, noxy = 2)
```

### Simulate

I’m going to randomly generate probabilities between for binary variables and means between for continuous variables. All standard deviations will be set to 1. All twenty direct effects will be generated from a uniform distribution between .

```
complex_sim <- dag.sim(
# REQUIRED: PROVIDE DAG OBJECT AND DESIRED SAMPLE SIZE
complex_dag,
n = 100,
# REQUIRED: SPECIFY NODE TYPES (1 = BINARY, 0 = CTS)
binary = c(1, # X
0, # IGP
0, # TMA
1, # PI
1, # CS
1, # C
0, # PGP
0, # FL
1, # G
1, # NMF
1, # CTD
1, # TW
1), # O
# REQUIRED: SPECIFY AND MEAN AND SD OF CTS NODES
# (*) MU = PREVALENCE FOR BINARY NODES
mu = c(runif(n = 1, min = 0.1, max = 0.9), # X
runif(n = 2, min = 10, max = 100), # IGP, TMA
runif(n = 3, min = 0.1, max = 0.9), # PI, CS, C
runif(n = 2, min = 10, max = 100), # PGP, FL
runif(n = 5, min = 0.1, max = 0.9)), # G, NMF, CTD, TW, O
stdev = rep(1, 13),
# REQUIRED: SPECIFY DIRECT EFFECTS OF 20 ARCS
b = runif(n = 20, min = -1, max = 1),
# OPTIONAL: SET SEED FOR REPRODUCIBILITY
seed = 12345, # set seed for reproducibility
naming = 2, # use custom labels for output variable names
verbose = F # suppress alg. output
)
```

## Note: undirected arcs are discarded, an xy arc is appended. ## simulating node 5 from Binomial, with prob 0.465184768103063 ## simulating node 6 from Binomial, with prob 0.23309742808342 ## simulating node 9 from Binomial, with prob 0.682164203003049 ## node 11 : calculated from nodes 9 via arcs 14 with exp(b) 0.836072462190722 , shifted to probability 0.127628348022699 (added noise SD: 1 ); ## node 12 : calculated from nodes 11 via arcs 16 with exp(b) 0.854790032790292 , shifted to probability 0.221898792125285 (added noise SD: 1 ); ## node 3 : calculated from nodes 6 via arcs 8 with b 0.524216595105827 , shifted to mean 78.4884095494635 (added noise SD: 1 ); ## node 4 : calculated from nodes 3 5 via arcs 4 5 with exp(b) 0.444515879520425 2.66794269837119 , shifted to probability 0.808899652957916 (added noise SD: 1 ); ## node 8 : calculated from nodes 6 9 via arcs 9 12 with b -0.336310201324522 0.907394020818174 , shifted to mean 55.8301902096719 (added noise SD: 1 ); ## node 10 : calculated from nodes 8 9 11 via arcs 11 13 15 with exp(b) 1.4709273414126 2.0516112877568 2.16855062119996 , shifted to probability 0.891789550334215 (added noise SD: 1 ); ## node 7 : calculated from nodes 8 via arcs 10 with b 0.697887991555035 , shifted to mean 39.2585848039016 (added noise SD: 1 ); ## node 1 : calculated from nodes 3 7 via arcs 3 7 with exp(b) 1.78283463041697 1.98592216989379 , shifted to probability 0.676723117008805 (added noise SD: 1 ); ## node 2 : calculated from nodes 1 5 10 via arcs 1 6 18 with b 0.732577715069056 0.576788076665252 0.703376833815128 , shifted to mean 88.8195873773657 (added noise SD: 1 ); ## node 13 : calculated from nodes 2 12 10 2 1 via arcs 2 17 19 20 21 with exp(b) 0.669151770139957 1.78079431432041 0.543818744665322 0.710830900554524 1 , shifted to probability 0.688547961972654 (added noise SD: 1 );

```
head(complex_sim)
```

## X IGP TMA PI CS C PGP FL G NMF CTD TW O ## 1 0 88.64873 78.96827 1 0 0 39.65966 58.32906 1 1 0 0 0 ## 2 1 89.27663 79.58407 0 0 1 40.24144 54.95490 0 1 0 0 0 ## 3 0 88.73670 78.42844 0 0 1 39.05994 55.80366 1 0 0 1 1 ## 4 0 87.55727 78.42953 1 0 1 38.62691 55.84037 1 1 0 0 0 ## 5 0 87.70700 78.87595 1 0 0 38.22488 55.41399 1 1 0 1 1 ## 6 1 89.25711 78.44596 0 0 0 39.28925 54.65930 1 1 0 0 1

Beauty!

Use this code for your own research by downloading this blogpost as an RMarkdown file here.