-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathreg_1_table_3_4_5.do
290 lines (219 loc) · 9.17 KB
/
reg_1_table_3_4_5.do
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
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
//Neil Davies 03/07/15
//This creates an indicator for birth month in the UK Biobank data
cap prog drop reduced_form
prog def reduced_form
use "working data/cleaned_biobank_outcomes_ENGLISH" if male!=`1',clear
//Run baseline regressions just using actual schooling (more than 15 years and eduyears):
reg out_phys_v_act more_edu_15 i.mob cov_male if bw12!=. [pweight=weight], cluster(mobi)
regsave more_edu_15 using "results/outcome_act_exposure_`2'", replace detail(all) pval ci
reg out_phys_v_act eduyears i.mob cov_male if bw12!=. [pweight=weight], cluster(mobi)
regsave eduyears using "results/outcome_act_exposure2_`2'", replace detail(all) pval ci
//On the policy change and the two negative control populations
foreach j in bw12 N1_bw12 N2_bw12{
reg out_phys_v_act `j' i.mob cov_male if more_edu_15!=. [pweight=weight], cluster(mobi)
regsave `j' using "results/outcome_`j'_`2'", replace detail(all) pval ci
}
//Instrumental variable regressions
ivreg2 out_phys_v_act (more_edu_15 =post_reform) imob_* cov_male if bw12!=. [pweight=weight], endog(more_edu_15) cluster(mobi) partial(imob_* cov_male)
regsave more_edu_15 using "results/outcome_IV_`2'", replace detail(all) pval ci
ivreg2 out_phys_v_act (eduyears =post_reform) imob_* cov_male if bw12!=. [pweight=weight], endog(eduyears) cluster(mobi) partial(imob_* cov_male)
regsave eduyears using "results/outcome_IV2_`2'", replace detail(all) pval ci
//Repeat these regression for all the other outcomes:
ds out_*
foreach i in `r(varlist)'{
reg `i' more_edu_15 i.mob cov_male if bw12!=. [pweight=weight], cluster(mobi)
regsave more_edu_15 using "results/outcome_act_exposure_`2'", append detail(all) pval ci
reg `i' eduyears i.mob cov_male if bw12!=. [pweight=weight], cluster(mobi)
regsave eduyears using "results/outcome_act_exposure2_`2'", append detail(all) pval ci
foreach j in bw12 N1_bw12 N2_bw12{
reg `i' `j' i.mob cov_male if more_edu_15!=. & yob>1947 [pweight=weight], cluster(mobi)
regsave `j' using "results/outcome_`j'_`2'", append detail(all) pval ci
}
ivreg2 `i' (more_edu_15 =post_reform) imob_* cov_male if bw12!=. [pweight=weight] , endog(more_edu_15) cluster(mobi) partial(imob_* cov_male)
regsave more_edu_15 using "results/outcome_IV_`2'", append detail(all) pval ci
ivreg2 `i' (eduyears =post_reform) imob_* cov_male if bw12!=. [pweight=weight], endog(eduyears) cluster(mobi) partial(imob_* cov_male)
regsave eduyears using "results/outcome_IV2_`2'", append detail(all) pval ci
}
end
reduced_form 2 all
reduced_form 1 female
reduced_form 0 male
//Clean results:
cap prog drop clean
prog def clean
use "`2'",clear
keep depvar N coef ci_lower ci_upper pval
gen n=_n
order depvar N coef ci_lower ci_upper pval
drop if _n==_N
foreach j in N coef ci_lower ci_upper pval{
rename `j' `j'_`1'
}
save "working data/temp`1'.dta",replace
rm "`2'.dta"
end
foreach j in male female all{
clean 1 "results/outcome_act_exposure_`j'"
clean 2 "results/outcome_bw12_`j'"
clean 3 "results/outcome_N1_bw12_`j'"
clean 4 "results/outcome_N2_bw12_`j'"
clean 5 "results/outcome_act_exposure2_`j'"
use "working data/temp1.dta",clear
rm "working data/temp1.dta"
forvalues i=2(1)5{
joinby n using "working data/temp`i'.dta", unmatched(both)
tab _m
drop _m
rm "working data/temp`i'.dta"
}
drop n
sort n
drop n
save "results/table3_`j'_reduced_form",replace
}
use "results/table3_all_reduced_form"
append using "results/table3_male_reduced_form"
append using "results/table3_female_reduced_form"
save "results/table3",replace
use "results/table3",clear
gen analysis="all" if _n<26
replace analysis="males" if analysis=="" & _n<51
replace analysis="females" if analysis==""
joinby depvar using "working data/order",unmatched(master)
sort analysis n
rm "results/table3_all_reduced_form.dta"
rm "results/table3_male_reduced_form.dta"
rm "results/table3_female_reduced_form.dta"
//Clean IV results for Table 4
cap prog drop clean_iv
prog def clean_iv
use "`2'",clear
keep depvar N coef ci_lower ci_upper pval estatp cdf
gen n=_n
order depvar N coef ci_lower ci_upper pval
drop if _n==_N
save "working data/temp`1'.dta",replace
rm "`2'.dta"
end
clean_iv 1 "results/outcome_IV_all"
clean_iv 2 "results/outcome_IV_male"
clean_iv 3 "results/outcome_IV_female"
use "working data/temp1.dta",clear
append using "working data/temp2"
append using "working data/temp3"
order depvar N coef ci_lower ci_upper pval estatp cdf
save "results/table4",replace
use "results/table4",clear
drop n
gen analysis="all" if _n<26
replace analysis="males" if analysis=="" & _n<51
replace analysis="females" if analysis==""
joinby depvar using "working data/varname",unmatched(master)
drop _m
joinby depvar using "working data/order",unmatched(master)
sort analysis n
//Create indicator for binary or continious measure:
#delimit ;
gen binary=(depvar=="out_highbloodpressure"
|depvar=="out_diabetes"
|depvar=="out_stroke"
|depvar=="out_heartattack"
|depvar=="out_depression"
|depvar=="out_cancer"
|depvar=="out_dead"
|depvar=="out_exsmoker"
|depvar=="out_smoker"
|depvar=="out_income_under_18k"
|depvar=="out_income_over_31k"
|depvar=="out_income_over_52k"
|depvar=="out_income_over_100k");
#delimit cr
//Create Figure 5
mkmat coef ci_lower ci_upper if analysis=="all" & binary==1, mat(coef1_bin) rownames(depvar)
mkmat coef ci_lower ci_upper if analysis=="all" & binary==0, mat(coef1_con) rownames(depvar)
use "results/table3",clear
gen analysis="all" if _n<26
replace analysis="males" if analysis=="" & _n<51
replace analysis="females" if analysis==""
joinby depvar using "working data/order",unmatched(master)
drop _m
joinby depvar using "working data/varname",unmatched(master)
sort analysis n
//Create indicator for binary or continious measure:
#delimit ;
gen binary=(depvar=="out_highbloodpressure"
|depvar=="out_diabetes"
|depvar=="out_stroke"
|depvar=="out_heartattack"
|depvar=="out_depression"
|depvar=="out_cancer"
|depvar=="out_dead"
|depvar=="out_exsmoker"
|depvar=="out_smoker"
|depvar=="out_income_under_18k"
|depvar=="out_income_over_31k"
|depvar=="out_income_over_52k"
|depvar=="out_income_over_100k");
#delimit cr
mkmat coef_1 ci_lower_1 ci_upper_1 if analysis=="all" & binary==1, mat(coef2_bin) rownames(depvar)
mkmat coef_1 ci_lower_1 ci_upper_1 if analysis=="all" & binary==0, mat(coef2_con) rownames(depvar)
//Binary outcomes
#delimit ;
coefplot (matrix(coef1[,1]), ms(T) mc(edkblue) lc() ciopts(lc(edkblue)) lc(edkblue) ci((coef1[,2] coef1[,3])) offset(.05)) ///
(matrix(coef2[,1]), ms(S) mc(red) ciopts(lc(red)) lstyle(p1) lc(red) ci((coef2[,2] coef2[,3])) offset(-.05)) , ///
xline(0) leg(off) ylabel(,format(%9.1f)) plotregion(lc(white)) xtitle(Risk difference) grid(none)
coeflabels(out_highbloodpressure = "Hypertension"
out_diabetes = "Diabetes"
out_stroke = "Stroke"
out_heartattack = "Heart attack"
out_depression = "Depression"
out_cancer = "Cancer"
out_dead = "Died"
out_exsmoker = "Ever smoked"
out_smoker = "Current smoker"
out_income_under_18k = "Income over £18k"
out_income_over_31k = "Income over £31k"
out_income_over_52k = "Income over £52k"
out_income_over_100k = "Income over £100k");
#delimit cr
graph export "results/figure_5.png", replace as(png)
//Continious outcomes
#delimit ;
coefplot (matrix(coef1_con[,1]), ms(T) mc(edkblue) lc() ciopts(lc(edkblue)) lc(edkblue) ci((coef1_con[,2] coef1_con[,3])) offset(.05)) ///
(matrix(coef2_con[,1]), ms(S) mc(red) ciopts(lc(red)) lstyle(p1) lc(red) ci((coef2_con[,2] coef2_con[,3])) offset(-.05)) , ///
xline(0) leg(off) ylabel(,format(%9.1f)) plotregion(lc(white)) xtitle(Risk difference) grid(none)
coeflabels(out_gripstrength = "Grip strength (kg)"
out_arterial_stiffness = "Arterial Stiffness"
out_height = "Height (cm)"
out_bmi = "BMI (kg/m2)"
out_dia_bp = "Diastolic blood pressure (mmHg)"
out_sys_bp = "Systolic blood pressure (mmHg)"
out_intell = "Intelligence (0 to 13)"
out_happiness = "Happiness (0 to 5 Likert)"
out_alcohol = "Alcohol consumption (1 low, 5 high)"
out_sedentary = "Hours watching television"
out_phys_m_act = "Moderate exercise (days/week)"
out_phys_v_act = "Vigorous exercise (days/week)");
#delimit cr
graph export "results/figure_6.png", replace as(png)
/*
*/
rm "working data/temp1.dta"
rm "working data/temp2.dta"
rm "working data/temp3.dta"
//Clean the results for Table 4a - the IV results with edu years
clean_iv 1 "results/outcome_IV2_all"
clean_iv 2 "results/outcome_IV2_male"
clean_iv 3 "results/outcome_IV2_female"
use "working data/temp1.dta",clear
append using "working data/temp2"
append using "working data/temp3"
order depvar N coef ci_lower ci_upper pval estatp cdf
save "results/table4a",replace
rm "working data/temp1.dta"
rm "working data/temp2.dta"
rm "working data/temp3.dta"
//Clean the observational association of eduyears and the outcome
use "results/outcome_act_exposure2_0",
append using "results/outcome_act_exposure2_1",
append using "results/outcome_act_exposure2_2",