forked from remlapmot/cibookex-r
-
Notifications
You must be signed in to change notification settings - Fork 0
/
16-iv-r.Rmd
93 lines (72 loc) · 2.81 KB
/
16-iv-r.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
# 16. Instrumental variables estimation{-}
```{r setup, include=FALSE}
knitr::opts_chunk$set(collapse = TRUE, comment = '#>')
```
## Program 16.1
- Estimating the average causal using the standard IV estimator via the calculation of sample averages
- Data from NHEFS
```{r, results='hide', message=FALSE, warning=FALSE}
library(here)
```
```{r}
#install.packages("readxl") # install package if required
library("readxl")
nhefs <- read_excel(here("data", "NHEFS.xls"))
# some preprocessing of the data
nhefs$cens <- ifelse(is.na(nhefs$wt82), 1, 0)
summary(nhefs$price82)
# for simplicity, ignore subjects with missing outcome or missing instrument
nhefs.iv <- nhefs[which(!is.na(nhefs$wt82) & !is.na(nhefs$price82)),]
nhefs.iv$highprice <- ifelse(nhefs.iv$price82>=1.5, 1, 0)
table(nhefs.iv$highprice, nhefs.iv$qsmk)
t.test(wt82_71 ~ highprice, data=nhefs.iv)
```
## Program 16.2
- Estimating the average causal effect using the standard IV estimator via two-stage-least-squares regression
- Data from NHEFS
```{r}
#install.packages ("sem") # install package if required
library(sem)
model1 <- tsls(wt82_71 ~ qsmk, ~ highprice, data = nhefs.iv)
summary(model1)
confint(model1) # note the wide confidence intervals
```
## Program 16.3
- Estimating the average causal using the standard IV estimator via additive marginal structural models
- Data from NHEFS
- G-estimation: Checking one possible value of psi
- See Chapter 14 for program that checks several values and computes 95% confidence intervals
```{r}
nhefs.iv$psi <- 2.396
nhefs.iv$Hpsi <- nhefs.iv$wt82_71-nhefs.iv$psi*nhefs.iv$qsmk
#install.packages("geepack") # install package if required
library("geepack")
g.est <- geeglm(highprice ~ Hpsi, data=nhefs.iv, id=seqn, family=binomial(),
corstr="independence")
summary(g.est)
beta <- coef(g.est)
SE <- coef(summary(g.est))[,2]
lcl <- beta-qnorm(0.975)*SE
ucl <- beta+qnorm(0.975)*SE
cbind(beta, lcl, ucl)
```
## Program 16.4
- Estimating the average causal using the standard IV estimator with altnerative proposed instruments
- Data from NHEFS
```{r}
summary(tsls(wt82_71 ~ qsmk, ~ ifelse(price82 >= 1.6, 1, 0), data = nhefs.iv))
summary(tsls(wt82_71 ~ qsmk, ~ ifelse(price82 >= 1.7, 1, 0), data = nhefs.iv))
summary(tsls(wt82_71 ~ qsmk, ~ ifelse(price82 >= 1.8, 1, 0), data = nhefs.iv))
summary(tsls(wt82_71 ~ qsmk, ~ ifelse(price82 >= 1.9, 1, 0), data = nhefs.iv))
```
## Program 16.5
- Estimating the average causal using the standard IV estimator
- Conditional on baseline covariates
- Data from NHEFS
```{r}
model2 <- tsls(wt82_71 ~ qsmk + sex + race + age + smokeintensity + smokeyrs +
as.factor(exercise) + as.factor(active) + wt71,
~ highprice + sex + race + age + smokeintensity + smokeyrs + as.factor(exercise) +
as.factor(active) + wt71, data = nhefs.iv)
summary(model2)
```