-
Notifications
You must be signed in to change notification settings - Fork 0
Expand file tree
/
Copy path8_1c_Cox_regression.do
More file actions
468 lines (341 loc) · 16.4 KB
/
Copy path8_1c_Cox_regression.do
File metadata and controls
468 lines (341 loc) · 16.4 KB
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
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
******************************************
** 8_1c cox regression
**
** Cox regression
******************************************
*==== define key date of wave
global wave1 = date("15Dec2020" , "DMY" ) // when wild-type ends
global wave2 = date("15May2021" , "DMY") // when alpha ends
global wave3 = date("15Dec2021" , "DMY") // when delta ends
global lockdown1start = date("23Mar2020" , "DMY") // first lockdown starts
global lockdown1end = date("23Jun2020" , "DMY") // first lockdown ends
global lockdown2start = date("31Oct2020" , "DMY") // second lockdown starts
global lockdown2end = date("2Dec2020" , "DMY") // second lockdown ends
global cox_option "level(95)"
global cox_option2 "level(99)"
******************************************
* Planned admission
******************************************
use "$derived_data_path\8_1a2_planned_admission_birth_conceived_feb2020jul2020.dta" , clear
*==== keep the first planned admission only , look at time to the first admission
* The dataset include only the first episode of each admission
* flag the first planned admission for each child
bysort tokenid (adm_no) : gen planadm_order = _n
keep if planadm_order==1 // keep only the first admission or the record of the child (if the child does not have admission)
*==== clean birth day and date of death
gen bday_raw = bday
format bday_raw %td
capture drop tmp*
gen tmp1 = string(bday , "%td")
gen tmp2= substr(tmp1 , 1 , 2 )
gen tmp3 = mofd(admidate) // month of admission date
format tmp3 %tm
// some children died on the same day of birthday or before birthday, so have 0 follow-up time and would be excluded automatically from calculating rates and cox-regressions. For these children, change the date of death to one day after birthday, so that they will be included in the anlaysis
gen tmp4 = mofd(dod)
format tmp4 %tm
replace dod = bday+1 if dod==bday_raw & death==1
replace bday = dod-1 if death==1 & dod<bday & tmp2=="15" & tmp4==month_bday
capture drop tmp*
*==== set survival analysis variable
gen censordate = round(bday + 365.25*2.5 , 1)
format censordate %td
**---- end time point
* date of the first admission
gen finalcensor1=censordate
replace finalcensor1 = admidate if admidate!=. & admidate<=censordate
format finalcensor1 %td
* date of death
gen finalcensor2 = censordate
replace finalcensor2 = dod if dod!=. & dod<=censordate
format finalcensor2 %td
* final censor date
gen finalcensor = min(finalcensor1 , finalcensor2)
format finalcensor %td
label variable finalcensor "2.5y birthday, first admission, death, whichever the earliest"
drop finalcensor1 finalcensor2 censordate
**---- outcome
gen outcome=0
replace outcome=1 if admidate!=. & admidate<=finalcensor
**---- time to event (exposed time)
gen time_days=finalcensor-bday
**---- declaire
stset time_days , failure(outcome) scale(365.25) id(tokenid)
* recode matage_compl to small groups
recode matage_compl (min/24=1 "<25") (25/29=2 "25-29") (30/34=3 "30-34") (35/39=4 "35-39") (40/max=5 ">=40") (.=.) , gen(matage_catd)
label variable matage_catd "maternal age"
**---- save data
save "$derived_data_path/8_1c1_first_episode_first_non-birth_planned_admission.dta" , replace
use "$derived_data_path/8_1c1_first_episode_first_non-birth_planned_admission.dta" , clear
**==== KM curve
* England
sts graph , by(infect_variant) failure legend(pos(6) row(1)) xtitle("Age in years") title("Planned admission") legend(label(1 "Negative") label(2 "Wild-type") label(3 "Alpha") label(4 "Untested")) plot1opts(lcolor(green)) plot2opts(lcolor(orange)) plot3opts(lcolor(red)) plot4opts(lcolor(blue)) name(planned , replace) saving(planned_england , replace) ylabel(0 0.05 0.1 ) ytitle("Proportion having a planned admission")
graph export "$output_path\8_1c1_KM_planned_admission_2halfy_England.png" , replace
*==== Main analysis - Cox regression
**---- models
* sample size for cox regression
count if ethgroup_mother_compl!=. & imd04decile_mother_compl!=. & mcc_yorn!=. & matage_compl!=.
* crude model
stcox i.infect_variant ///
if ethgroup_mother_compl!=. & imd04decile_mother_compl!=. & mcc_yorn!=. & matage_compl!=. & londonres!=. , $cox_option2
est store plan0
* + mother ethnicity
stcox i.infect_variant ///
ib5.ethgroup_mother_compl ///
if ethgroup_mother_compl!=. & imd04decile_mother_compl!=. & mcc_yorn!=. & matage_compl!=. , $cox_option2
est store plan1
* + baby imd04 decile
stcox i.infect_variant ///
ib5.ethgroup_mother_compl ///
i.imd04decile_mother_compl ///
if ethgroup_mother_compl!=. & imd04decile_mother_compl!=. & mcc_yorn!=. & matage_compl!=. , $cox_option2
est store plan2
* + mother's chronic condition
stcox i.infect_variant ///
ib5.ethgroup_mother_compl ///
i.imd04decile_mother_compl ///
i.mcc_yorn ///
if ethgroup_mother_compl!=. & imd04decile_mother_compl!=. & mcc_yorn!=. & matage_compl!=. , $cox_option2
est store plan3
* + maternal age
stcox i.infect_variant ///
ib5.ethgroup_mother_compl ///
i.imd04decile_mother_compl ///
i.mcc_yorn ///
i.matage_catd ///
if ethgroup_mother_compl!=. & imd04decile_mother_compl!=. & mcc_yorn!=. & matage_compl!=. , $cox_option2
est store plan4
* + london residence
stcox i.infect_variant ///
ib5.ethgroup_mother_compl ///
i.imd04decile_mother_compl ///
i.mcc_yorn ///
i.matage_catd ///
i.londonres ///
if ethgroup_mother_compl!=. & imd04decile_mother_compl!=. & mcc_yorn!=. & matage_compl!=. , $cox_option2
est store plan5
* + month_conception
stcox i.infect_variant ///
ib5.ethgroup_mother_compl ///
i.imd04decile_mother_compl ///
i.mcc_yorn ///
i.matage_catd ///
i.londonres ///
i.month_conception ///
if ethgroup_mother_compl!=. & imd04decile_mother_compl!=. & mcc_yorn!=. & matage_compl!=. , $cox_option2 cformat(%9.2f)
est store plan6
**---- proportional hazard assumption
* time-dependent variant
stcox i.infect_variant ///
ib5.ethgroup_mother_compl ///
i.imd04decile_mother_compl ///
i.mcc_yorn ///
i.matage_catd ///
i.londonres ///
i.month_conception ///
if ethgroup_mother_compl!=. & imd04decile_mother_compl!=. & mcc_yorn!=. & matage_compl!=. , $cox_option2 nohr
est store planned_notime
stcox i.infect_variant ///
ib5.ethgroup_mother_compl ///
i.imd04decile_mother_compl ///
i.mcc_yorn ///
i.matage_catd ///
i.londonres ///
i.month_conception ///
if ethgroup_mother_compl!=. & imd04decile_mother_compl!=. & mcc_yorn!=. & matage_compl!=. , $cox_option2 nohr tvc(i.infect_variant) texp(ln(_t))
est store planned_time
lrtest planned_notime planned_time
*---- impact of trimester of infection on child admission
gen tmp1 = (positive_term1==1)
gen tmp2 = (positive_term2==1)
gen tmp3 = (positive_term3==1)
gen tmp4 = tmp1 + tmp2 + tmp3
gen positive_term = 1 if tmp1==1
replace positive_term = 2 if tmp2==1
replace positive_term = 3 if tmp3==1
replace positive_term = . if tmp4>1
capture drop tmp*
stcox i.infect_variant ib2.positive_term if positive_term!=. & ethgroup_mother_compl!=. & imd04decile_mother_compl!=. & mcc_yorn!=. & matage_compl!=. , $cox_option2
est store tri_plan0
stcox i.infect_variant ib2.positive_term ///
ib5.ethgroup_mother_compl ///
i.imd04decile_mother_compl ///
i.mcc_yorn ///
i.matage_catd ///
i.londonres ///
i.month_conception ///
if positive_term!=. & ethgroup_mother_compl!=. & imd04decile_mother_compl!=. & mcc_yorn!=. & matage_compl!=. , $cox_option2
est store tri_plan1
******************************************
* emergency admission
******************************************
use "$derived_data_path\8_1a2_emergency_admission_birth_conceived_feb2020jul2020.dta" , clear
*==== keep the first planned admission only , look at time to the first admission
* The dataset include only the first episode of each admission
* flag the first planned admission for each child
bysort tokenid (adm_no) : gen emergency_order = _n
keep if emergency_order==1
*==== clean birth day and date of death
gen bday_raw = bday
format bday_raw %td
// some children with no birth admission were identified in HES APC, so the date of birth is recorded as the 15th of the month of birth, but there are admissions with a admission date earlier than the 15th of the month of birth. For these children, change the birthday to one day before the first admission date
capture drop tmp*
gen tmp1 = string(bday , "%td")
gen tmp2= substr(tmp1 , 1 , 2 )
gen tmp3 = mofd(admidate) // month of admission date
format tmp3 %tm
replace bday = admidate-1 if admidate!=. & tmp3==month_bday & tmp2=="15" & admidate<bday_raw
replace bday = admidate-1 if admidate==bday_raw & month_bday==tmp3
// some children died on the same day of birthday or before birthday, so have 0 follow-up time and would be excluded automatically from calculating rates and cox-regressions. For these children, change the date of death to one day after birthday, so they will be included in the anlaysis
gen tmp4 = mofd(dod)
format tmp4 %tm
replace dod = bday+1 if dod==bday_raw & death==1
replace bday = dod-1 if death==1 & dod<bday & tmp2=="15" & tmp4==month_bday
capture drop tmp*
*==== set survival analysis variable
gen censordate = bday + 365.25*2.5
format censordate %td
**---- end time point
* date of the first admission
gen finalcensor1=censordate
replace finalcensor1 = admidate if admidate!=. & admidate<=censordate
format finalcensor1 %td
* date of death
gen finalcensor2 = censordate
replace finalcensor2 = dod if dod!=. & dod<=censordate
format finalcensor2 %td
* final censor date
gen finalcensor = min(finalcensor1 , finalcensor2)
format finalcensor %td
label variable finalcensor "2.5y birthday, first admission, death, whichever the earliest"
drop finalcensor1 finalcensor2 censordate
**---- outcome
gen outcome=0
replace outcome=1 if admidate!=. & admidate<=finalcensor
**---- time to event (exposed time)
gen time_days=finalcensor-bday
**---- declaire
stset time_days , failure(outcome) scale(365.25) id(tokenid)
** categorise maternal age to small groups
recode matage_compl (min/24=1 "<25") (25/29=2 "25-29") (30/34=3 "30-34") (35/39=4 "35-39") (40/max=5 ">=40") (.=.) , gen(matage_catd)
label variable matage_catd "maternal age"
**---- save data
save "$derived_data_path/8_1c2_first_episode_first_non-birth_emergency_admission.dta" , replace
*==== Hospital admission rate
use "$derived_data_path/8_1c2_first_episode_first_non-birth_emergency_admission.dta" , clear
**---- KM curve
* England
sts graph , by(infect_variant) failure legend(pos(6) row(1)) xtitle("Age in years") title("Emergency admission") legend(label(1 "Negative") label(2 "Wild-type") label(3 "Alpha") label(4 "Untested")) plot1opts(lcolor(green)) plot2opts(lcolor(orange)) plot3opts(lcolor(red)) plot4opts(lcolor(blue)) name(planned , replace) saving(emergency_england , replace) ylabel(0 0.1 0.2 0.3 0.4) ytitle("Proportion having an emergency admission")
graph export "$output_path\8_1c1_KM_emergency_admission_2halfy_England.png" , replace
*==== Main analysis: Cox regression
**---- models
* sample size for cox regression
count if ethgroup_mother_compl!=. & imd04decile_mother_compl!=. & mcc_yorn!=. & matage_compl!=.
* crude model
stcox i.infect_variant ///
if ethgroup_mother_compl!=. & imd04decile_mother_compl!=. & mcc_yorn!=. & matage_compl!=. , $cox_option2
est store emergency0
* + mother ethnicity
stcox i.infect_variant ///
ib5.ethgroup_mother_compl ///
if ethgroup_mother_compl!=. & imd04decile_mother_compl!=. & mcc_yorn!=. & matage_compl!=. , $cox_option2
est store emergency1
* + baby imd04 decile
stcox i.infect_variant ///
ib5.ethgroup_mother_compl ///
i.imd04decile_mother_compl ///
if ethgroup_mother_compl!=. & imd04decile_mother_compl!=. & mcc_yorn!=. & matage_compl!=. , $cox_option2
est store emergency2
* + mother's chronic condition
stcox i.infect_variant ///
ib5.ethgroup_mother_compl ///
i.imd04decile_mother_compl ///
i.mcc_yorn ///
if ethgroup_mother_compl!=. & imd04decile_mother_compl!=. & mcc_yorn!=. & matage_compl!=. , $cox_option2
est store emergency3
* + maternal age
stcox i.infect_variant ///
ib5.ethgroup_mother_compl ///
i.imd04decile_mother_compl ///
i.mcc_yorn ///
i.matage_catd ///
if ethgroup_mother_compl!=. & imd04decile_mother_compl!=. & mcc_yorn!=. & matage_compl!=. , $cox_option2
est store emergency4
* + london residence
stcox i.infect_variant ///
ib5.ethgroup_mother_compl ///
i.imd04decile_mother_compl ///
i.mcc_yorn ///
i.matage_catd ///
i.londonres ///
if ethgroup_mother_compl!=. & imd04decile_mother_compl!=. & mcc_yorn!=. & matage_compl!=. , $cox_option2
est store emergency5
* + month of conception
stcox i.infect_variant ///
ib5.ethgroup_mother_compl ///
i.imd04decile_mother_compl ///
i.mcc_yorn ///
i.matage_catd ///
i.londonres ///
i.month_conception ///
if ethgroup_mother_compl!=. & imd04decile_mother_compl!=. & mcc_yorn!=. & matage_compl!=. , $cox_option2 cformat(%9.2f)
est store emergency6
**---- propotional hazard assumption
stcox i.infect_variant ///
ib5.ethgroup_mother_compl ///
i.imd04decile_mother_compl ///
i.mcc_yorn ///
i.matage_catd ///
i.londonres ///
i.month_conception ///
if ethgroup_mother_compl!=. & imd04decile_mother_compl!=. & mcc_yorn!=. & matage_compl!=. , $cox_option2
est store emergency_notime
stcox i.infect_variant ///
ib5.ethgroup_mother_compl ///
i.imd04decile_mother_compl ///
i.mcc_yorn ///
i.matage_catd ///
i.londonres ///
i.month_conception ///
if ethgroup_mother_compl!=. & imd04decile_baby_compl!=. & mcc_yorn!=. & matage_compl!=. , $cox_option2 tvc(i.infect_variant) texp(ln(_t))
est store emergency_time
lrtest emergency_time emergency_notime
*---- impact of trimester of infection on child admission
gen tmp1 = (positive_term1==1)
gen tmp2 = (positive_term2==1)
gen tmp3 = (positive_term3==1)
gen tmp4 = tmp1 + tmp2 + tmp3
gen positive_term = 1 if tmp1==1
replace positive_term = 2 if tmp2==1
replace positive_term = 3 if tmp3==1
replace positive_term = . if tmp4>1
catpure drip tmp*
stcox i.infect_variant ib2.positive_term if positive_term!=. & ethgroup_mother_compl!=. & imd04decile_mother_compl!=. & mcc_yorn!=. & matage_compl!=. , $cox_option2
est store tri_emergency0
stcox i.infect_variant ib2.positive_term ///
ib5.ethgroup_mother_compl ///
i.imd04decile_mother_compl ///
i.mcc_yorn ///
i.matage_catd ///
i.londonres ///
i.month_conception ///
if positive_term!=. & ethgroup_mother_compl!=. & imd04decile_mother_compl!=. & mcc_yorn!=. & matage_compl!=. , $cox_option2
est store tri_emergency1
******************************************
* plots of hazard ratio
******************************************
* exposure groups
coefplot ( emergency0 , label(Crude)) ( emergency6, label(Adjusted)) , bylabel("Emergency admission") || plan0 plan6 , bylabel("Planned admission") drop(_cons) byopt(row(1) legend(pos(6))) xline(1) keep(*.infect_variant) mlabel mlabsize(*1.35) mlabposition(2) format(%4.3f) ci(99 box) xlabel(0.8 "0.8" 1 "1" 1.2 "1.2") eform rename(0.infect_variant = "Test-negative" 1.infect_variant = "Positive: Wild-type" 2.infect_variant="Positive: Alpha" 10.infect_variant = "No-recorded-result") coeflabels(, labsize(medlarge)) legend(row(1)) baselevels
graph export "$output_path\8_1c2_coxreg_planned_emergency_2halfy_99ci.png" , replace width(1600) height(900)
*---- trimester and exposure groups
* plot coefficients
label define positive_term 1 "First trimester" 2 "Second trimester" 3 "Third trimester" , replace
label value positive_term positive_term
coefplot ( tri_emergency0 , label(Crude)) ( tri_emergency1, label(Adjusted)) , bylabel("Emergency admission") || tri_plan0 tri_plan1 , bylabel("Planned admission") ///
drop(_cons) byopt(row(1) legend(pos(6))) xline(1) ///
keep(*.infect_variant *.positive_term) ///
mlabel mlabsize(*1.35) mlabposition(2) format(%4.3f) ci(99 box) ///
xlabel(0.8 "0.8" 1 "1" 1.2 "1.2") eform ///
headings(1.infect_variant = "{bf:Variant of in utero exposure}" ///
1.positive_term = "{bf: Trimester at in utero exposure}" , nogap labsize(*0.9)) ///
coeflabels(, labsize(medlarge)) legend(row(1)) baselevels
graph export "$output_path\8_1c3_coxreg_trimester_planned_emergency_2halfy.png" , replace width(1600) height(900)