simulation-study-example/README.Rmd at main · ellessenne/simulation-study-example · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
---
output: github_document
editor_options:
  chunk_output_type: console
---

<!-- README.md is generated from README.Rmd. Please edit that file -->

```{r, include = FALSE}
library(knitr)
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  out.width = "70%",
  fig.width = 5,
  fig.height = 3,
  fig.align = "center",
  dpi = 600,
  message = FALSE,
  warning = FALSE
)
```

# Example Code for a Monte Carlo Simulation Study

Code in R and Stata to code, run, and analyse a simulation study is included in this repository.
The specific simulation study is described down below.

## Structure of this repository

This repository is organised as follows:

* R code is included in the folder named `R`, with three scripts: one for coding the simulation study, one for analysing the results, and one for creating plots and tables with the results;

* Stata code is included in the folder named `Stata`, including three analogous scripts.

* Simulated data with the results obtained by the above-mentioned code in R and Stata are included in the folder named `data`.
  The R dataset has extension `*.RDS`, while the Stata dataset has extension `*.DTA`.

## Aims of the simulation study

The aim of this simulation study is pedagogical.

Say we have a study which can be summarised by the following DAG:

```{r dag}
library(ggdag)
library(ggplot2)

confounder_triangle() %>%
  ggdag() +
  theme_dag_blank(base_size = 12) +
  coord_cartesian(xlim = c(-0.5, 2.5), ylim = c(-0.5, 1.5))
```

where X is the exposure, Y the outcome, and Z a confounder.

> Note that the exposure and the outcome are independent, there is no edge between X and Y.

We know that there should be no observed association between X and Y once we adjust for Z in our analysis.

However:

1. Is it true that once we adjust for Z we observe no association between X and Y?

1. What happens if we fail to adjust for a confounder in our analysis?

This is what we are trying to learn from this simulation study.

## Data-generating mechanisms

X is a binary treatment variable, Z is a continuous confounder, and Y is a continuous outcome.

Z is simulated from a standard normal distribution, while X depends on Z according to a logistic model:

$$
\log \frac{P(X = 1)}{P(X = 0)} = \gamma_0 + \gamma_1 Z
$$

Y depends on Z according to a linear model:

$$
Y = \alpha_0 + \alpha_1 Z + \varepsilon
$$

with $\varepsilon$ following a standard normal distribution.

We fix the following parameters:

* $\gamma_0 = 1$,

* $\gamma_1 = 3$,

* $\alpha_0 = 10$.

For $\alpha_1$, we actually simulate two scenarios:

1. $\alpha_1 = 5$: in this scenario, Z is actually a confounder;

1. $\alpha_1 = 0$: in the second scenario, we remove the edge between Z and Y. Here Z is not a counfounder anymore.

Finally, we generate N = 200 independent subjects per each simulated dataset.

## Estimand

We have a single estimand of interest here, the regression coefficient that quantifies the association between X and Y.

> Remember that, according to our data-generating mechanisms, we expect _no_ association between X and Y.

## Methods

We estimate the regression coefficient of interest using linear regression.
Specifically, we fit the following two models:

1. A model that fails to adjust for the observed confounder Z (model 1):

$$
Y = \beta_0 + \beta_1 X + \varepsilon
$$

1. A model that properly adjusts for Z (model 2):

$$
Y = \beta_0 + \beta_1 X + \beta_2 Z + \varepsilon
$$

The coefficient of interest is $\beta_1$ according to this notation.

## Performance measures

The key performance measure of interest is bias, defined as

$$
E(\hat{\beta_1})
$$

as we know that — according to our DGMs — the true $\beta_1 = 0$.

We might also be interested in estimating the power of a significance test for $\hat{\beta_1}$ at a given significance level, e.g. 0.05.
In our example, think of that as a test for the null hypothesis $\beta_1 = 0$.

# License

<a rel="license" href="http://creativecommons.org/licenses/by/4.0/">
    <img alt="Creative Commons License" style="border-width:0" src="https://i.creativecommons.org/l/by/4.0/88x31.png" />
</a>
<br />
This work is licensed under a <a rel="license" href="http://creativecommons.org/licenses/by/4.0/">Creative Commons Attribution 4.0 International License</a>.