-
Notifications
You must be signed in to change notification settings - Fork 61
/
14-g-est-snms-stata.Rmd
170 lines (137 loc) · 4.43 KB
/
14-g-est-snms-stata.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
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
# 14. G-estimation of Structural Nested Models: Stata{-}
```{r, results='hide', message=FALSE, warning=FALSE}
library(Statamarkdown)
```
```
/***************************************************************
Stata code for Causal Inference: What If by Miguel Hernan & Jamie Robins
Date: 10/10/2019
Author: Eleanor Murray
For errors contact: [email protected]
***************************************************************/
```
## Program 14.1
- Ranks of extreme observations
- Data from NHEFS
- Section 14.4
```{stata}
/*For Stata 15 or later, first install the extremes function using this code:*/
* ssc install extremes
*Data preprocessing***
use ./data/nhefs, clear
gen byte cens = (wt82 == .)
/*Ranking of extreme observations*/
extremes wt82_71 seqn
/*Estimate unstabilized censoring weights for use in g-estimation models*/
glm cens qsmk sex race c.age##c.age ib(last).education ///
c.smokeintensity##c.smokeintensity c.smokeyrs##c.smokeyrs ///
ib(last).exercise ib(last).active c.wt71##c.wt71 ///
, family(binomial)
predict pr_cens
gen w_cens = 1/(1-pr_cens)
replace w_cens = . if cens == 1
/*observations with cens = 1 contribute to censoring models but not outcome model*/
summarize w_cens
/*Analyses restricted to N=1566*/
drop if wt82 == .
summarize wt82_71
save ./data/nhefs-wcens, replace
```
## Program 14.2
- G-estimation of a 1-parameter structural nested mean model
- Brute force search
- Data from NHEFS
- Section 14.5
```{stata}
use ./data/nhefs-wcens, clear
/*Generate test value of Psi = 3.446*/
gen psi = 3.446
/*Generate H(Psi) for each individual using test value of Psi and
their own values of weight change and smoking status*/
gen Hpsi = wt82_71 - psi * qsmk
/*Fit a model for smoking status, given confounders and H(Psi) value,
with censoring weights and display H(Psi) coefficient*/
logit qsmk sex race c.age##c.age ib(last).education ///
c.smokeintensity##c.smokeintensity c.smokeyrs##c.smokeyrs ///
ib(last).exercise ib(last).active c.wt71##c.wt71 Hpsi ///
[pw = w_cens], cluster(seqn)
di _b[Hpsi]
/*G-estimation*/
/*Checking multiple possible values of psi*/
cap noi drop psi Hpsi
local seq_start = 2
local seq_end = 5
local seq_by = 0.1 // Setting seq_by = 0.01 will yield the result 3.46
local seq_len = (`seq_end'-`seq_start')/`seq_by' + 1
matrix results = J(`seq_len', 4, 0)
qui gen psi = .
qui gen Hpsi = .
local j = 0
forvalues i = `seq_start'(`seq_by')`seq_end' {
local j = `j' + 1
qui replace psi = `i'
qui replace Hpsi = wt82_71 - psi * qsmk
quietly logit qsmk sex race c.age##c.age ///
ib(last).education c.smokeintensity##c.smokeintensity ///
c.smokeyrs##c.smokeyrs ib(last).exercise ib(last).active ///
c.wt71##c.wt71 Hpsi ///
[pw = w_cens], cluster(seqn)
matrix p_mat = r(table)
matrix p_mat = p_mat["pvalue","qsmk:Hpsi"]
local p = p_mat[1,1]
local b = _b[Hpsi]
di "coeff", %6.3f `b', "is generated from psi", %4.1f `i'
matrix results[`j',1]= `i'
matrix results[`j',2]= `b'
matrix results[`j',3]= abs(`b')
matrix results[`j',4]= `p'
}
matrix colnames results = "psi" "B(Hpsi)" "AbsB(Hpsi)" "pvalue"
mat li results
mata
res = st_matrix("results")
for(i=1; i<= rows(res); i++) {
if (res[i,3] == colmin(res[,3])) res[i,1]
}
end
* Setting seq_by = 0.01 will yield the result 3.46
```
## Program 14.3
- G-estimation for 2-parameter structural nested mean model
- Closed form estimator
- Data from NHEFS
- Section 14.6
```{stata}
use ./data/nhefs-wcens, clear
/*create weights*/
logit qsmk sex race c.age##c.age ib(last).education ///
c.smokeintensity##c.smokeintensity c.smokeyrs##c.smokeyrs ///
ib(last).exercise ib(last).active c.wt71##c.wt71 ///
[pw = w_cens], cluster(seqn)
predict pr_qsmk
summarize pr_qsmk
/* Closed form estimator linear mean models **/
* ssc install tomata
putmata *, replace
mata: diff = qsmk - pr_qsmk
mata: part1 = w_cens :* wt82_71 :* diff
mata: part2 = w_cens :* qsmk :* diff
mata: psi = sum(part1)/sum(part2)
/*** Closed form estimator for 2-parameter model **/
mata
diff = qsmk - pr_qsmk
diff2 = w_cens :* diff
lhs = J(2,2, 0)
lhs[1,1] = sum( qsmk :* diff2)
lhs[1,2] = sum( qsmk :* smokeintensity :* diff2 )
lhs[2,1] = sum( qsmk :* smokeintensity :* diff2)
lhs[2,2] = sum( qsmk :* smokeintensity :* smokeintensity :* diff2 )
rhs = J(2,1,0)
rhs[1] = sum(wt82_71 :* diff2 )
rhs[2] = sum(wt82_71 :* smokeintensity :* diff2 )
psi = (lusolve(lhs, rhs))'
psi
psi = (invsym(lhs'lhs)*lhs'rhs)'
psi
end
```