-
Notifications
You must be signed in to change notification settings - Fork 2
/
README.Rmd
141 lines (97 loc) · 7.41 KB
/
README.Rmd
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
---
output: github_document
---
<!-- README.md is generated from README.Rmd. Please edit that file -->
# tipr: R tools for tipping point sensitivity analyses <img src="man/figures/logo.png" align="right" height="138" />
<!-- badges: start -->
[](https://github.com/r-causal/tipr/actions)
[](https://doi.org/10.21105/joss.04495)
<!-- badges: end -->
**Authors:** [Lucy D'Agostino McGowan](https://www.lucymcgowan.com/), [Malcolm Barrett](https://malco.io/)<br/>
**License:** [MIT](https://opensource.org/license/mit/)
## Installation
Install the CRAN version
```{r, eval = FALSE}
install.packages("tipr")
```
Or install the development version from GitHub:
```{r, eval = FALSE}
# install.packages(devtools)
devtools::install_github("r-causal/tipr")
```
```{r}
library(tipr)
```
## Usage
After fitting your model, you can determine the unmeasured confounder needed to tip your analysis. This unmeasured confounder is determined by two quantities, the relationship between the exposure and the unmeasured confounder (if the unmeasured confounder is continuous, this is indicated with `exposure_confounder_effect`, if binary, with `exposed_confounder_prev` and `unexposed_confounder_prev`), and the relationship between the unmeasured confounder and outcome `confounder_outcome_effect`. Using this `r emo::ji("package")`, we can fix one of these and solve for the other. Alternatively, we can fix both and solve for `n`, that is, how many unmeasured confounders of this magnitude would tip the analysis.
This package comes with a few example data sets. For this example, we will use `exdata_rr`. This data set was simulated such that there are two confounders, one that was "measured" (and thus usable in the main analysis, this is called `measured_confounder`) and one that is "unmeasured" (we have access to it because this is simulated data, but ordinarily we would not, this variable is called `.unmeasured_confounder`).
Using the example data `exdata_rr`, we can estimate the exposure-outcome relationship using the measured confounder as follows:
```{r}
mod <- glm(outcome ~ exposure + measured_confounder, data = exdata_rr,
family = poisson)
mod |>
broom::tidy(exponentiate = TRUE, conf.int = TRUE)
```
We see the above example, the exposure-outcome relationship is 1.5 (95% CI: 1.1, 2.1). Note, in practice when estimating the effect of an exposure on a binary outcome using a GLM with the Poisson distribution and log link function, it is important to use a sandwich estimator to appropriately estimate the variability (this can be done in R using the `sandwich` package), which in this case gives a very similar result (95% CI: 1.1, 2.0).
## Continuous unmeasured confounder example
We are interested in a continuous unmeasured confounder, so we will use the `tip_with_continuous()` function.
Let's assume the unmeasured confounder is normally distributed with a mean of 0.5 in the exposed group and 0 in the unexposed (and unit variance in both), resulting in a mean difference of 0.5 (`exposure_confounder_effect = 0.5`), let's solve for the relationship between the unmeasured confounder and outcome needed to tip the analysis (in this case, we are solving for `confounder_outcome_effect`).
```{r}
tip(effect_observed = 1.5, exposure_confounder_effect = 0.5)
```
A hypothetical unobserved continuous confounder a scaled mean difference between exposure groups of `0.5` would need a relationship of at least 2.25 with the outcome to tip this analysis at the point estimate.
```{r}
tip(effect_observed = 1.09, exposure_confounder_effect = 0.5)
```
A hypothetical unobserved continuous confounder a scaled mean difference between exposure groups of `0.5` would need a relationship of at least 1.19 with the outcome to tip this analysis at the 5% level, rendering it inconclusive.
Because this is simulated data, we can see what the *true* unmeasured confounder looked like. First we will calculate the difference in scaled means.
```{r}
exdata_rr |>
dplyr::group_by(exposure) |>
dplyr::summarise(m = mean(.unmeasured_confounder / sd(.unmeasured_confounder))) |>
tidyr::pivot_wider(names_from = exposure,
values_from = m,
names_prefix = "u_") |>
dplyr::summarise(estimate = u_1 - u_0)
```
Now we can refit the above model with this unmeasured confounder included. According to our tipping point result, as long as the risk ratio of the unmeasured confounder and outcome in the model is greater than 2.25, the result that we observed will be "tipped" (the point estimate will cross the null).
```{r}
mod_true <- glm(
outcome ~ exposure + measured_confounder + .unmeasured_confounder,
data = exdata_rr,
family = poisson)
mod_true |>
broom::tidy(exponentiate = TRUE, conf.int = TRUE)
```
Notice here the `.unmeasured_confounder` effect is 2.42 (which is greater than the 2.25 we calculated that would be needed to render our result null) and, as expected, the point estimate for the `exposure` has crossed the null (and now is less than 1).
## Binary unmeasured confounder example
Now we are interested in the binary unmeasured confounder, so we will use the `tip_with_binary()` function.
Let's assume the unmeasured confounder is prevalent in 25% of the exposed population (`exposed_confounder_prev = 0.25`) and in 10% of the unexposed population (`unexposed_confounder_prev = 0.10`) -- let's solve for the relationship between the unmeasured confounder and the outcome needed to tip the analysis (`confounder_outcome_effect`).
```{r}
tip_with_binary(effect_observed = 1.09,
exposed_confounder_prev = 0.25,
unexposed_confounder_prev = 0.10)
```
A hypothetical unobserved binary confounder that is prevalent in 10% of the unexposed population and 25% of the exposed population would need to have a relationship with the outcome of 1.64 to tip this analysis at the 5% level, rendering it inconclusive.
## Many unmeasured confounders
Suppose we are concerned that there are many small, independent, continuous, unmeasured confounders present.
```{r}
tip(effect_observed = 1.09,
exposure_confounder_effect = 0.25,
confounder_outcome_effect = 1.05)
```
It would take about `7` independent standardized Normal unmeasured confounders with a mean difference between exposure groups of 0.25 and a relationship with the outcome of 1.05 tip the observed analysis at the 5% level, rendering it inconclusive.
## Integration with broom
These functions were created to easily integrate with models tidied using the **broom** package. This is not _necessary_ to use these functions, but a nice feature if you choose to do so. Here is an example of a logistic regression fit with `glm` and tidied with the `tidy` function **broom** that can be directly fed into the `tip()` function.
```{r}
if (requireNamespace("broom", quietly = TRUE) && requireNamespace("dplyr", quietly = TRUE)) {
glm(outcome ~ exposure + measured_confounder, data = exdata_rr,
family = poisson) |>
broom::tidy(conf.int = TRUE, exponentiate = TRUE) |>
dplyr::filter(term == "exposure") |>
dplyr::pull(conf.low) |>
tip(confounder_outcome_effect = 2.5)
}
```
## Code of Conduct
Please note that the tipr project is released with a [Contributor Code of Conduct](https://contributor-covenant.org/version/2/0/CODE_OF_CONDUCT.html). By contributing to this project, you agree to abide by its terms.