Targeted Inference.
Targeted Inference in R: targeted
Various methods for targeted and semiparametric inference including augmented inverse probability weighted estimators for missing data and causal inference (Bang and Robins (2005) <10.1111/j.1541-0420.2005.00377.x>) and estimators for risk differences and relative risks (Richardson et al. (2017) <10.1080/01621459.2016.1192546>).
Installation
You can install the released version of targeted from CRAN with:
install.packages("targeted")
And the development version from GitHub with:
# install.packages("devtools")
remotes::install_github("kkholst/targeted")
Examples
Targeted risk regression
library(targeted)
#> Loading required package: lava
library(magrittr)
Simulate some data:
m <- lvm() %>%
regression(a ~ x+z) %>%
regression(lp.target ~ 1) %>%
regression(lp.nuisance ~ x + z) %>%
distribution('a', binomial.lvm("logit")) %>%
binomial.rr('y', 'a', 'lp.target', 'lp.nuisance')
par <- c('a'=-2, 'lp.target'=1, 'lp.nuisance'=-1, 'lp.nuisance~x'=2)
d <- lava::sim(m, n=1e4, seed=1, p=par) %>%
subset(select=c('y', 'a','x','z'))
head(d)
#> y a x z
#> 1 0 0 -0.6264538 -0.8043316
#> 2 0 0 0.1836433 -1.0565257
#> 3 0 0 -0.8356286 -1.0353958
#> 4 0 0 1.5952808 -1.1855604
#> 5 1 0 0.3295078 -0.5004395
#> 6 0 0 -0.8204684 -0.5249887
fit <- riskreg(y ~ a, nuisance=~x+z, data=d, type="rr")
fit
#> Estimate Std.Err 2.5% 97.5% P-value
#> (Intercept) 0.9722 0.02896 0.9155 1.029 4.281e-247
Here the same design matrix is used for both the propensity model and the nuisance parameter (odds-product) model
summary(fit)
#> riskreg(formula = y ~ a, nuisance = ~x + z, data = d, type = "rr")
#>
#> Relative risk model
#> Response: y
#> Exposure: a
#>
#> Estimate Std.Err 2.5% 97.5% P-value
#> log(RR):
#> (Intercept) 0.9722 0.02896 0.9155 1.0290 4.281e-247
#> log(OP):
#> (Intercept) -0.9636 0.06603 -1.0930 -0.8342 3.072e-48
#> x 2.0549 0.07901 1.9000 2.2098 4.182e-149
#> z 1.0329 0.06728 0.9010 1.1648 3.468e-53
#> logit(Pr):
#> (Intercept) -1.9753 0.05631 -2.0856 -1.8649 1.284e-269
#> x 0.9484 0.04186 0.8664 1.0305 1.235e-113
#> z 1.0336 0.04878 0.9380 1.1292 1.187e-99
Double-robustness illustrated by using a wrong propensity model but a correct nuisance paramter (odds-product) model:
riskreg(y ~ a, nuisance=~x+z, propensity=~z, data=d, type="rr")
#> Estimate Std.Err 2.5% 97.5% P-value
#> (Intercept) 0.9709 0.02893 0.9142 1.028 7.053e-247
Or vice-versa
riskreg(y ~ a, nuisance=~z, propensity=~x+z, data=d, type="rr")
#> Estimate Std.Err 2.5% 97.5% P-value
#> (Intercept) 0.9931 0.03597 0.9226 1.064 8.286e-168
Whereas the MLE yields a biased estimate of the relative risk:
fit_mle <- with(d, riskreg_mle(y, a, x1=model.matrix(~1,d), x2=model.matrix(~z, d)))
estimate(fit_mle, 1)
#> Estimate Std.Err 2.5% 97.5% P-value
#> [p1] 1.289 0.02855 1.233 1.345 0
#>
#> Null Hypothesis:
#> [p1] = 0
To obtain an estimate of the risk-difference (here wrong model) we simply chance the type
argument
riskreg(y ~ a, nuisance=~x+z, data=d, type="rd")
#> Estimate Std.Err 2.5% 97.5% P-value
#> (Intercept) 0.5102 0.01613 0.4786 0.5418 1.135e-219
Interactions with the exposure can be examined with the target
argument
riskreg(y ~ a, target=a~x, nuisance=~x+z, data=d, type="rr")
#> Estimate Std.Err 2.5% 97.5% P-value
#> (Intercept) 1.0241 0.03659 0.9524 1.09584 1.986e-172
#> x -0.0825 0.03469 -0.1505 -0.01451 1.739e-02