-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy path7_sensitivity_analysis.ado
205 lines (187 loc) · 7.36 KB
/
7_sensitivity_analysis.ado
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
*------------------------------------------------------------------------------*
* This script performs sensitivity analysis for the degrees of freedom in the--*
* flexible parametric survival model-------------------------------------------*
*------------------------------------------------------------------------------*
*-------------------------------*
*** Install relevant packages ***
*-------------------------------*
* tell Stata to never pause
set more off, perm
* style of decimal point is a comma
set dp comma, perm
capture program drop sensitivity_analysis
program sensitivity_analysis
args baseline_min baseline_max ///
spline_min spline_max ///
TD_min TD_max ///
df_spline_year df_spline_age ///
df_baseline df_TD ///
max_CBS_age max_fupdat bool_CI bool_stage ///
strs_path results_data_path SA_path path_ado_files ///
disease age1 age2 age3 age4 start_year end_year
capture drop *
/* Estimate the loss in expectation of life */
quietly cd "`strs_path'"
/* Vary degrees of freedom of baseline hazard */
file open SA_`disease' using "`SA_path'\SA_`bool_stage'_`disease'.txt", write replace
local min_BIC = 1e+100
local chosen_df_baseline = 0
forval try_df_baseline = `baseline_min'/`baseline_max' {
/* Load data */
use "`results_data_path'\\`disease'_year`df_spline_year'_age`df_spline_age'_`bool_stage'.dta", clear
if ("`bool_stage'" == "with_stage"){
quietly keep if stage4 == 0
quietly drop stage4
}
/* Train model */
if (`df_TD' == 0 & "`bool_stage'" == "without_stage"){
capture noisily: quietly stpm2 syr* fem* sag*, scale(hazard) df(`try_df_baseline') bhazard(rate)
}
else if (`df_TD' == 0 & "`bool_stage'" == "with_stage"){
capture noisily: quietly stpm2 syr* fem* sag* stage*, scale(hazard) df(`try_df_baseline') bhazard(rate)
}
else if (`df_TD' > 0 & "`bool_stage'" == "without_stage"){
capture noisily: quietly stpm2 syr* fem* sag*, scale(hazard) df(`try_df_baseline') bhazard(rate) tvc(syr* fem* sag*) dftvc(`df_TD')
}
else if (`df_TD' > 0 & "`bool_stage'" == "with_stage"){
capture noisily: quietly stpm2 syr* fem* sag* stage*, scale(hazard) df(`try_df_baseline') bhazard(rate) tvc(syr* fem* sag* stage*) dftvc(`df_TD')
}
// continue if convergence not achieved
if (_rc+(e(converge)==0) > 0){
continue
}
else{
/* If BIC smaller, store as best model */
estimates store baseline_`try_df_baseline'
estat ic
mat s=r(S)
local new_BIC = (s[1,6])
if (`new_BIC' < `min_BIC'){
local min_BIC = `new_BIC'
local chosen_df_baseline = `try_df_baseline'
}
}
}
file write SA_`disease' "df baseline: `chosen_df_baseline'" _n
/* Vary degrees of freedom of spline of age at diagnosis */
local min_BIC = 1e+100
local chosen_df_spline_age = 0
forval try_df_spline_age = `spline_min'/`spline_max' {
/* Load data */
use "`results_data_path'\\`disease'_year`df_spline_year'_age`try_df_spline_age'_`bool_stage'.dta", clear
if ("`bool_stage'" == "with_stage"){
keep if stage4 == 0
drop stage4
}
/* Train model */
if (`df_TD' == 0 & "`bool_stage'" == "without_stage"){
capture noisily: quietly stpm2 syr* fem* sag*, scale(hazard) df(`chosen_df_baseline') bhazard(rate)
}
else if (`df_TD' == 0 & "`bool_stage'" == "with_stage"){
capture noisily: quietly stpm2 syr* fem* sag* stage*, scale(hazard) df(`chosen_df_baseline') bhazard(rate)
}
else if (`df_TD' > 0 & "`bool_stage'" == "without_stage"){
capture noisily: quietly stpm2 syr* fem* sag*, scale(hazard) df(`chosen_df_baseline') bhazard(rate) tvc(syr* fem* sag*) dftvc(`df_TD')
}
else if (`df_TD' > 0 & "`bool_stage'" == "with_stage"){
capture noisily: quietly stpm2 syr* fem* sag* stage*, scale(hazard) df(`chosen_df_baseline') bhazard(rate) tvc(syr* fem* sag* stage*) dftvc(`df_TD')
}
// continue if convergence not achieved
if (_rc+(e(converge)==0) > 0){
continue
}
else{
/* If BIC smaller, store as best model */
estimates store spline_age_`try_df_spline_age'
estat ic
mat s=r(S)
local new_BIC = (s[1,6])
if (`new_BIC' < `min_BIC'){
local min_BIC = `new_BIC'
local chosen_df_spline_age = `try_df_spline_age'
}
}
}
file write SA_`disease' "df spline age: `chosen_df_spline_age'" _n
/* Vary degrees of freedom of spline of year at diagnosis */
local min_BIC = 1e+100
local chosen_df_spline_year = 0
forval try_df_spline_year = `spline_min'/`spline_max' {
/* Load data */
use "`results_data_path'\\`disease'_year`try_df_spline_year'_age`chosen_df_spline_age'_`bool_stage'.dta", clear
if ("`bool_stage'" == "with_stage"){
keep if stage4 == 0
drop stage4
}
/* Train model */
if (`df_TD' == 0 & "`bool_stage'" == "without_stage"){
capture noisily: quietly stpm2 syr* fem* sag*, scale(hazard) df(`chosen_df_baseline') bhazard(rate)
}
else if (`df_TD' == 0 & "`bool_stage'" == "with_stage"){
capture noisily: quietly stpm2 syr* fem* sag* stage*, scale(hazard) df(`chosen_df_baseline') bhazard(rate)
}
else if (`df_TD' > 0 & "`bool_stage'" == "without_stage"){
capture noisily: quietly stpm2 syr* fem* sag*, scale(hazard) df(`chosen_df_baseline') bhazard(rate) tvc(syr* fem* sag*) dftvc(`df_TD')
}
else if (`df_TD' > 0 & "`bool_stage'" == "with_stage"){
capture noisily: quietly stpm2 syr* fem* sag* stage*, scale(hazard) df(`chosen_df_baseline') bhazard(rate) tvc(syr* fem* sag* stage*) dftvc(`df_TD')
}
// continue if convergence not achieved
if (_rc+(e(converge)==0) > 0){
continue
}
else{
/* If BIC smaller, store as best model */
estimates store spline_year_`try_df_spline_year'
estat ic
mat s=r(S)
local new_BIC = (s[1,6])
if (`new_BIC' < `min_BIC'){
local min_BIC = `new_BIC'
local chosen_df_spline_year = `try_df_spline_year'
}
}
}
file write SA_`disease' "df spline year: `chosen_df_spline_year'" _n
/* Vary degrees of freedom to model time-dependent effects */
local min_BIC = 1e+100
local chosen_df_TD = 0
forval try_df_TD = `TD_min'/`TD_max' {
/* Load data */
use "`results_data_path'\\`disease'_year`chosen_df_spline_year'_age`chosen_df_spline_age'_`bool_stage'.dta", clear
if ("`bool_stage'" == "with_stage"){
keep if stage4 == 0
drop stage4
}
/* Train model */
if (`try_df_TD' == 0 & "`bool_stage'" == "without_stage"){
capture noisily: quietly stpm2 syr* fem* sag*, scale(hazard) df(`chosen_df_baseline') bhazard(rate)
}
else if (`try_df_TD' == 0 & "`bool_stage'" == "with_stage"){
capture noisily: quietly stpm2 syr* fem* sag* stage*, scale(hazard) df(`chosen_df_baseline') bhazard(rate)
}
else if (`try_df_TD' > 0 & "`bool_stage'" == "without_stage"){
capture noisily: quietly stpm2 syr* fem* sag*, scale(hazard) df(`chosen_df_baseline') bhazard(rate) tvc(syr* fem* sag*) dftvc(`try_df_TD')
}
else if (`try_df_TD' > 0 & "`bool_stage'" == "with_stage"){
capture noisily: quietly stpm2 syr* fem* sag* stage*, scale(hazard) df(`chosen_df_baseline') bhazard(rate) tvc(syr* fem* sag* stage*) dftvc(`try_df_TD')
}
if (_rc+(e(converge)==0) > 0){
// continue if convergence not achieved
continue
}
else{
/* If BIC smaller, store as best model */
estimates store df_TD_`try_df_TD'
estat ic
mat s=r(S)
local new_BIC = (s[1,6])
if (`new_BIC' < `min_BIC'){
local min_BIC = `new_BIC'
local chosen_df_TD = `try_df_TD'
}
}
}
file write SA_`disease' "df TD: `chosen_df_TD'" _n
file close SA_`disease'
end