-
Notifications
You must be signed in to change notification settings - Fork 0
/
strbee.ado
1254 lines (1146 loc) · 45.8 KB
/
strbee.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
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
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
*! version 1.9.0 Ian White 7jan2022
/*************************** NOTES ********************************
v1.9.0 7jan2022 - ON SSC
small changes to km option
never show counterfactual graphs when they match observed
default is observed and untreated
change "if observed" to "as observed"
v1.8.9 13feb2020
added "if" to legends in kmgraph option
version 1.8.8 30may2018 - ON MY UCL WEBSITE
gen() now outputs correct recensoring variable
was missing
internal calculation was already correct
test now defaults to cox if adjvars is specified (but not ipe)
version 1.8.7 12feb2018 - ON SSC
help file updated to UCL
undocumented changes:
allow kmgraph(show(#...)) to control which curves are graphed
allow hr(x x) to control which curves are compared
where x=fully|observed|untreated or abbreviation
version 1.8.6 17jan2018 - ON MY UCL WEBSITE
kmgraph() and zgraph() fixed to allow quotes
version 1.8.5 27mar2017 - ON SSC
no changes to code; small changes to help file
version 1.8.4 22apr2016 - ON MY UCL WEBSITE
kmgraph greatly expanded to improve patterns and colours
version 1.8.3 31mar2016 - ON MY UCL WEBSITE
kmgraph(showall) shows all 6 graphs
kmgraph(untreated) shows just the untreated counterfactual graphs
fixed bug in hr, kmgraph and gen() options with ton()/toff() syntax
experimental psimult() option
corrected correction made in v1.8.1
version 1.8.2 24mar2016
near-zero values of psi rounded to zero
solves the "ITT P-value not found" problem which occurred e.g. with options psimin(-1) psimax(0) psistep(0.01)
version 1.8.1 22mar2016 (on beta website)
corrected error in ipe+weibull (used log HR not log TR)
version 1.8 6mar2014 (on beta website 16jun2015)
new syntax using ton and/or toff
suboptions to improve the recensoring such as ton(var, max0(1) min1(0.1))
- these replace mincens
more tests - all options for either -sts test- or -streg-
version 1.7.2 17dec2013
minor changes to help file and documentation
version 1.7.1 1oct2013
spots if current directory is not writeable & gives helpful error message
also split parsing from data checking section
version 1.7 20-21sep2012
Program crashed if no events after recensoring: now if no events it directly sets Z missing.
Behaviour after Z=. in stages 1 and 2: before it shifted psi and entered an endless loop; now it moves to next stage
New mincens option assumes that switching cannot occur before minimum observed switching times when computing the recensoring time as the minimum over all possible treatment history - potentially yields less recensoring
Hr option:
now uses strata() option if specified
if counterfactual is very close to observed, they are set exactly equal - this makes unadjusted and adjusted HR agree exactly
Trace & debug output improved
Renumbered stages
version 1.6 21may2012
new strata() option
adjvars and strata in output
corrected bug that stopped test(cox) and test(other_sts_option) working
version 1.5.5 21may2012
version 1.5.4 7feb2012
corrected level() option to allow level(99.9) etc.
version 1.5.3 31jan2012
warns if xo0, xo1 options appear to have timevar eventvar instead of eventvar timevar
version 1.5.2 27jan2012
corrected naming of gen() variables
improved their description in the help file
version 1.5.1 31oct2011
minor tidying up; help file up to date.
version 1.5 14oct2011
savedta defaults to _strbee_savedta
new replay ability (except after IPE)
new zgraph option (graph still works)
zgraph now uses up-to-date graphics; graph options updated; symbol replaced with msymbol
kmgraph() and zgraph() now take (graph) options
results only tabulated in replay stage
hr works at replay stage - risky as data mustn't have changed!
program separated into estimate (creates the "savedta" results file) and output (returns r())
version 1.4.1 30sep2011
new ipestore option stores last regression under IPE
version 1.4 26aug2011
trace option no longer saves data (can do the same with gen)
data not changed in program - preserve only used when psi-Z results are loaded
new hr option
output made neater
better behaviour with increasing Z vs psi graph - suggests back-to-front treatment
version 1.3.1 25aug2011
upgrade to Stata 10
psimin and psimax are pushed out to +/-10 if necessary
preserve moved to top (handles program crash)
new psistart() option with ipe
general clean-up
makeu subroutine used - now gen works
corrected error with XO times > T
corrected small bug: program could change type of _t
version 1.2.5 24aug2011
added warning that ipe se & CI are too small
version 1.2.4 19may2011
moved preserve from top so that gen() works - but not fully tested!
version 1.2.3 24mar2011
allow adjustment for covariates using adjvars(varlist)
postfile label is a tempname
tidied up output a little
version 1.2.2 6aug2010
correct sign of wilcoxon test
version 1.2.1 11dec2002
correct label of Z in output data
version 1.2 9 May 2002
upgrade to version 7
gen() option
maxiter(#) option for ipe
correct results when Z=. at bounds
version 1.1 12feb2002
trace option outputs recensored events as well as subjects
allow savedta(x,replace) even if x doesn't exist
new option noci (useful with trace for rbeetrace.dta)
store results in char() not globals
version 1.0 5 December 2000
Certification script: H:\trtchg\strbee\update2011\strbee_cscript.do
******************************************************************/
********************* WRAPPER PROGRAM **********************
* options without arguments get "run" as argument
program define strbee,
version 10
syntax [anything] [if] [in] [using], [ZGRaph GRaph KMgraph debug HR *]
foreach opt in zgraph graph kmgraph hr {
if `"``opt''"' != "" local options `opt'(run) `options'
}
local options `options' `debug'
if "`anything'"=="" local cmd Output `anything' `if' `in' `using', `options'
else local cmd Gestimate `anything' `if' `in', `options'
if "`debug'"=="debug" {
if "`using'"=="" di as text "Estimation and output syntax"
else di as text "Output syntax"
di as input `"`cmd'"'
}
`cmd'
end
********************* START OF G-ESTIMATE PROGRAM **********************
program define Gestimate
********************* PARSING **********************
syntax varlist(max=1) [if] [in], [ ///
xo0(string) xo1(string) ton(string) toff(string) /// model options
ENDstudy(varname) psimult(string) /// alternative to xo*
TEst(string) ADJvars(varlist) STrata(varlist) /// estimator options
ipe ipecens /// estimator options
psimin(real -10) psimax(real 10) psistep(real 0) /// search options
tol(int 3) noci level(cilevel) /// search options
trace psistart(real 0) maxiter(int 100) /// search options
SAVEdta(string) IPESTore(string) /// storage options
debug noOUTput /// undocumented options
*]
local zcrit = invnorm(1-(100-`level')/200)
if "`strata'"!="" local strataopt strata(`strata')
local treat `varlist'
local Gestimateoptions `options'
* SORT OUT SAVE FILE
if "`savedta'"~="" {
tokenize "`savedta'", parse(",")
if "`3'"=="append" local append append
else if "`3'"=="replace" local replace replace
else if "`3'"~="" {
di as error "savedta(`savedta') not allowed"
exit 498
}
local savedta `1'
tokenize "`savedta'", parse(".")
if "`3'"=="" local savedta `1'.dta
if "`append'"=="append" confirm file `savedta'
else if "`replace'"=="replace" {
cap confirm file `savedta'
if _rc>0 {
di as error "Warning: file `savedta' not found"
local replace
}
}
}
else {
local savedta _strbee_savedta
cap confirm file `savedta'
if _rc==0 local replace replace
}
local new = cond(mi("`append'`replace'"),"new","")
cap noi confirm `new' file `savedta'
if _rc==603 di as error " (you may not have write-access to your current directory `c(pwd)')"
if _rc exit _rc
* PARSE TEST
if "`test'"=="" {
if !mi("`ipe'") local assumedtest weibull
else if "`adjvars'"!="" local assumedtest cox
else local assumedtest logrank
di as text "test() not specified: assumed test(" as result "`assumedtest'" as text ")"
local test `assumedtest'
}
if "`test'"=="cox" local test stcox
parsetest, `test'
local test = r(test)
local testtype = r(testtype)
if "`ipe'"=="ipe" & !inlist("`test'","exponential","weibull") {
di as error "Only test(exponential|weibull) is allowed with IPE option"
exit 498
}
if "`testtype'"=="sts_test" & !mi("`adjvars'") {
di as error "adjvars(`adjvars') ignored with test(`test')"
local adjvars
}
* STORE METHODS
if "`ipe'"=="ipe" {
local method Branson's IPE
local search Branson's IPE
}
else {
local method Robins-Tsiatis method
if `psistep'==0 local search interval bisection
else local search grid search
local ipe noipe
}
if "`ipecens'"=="" local ipecens noipecens
if "`trace'"=="" local trace notrace
* CHECK VARIABLES EXIST
confirm var `treat' `endstudy' `xo0' `xo1'
* CHECK OTHER OPTIONS
if "`gen'"~="" confirm new var `gen' d`gen'
if `psistep' != 0 & `psistep' <= epsfloat() di as error "Warning: psistep is dangerously small"
if !mi("`psimult'") {
cap assert `psimult'>=0
if _rc di as error "Error in psimult(`psimult') option: must be non-negative"
if _rc exit 498
cap assert !mi("`psimult'")
if _rc di as error "Error in psimult(`psimult') option: must be non-missing"
if _rc exit 498
local psimultopt psimult(`psimult')
}
* WHICH SYNTAX - xo0() xo1() OR ton() toff()
if !mi("`ton'`toff'") & !mi("`xo0'`xo1'") {
di as error "Error: can't combine xo0() xo1() with ton() toff()"
exit 198
}
local tsyntax = !mi("`ton'`toff'") // 1 for ton toff, 0 for xo0 xo1
********************* CHECK DATA **********************
* CHECK ST DATA
st_is 2 analysis
cap assert _t0==0
if _rc>0 {
di as error "Sorry, data must be stset with 1 record per subject and _t0=0"
exit 498
}
* CHECK CENSORING
cap assert _d==1
if _rc>0 & "`endstudy'"=="" noi di as error "Warning: you have censored data - specify endstudy() to avoid bias"
************ END OF CHECKING SYNTAX AND DATA. START OF ANALYSIS ******************
preserve
* SET UP POSTFILE: TEMPORARILY STORE RESULTS IN `results'
tempfile results
tempname postname
if "`ipe'"=="ipe" local ipepostvars psise stage
qui postfile `postname' psi Z `ipepostvars' using `results', replace double
qui {
* SELECT RECORDS
tempvar touse
mark `touse' `if' `in'
markout `touse' `varlist'
keep if `touse' & _st
drop `touse'
* CHECK TREATMENT VARIABLE
summ `treat'
local error = r(min)~=0|r(max)~=1
count if `treat'~=0 & `treat'~=1
if r(N)>0 | `error'==1 {
di as error "`treat' is not a 0/1 variable"
exit 498
}
* INITIALISE COUNTERFACTUALS
tempvar t dt u du
local type : type _t
gen `type' `t' = _t /* this is needed as `u' gets stset later */
gen byte `dt' = _d
gen double `u' = .
gen byte `du' = .
label var `t' "t"
label var `dt' "dt"
label var `u' "u, created by strbee"
label var `du' "du, created by strbee"
if "`endstudy'"!="" {
tempvar recens
gen double `recens'=.
label var `recens' "recensoring time, created by strbee"
}
else local endstudy .
tempvar isrecens
gen byte `isrecens'=.
label var `isrecens' "whether recensored, created by strbee"
* COUNT
forvalues i = 0/1 {
count if `treat'==`i'
local n`i' = r(N)
count if `treat'==`i' & _d
local d`i' = r(N)
}
* IF XO SYNTAX: CHECK SWITCH VARIABLES AND CONVERT TO TON/TOFF
tempvar toffvar tonvar
if `tsyntax'==0 { // xo0 xo1 syntax
* CONVERT XO'S TO TON, TOFF
maketonoff, t(`t') treat(`treat') toff(`toffvar') ton(`tonvar') xo0(`xo0') xo1(`xo1')
}
else { // ton toff syntax
local 0 `ton'
syntax [anything], [max0(string) max1(string) min0(string) min1(string)]
if !mi("`anything'") local tonexp `anything'
foreach opt in max0 max1 min0 min1 {
local ton`opt' ``opt''
}
local 0 `toff'
syntax [anything], [max0(string) max1(string) min0(string) min1(string)]
if !mi("`anything'") local toffexp `anything'
foreach opt in max0 max1 min0 min1 {
local toff`opt' ``opt''
}
if mi("`tonexp'") {
if mi("`toffexp'") {
di as error "Must specify tonexp or toffexp"
exit 198
}
gen double `toffvar' = `toffexp'
gen double `tonvar' = _t - `toffvar'
}
else {
gen double `tonvar' = `tonexp'
if !mi("`toffexp'") {
cap assert abs(((`toffexp') + (`tonexp') - _t ) / _t) < epsfloat()
if _rc di as error "`toffexp' and `tonexp' do not sum to _t"
gen double `toffvar' = `toffexp'
}
else gen double `toffvar' = _t - `tonvar'
}
}
* SUMMARISE SWITCHES
count if `treat'==0 & `tonvar'>0
local nxo0 = r(N)
count if `treat'==1 & `toffvar'>0
local nxo1 = r(N)
* EFFICIENT RECENSORING OPTIONS: CASE OF NO SWITCH IN ONE ARM
if `tsyntax'==0 {
if mi("`xo0'") local tonmax0 0 // no xo in arm 0
if mi("`xo1'") local toffmax1 0 // no xo in arm 1
}
else {
cap assert `tonvar'==0 if `treat'==0
if !_rc local tonmax0 0 // no xo in arm 0
cap assert `toffvar'==0 if `treat'==1
if !_rc local toffmax1 0 // no xo in arm 1
}
* EFFICIENT RECENSORING OPTIONS: USER-SPECIFIED
foreach opt in toffmin0 toffmin1 tonmin0 tonmin1 {
if !mi("``opt''") local recensopts `recensopts' `opt'(``opt'')
}
foreach opt in toffmax0 toffmax1 tonmax0 tonmax1 {
if !mi("``opt''") local recensopts `recensopts' `opt'(``opt'')
}
* `recensopts' is passed to makeu
* SET UP LOOP
if (`psimin' >= `psimax') {
di as error "psimin must be < psimax"
exit 498
}
if `psistep' > 0 local stage 4 /* 4=grid-searching */
else local stage 1 /* 1=starting/looking for lower bound,
2=looking for upper bound,
3=found bounds and searching for solutions,
9=stop */
if "`ipe'"=="ipe" local psi `psistart' // IPE
else if `psistep' == 0 local psi 0 // Bisection
else local psi `psimin' // Grid search
local estlo = .
local esthi = .
local lowlo = .
local lowhi = .
local upplo = .
local upphi = .
local iter 0
local oldstage .
noi di as text "Searching by `search'" _c
if "`trace'"=="trace" & "`endstudy'"!="." { // PRINT TOTAL NUMBERS
noi di _new as text "Total in arm 0/1: " _col(32) as result `n0' as text "/" as result `n1' as text " subjects" _c
noi di as text " (" as result `d0' as text "/" as result `d1' as text " events)"
foreach var in n0 n1 d0 d1 {
local l`var' = length("``var''")
}
}
else if "`trace'"=="trace" | "`debug'"=="debug" noi di
****** START LOOP ******
while `stage' ~= 9 {
local ++iter
if "`trace'"=="trace" | "`debug'"=="debug" {
if `stage'!=`oldstage' & "`ipe'"=="" noi di as text "Starting stage `stage'"
noi di as text "psi = " as result %10.7g `psi' " " _c
}
local oldstage `stage'
* GENERATE U, DU
noi makeu, u(`u') du(`du') t(`t') dt(`dt') treat(`treat') psi(`psi') `psimultopt' ///
toff(`toffvar') ton(`tonvar') `recensopts' `ipe' `ipecens' ///
`trace' endstudy(`endstudy') isrecens(`isrecens') `debug'
if "`trace'"=="trace" & "`endstudy'"!="." { // PRINT NUMBERS RECENSORED
count if `isrecens' & `treat'==0
noi di as text " recensoring " as result %`ln0'.0f r(N) _c
count if `isrecens' & `treat'==1
noi di as text "/" as result %`ln1'.0f `r(N)' _c
count if `dt'==1 & `isrecens' & `treat'==0
noi di as text " subjects (" as result %`ld0'.0f `r(N)' _c
count if `dt'==1 & `isrecens' & `treat'==1
noi di as text "/" as result %`ld1'.0f `r(N)' as text " events) " _c
}
* TEST
count if `du'
if r(N) {
stset `u' `du'
if "`testtype'"=="stcox" {
stcox `treat' `adjvars', `strataopt'
local psihat = _b[`treat']
local Z = _b[`treat']/_se[`treat']
}
else if "`testtype'"=="streg" {
if inlist("`test'","exponential","weibull") local time time
streg `treat' `adjvars', dist(`test') `time' `strataopt'
local psihat = -_b[`treat']
if inlist("`test'","gompertz") local psihat = -`psihat' // gompertz is the only PH model
local Z = `psihat'/_se[`treat']
}
else if "`testtype'"=="sts_test" {
sts test `treat', mat(U V) `test' `strataopt'
local Z = -U[1,1]/sqrt(V[1,1])
if "`test'"=="logrank" local Z = -`Z' // logrank sign is wrong
}
}
else local Z = .
if "`trace'"=="trace" | "`debug'"=="debug" {
if `Z'==. noi di as error "Z = missing " _c
else noi di as text "Z = " as result %10.7f `Z' " " _c
if inlist("`test'","exponential","weibull","cox") {
noi di _skip(3) as text "psi-hat= " as result %10.7f `psihat' _c
}
}
else if `Z'==. noi di as error "Warning: Z=. at psi=`psi' " _c
else noi di as text "." _continue
* Post results, and finish IPE
if "`ipe'"=="ipe" {
if "`trace'"=="trace" | "`debug'"=="debug" noi di
if abs(`psihat'-`psi') <= 10^(-`tol') {
noi di as text "Convergence achieved with step < 10^-`tol'"
local stage 9
}
local psi = `psihat'
local psise = _se[`treat']
if `iter'>=`maxiter' {
di as error "`maxiter' iterations exceeded"
local stage 9
}
post `postname' (`psi') (`Z') (_se[`treat']) (`stage')
continue
}
else post `postname' (`psi') (`Z')
* update results
if `Z'<0 local esthi = min(`esthi',`psi')
else local estlo = max(`estlo',`psi')
if `Z'<-`zcrit' local upphi = min(`upphi',`psi')
else local upplo = max(`upplo',`psi')
if `Z'<`zcrit' local lowhi = min(`lowhi',`psi')
else local lowlo = max(`lowlo',`psi')
if "`debug'"=="debug" noi di as text ///
" esthi=" `esthi' " estlo=`estlo' upphi=`upphi' upplo=`upplo'" ///
" lowhi=`lowhi' lowlo=`lowlo' " _c
* CHOOSE NEXT PSI
if `psistep' == 0 {
* interval bisection
local estOK = (`esthi'-`estlo')<10^(-`tol') | `esthi'==. | `estlo'==.
local hiOK = (`upphi'-`upplo')<10^(-`tol') | `upphi'==. | `upplo'==.
local lowOK = (`lowhi'-`lowlo')<10^(-`tol') | `lowhi'==. | `lowlo'==.
if `stage' == 1 { // looking for psimin
if `lowlo'==. & `psi'>`psimin' & `Z'!=. {
if "`debug'"=="debug" noi di as text "Extending psi range to left " _c
local psi = `psi'-1
}
else if `estlo'==. {
noi di as error "Failed to find a psi for which Z>=0 " _c
noi di as error "(this is often because there is more treatment in the control arm) " _c
local estOK 0
local stage 9
}
else {
local psi 1
local stage 2
}
if `Z'==. & ("`debug'"=="debug") di as error "Terminating stage 1 " _c
}
else if `stage' == 2 { // looking for psimax
if `upphi'==. & `psi'<=`psimax'-1 & `Z'!=. {
if "`debug'"=="debug" noi di as text "Extending psi range to right " _c
local psi = `psi'+1
}
else {
local psi = (`esthi'+`estlo')/2
local stage 3
}
if `Z'==. & ("`debug'"=="debug") di as error "Terminating stage 2" _c
}
else if ~`estOK' {
if "`debug'"=="debug" noi di as text "Improving estimate " _c
local psi = (`esthi'+`estlo')/2
}
else if ~`hiOK' & "`ci'"~="noci" {
if "`debug'"=="debug" noi di as text "Improving upper CL " _c
local psi = (`upphi'+`upplo')/2
}
else if ~`lowOK' & "`ci'"~="noci" {
if "`debug'"=="debug" noi di as text "Improving lower CL " _c
local psi = (`lowhi'+`lowlo')/2
}
else {
if "`debug'"=="debug" noi di as text _new "FINISH!!! " _c
local stage 9
}
} /* END OF INTERVAL BISECTION LOOP */
else {
* grid search
local psi = `psi' + `psistep'
if abs(`psi')<epsfloat() local psi 0 // 24mar2016: avoids occasional "ITT P-value not found"
if `psi'>`psimax'+epsfloat() local stage 9
local estOK = !mi(`estlo',`esthi')
}
if "`trace'"=="trace" | "`debug'"=="debug" | `Z'==. noi di
}
postclose `postname'
stset `t' `dt'
if "`ipe'"=="ipe" {
local psilow = `psi' - `zcrit'*`psise'
local psiupp = `psi' + `zcrit'*`psise'
if "`ipestore'"!="" est store `ipestore'
}
} /* end of qui */
*** SORT OUT SAVE-DATA FILE
if "`append'"=="append" {
* check method chars match
qui use `savedta', clear
if "`_dta[type]'" ~= "strbee" {
di as error "Data `savedta' wasn't saved by strbee"
exit 498
}
if "`_dta[method]'" ~= "`method'" {
di as error "Can't append: `savedta' used method `_dta[method]'"
exit 498
}
if "`_dta[test]'" ~= "`test'" {
di as error "Can't append: `savedta' used test `_dta[test]'"
exit 498
}
if "`_dta[treat]'" ~= "`treat'" {
di as error "Can't append: `savedta' used treatment=`_dta[treat]'"
exit 498
}
if "`_dta[xo0]'" ~= "`xo0'" {
di as error "Can't append: `savedta' used xo0=`_dta[xo0]'"
exit 498
}
if "`_dta[xo1]'" ~= "`xo1'" {
di as error "Can't append: `savedta' used xo1=`_dta[xo10]'"
exit 498
}
if "`_dta[endstudy]'" ~= "`endstudy'" {
di as error "Can't append: `savedta' used endstudy=`_dta[endstudy]'"
exit 498
}
noi di as text "Appending results from `savedta'"
qui use `results', clear
append using `savedta'
local estOK .
* check for duplicates
sort psi
di
tempvar dup dupok
gen byte `dup' = psi==psi[_n-1]
gen byte `dupok' = `dup' & Z==Z[_n-1]
qui count if `dupok'
if r(N)>0 noi di as error r(N) " duplicate records deleted"
qui drop if `dupok'
qui count if `dup' & ~`dupok'
if r(N)>0 noi di as error "Save file has duplicate psi with different Z"
drop `dup' `dupok'
}
else use `results', clear
* copy estimation summaries from locals to characteristics
local charlist method test treat xo0 xo1 endstudy estOK n0 n1 d0 d1 nxo0 nxo1 ///
adjvars strata recensopts tsyntax toffexp tonexp psimult
foreach char in charlist `charlist' {
char _dta[`char'] ``char''
}
char _dta[type] strbee
di
noi save `savedta', replace
di
restore
********************* NOW CALL THE OUTPUT PROGRAM **********************
if "`output'"!="nooutput" {
local cmd Output using `savedta', level(`level') `debug' `Gestimateoptions'
if "`debug'"=="debug" di as input `"strbee.`cmd'"'
`cmd'
}
end
********************* END OF G-ESTIMATE PROGRAM **********************
********************* START OF OUTPUT PROGRAM **********************
program define Output, rclass
syntax [using/], [list ZGRaph(string asis) psimin(real -10) psimax(real 10) /// output options
hr(string) KMgraph(string asis) gen(string) /// output options
GRaph(string asis) /// undocumented: equivalent to zgraph(string)
level(cilevel) debug /// undocumented options
]
if "`using'"=="" {
local using `r(savedta)'
if "`using'"=="" {
di as error "Please specify the 'using' file - no default found"
exit 498
}
}
if `"`graph'"'!="" & `"`zgraph'"'!="" {
di as error "Can't have both graph and zgraph"
exit 198
}
if `"`graph'"'!="" local zgraph `"`graph'"'
local zcrit = invnorm(1-(100-`level')/200)
local zcri = round(`zcrit',.01)
******************** LOAD RESULTS (Z vs PSI) ********************
preserve
if "`debug'"=="debug" noi di as text "Loading results file `using'"
use `using', clear
* copy estimation summaries from characteristics to locals
foreach char in `_dta[charlist]' {
local `char' `_dta[`char']'
}
if "`method'"=="Branson's IPE" local ipe ipe
else if "`method'"!="Robins-Tsiatis method" exit 497
if !mi("`psimult'") {
local psimultopt psimult(`psimult')
}
******************** HEADER FOR RESULTS ********************
* SET UP OUTPUT
local colb as result _col(24)
local col2 _col(17)
local col3 _col(31)
local col4 _col(42)
local col5 _col(54)
local line as text _dup(64) "{c -}"
* SUMMARY OF METHODS USED
di `line'
di as text "Estimating psi from accelerated life model"
di `line'
forvalues i = 0/1 {
noi di as text "`treat'=`i' " `colb' `n`i'' as text " subjects, " _c
noi di as text `colb' `d`i'' as text " events, " _c
noi di as result `nxo`i'' as text " switches"
}
di as text "Method " `colb' "`method'"
di as text "Randomised group " `colb' "`treat'"
di as text "Test " `colb' "`test'"
di as text " adjusted for " `colb' cond("`adjvars'"=="","-","`adjvars'")
di as text " stratified by " `colb' cond("`strata'"=="","-","`strata'")
if `tsyntax'==0 {
di as text "Switch in arm 0 " _c
if "`xo0'"~="" di `colb' "`xo0'"
else di `colb' "-"
di as text "Switch in arm 1 " _c
if "`xo1'"~="" di `colb' "`xo1'"
else di `colb' "-"
}
else {
if !mi("`tonexp'") di as text "Time on treatment" `colb' as result "`tonexp'"
if !mi("`toffexp'") di as text "Time off treatment" `colb' as result "`toffexp'"
}
if !mi("`psimult'") di as text "Psi multiplier" `colb' "`psimult'"
di as text "Treatment variable " `colb' "`treat'"
di as text "End of study variable " _c
if "`endstudy'"~="." di `colb' "`endstudy'"
else di `colb' "none"
di as text "Recensoring options" _c
if !mi("`recensopts'") di `colb' "`recensopts'"
else di `colb' "-"
di as text "Results file " `colb' "`using'"
******************** COMPUTE RESULTS (EST, LOW, UPP, ZITT, PITT) ********************
local ZITT .
local PITT .
if "`ipe'"=="" {
if "`list'"=="list" l psi Z, noobs
sort psi
qui results 0
local est = r(est)
local estlo = r(lo)
local esthi = r(hi)
qui results `zcrit'
local low = r(est)
local lowlo = r(lo)
local lowhi = r(hi)
qui results -`zcrit'
local upp = r(est)
local upplo = r(lo)
local upphi = r(hi)
qui summ Z if psi==0
if r(N)==0 di as error "ITT P-value not found"
else if r(N)>2 & r(sd)>0 di as error "Error: multiple ITT P-values found"
else {
local ZITT = r(mean)
local PITT = chi2tail(1,(`ZITT')^2)
}
if `esthi'==. noi di as error "No point estimate: Z(psi) > 0 at all points tested"
else if `estlo'==. noi di as error "No point estimate: Z(psi) < 0 at all points tested"
if `estlo'==. | `esthi'==. {
di as error "Check data / model or try a wider psi range"
exit 498
}
if `esthi'<`estlo' noi di as error "Z(psi) = 0 has multiple solutions"
if `lowhi'==. noi di as error "No lower CL: Z(psi) > `zcri' at all points tested"
else if `lowlo'==. noi di as error "No lower CL: Z(psi) < `zcri' at all points tested"
else if `lowhi'<`lowlo' noi di as error "Z(psi) = `zcri' has multiple solutions"
if `upphi'==. noi di as error "No upper CL: Z(psi) > -`zcri' at all points tested"
else if `upplo'==. noi di as error "No upper CL: Z(psi) < -`zcri' at all points tested"
else if `upphi'<`upplo' noi di as error "Z(psi) = -`zcri' has multiple solutions"
}
else {
qui summ psi if stage==9
assert r(N)==1
local est = r(mean)
qui summ psise if stage==9
assert r(N)==1
local estse = r(mean)
local low = `est'-`zcrit'*`estse'
local upp = `est'+`zcrit'*`estse'
}
*********************** RETURN RESULTS ********************
* must precede Z-graph and KM-graph which use -syntax-
return scalar psi = `est'
return scalar psi_low = `low'
return scalar psi_upp = `upp'
return scalar Z_ITT = `ZITT'
return scalar P_ITT = `PITT'
return local savedta `using'
return local method `method'
return local adjvars `adjvars'
return local strata `strata'
forvalues i=0/1 {
return scalar n`i' = `n`i''
return scalar d`i' = `d`i''
return scalar nxo`i' = `nxo`i''
}
*********************** Z-GRAPH RESULTS ********************
if "`ipe'"=="" & `"`zgraph'"'!=`""' {
local 0 , `zgraph'
syntax, [run TItle(passthru) YTItle(string) XLAbel(passthru) YLAbel(passthru) ///
Connect(passthru) MSYMbol(passthru) XLIne(passthru) YLIne(passthru) *]
sort psi
label var Z "`test' test statistic"
if `"`title'"'==`""' local title title(strbee results)
if `"`ytitle'"'==`""' local ytitle `test' Z
if `"`ylabel'"'==`""' local ylabel ylabel(-`zcri' 0 `zcri')
if `"`connect'"'==`""' local connect connect(l)
if `"`msymbol'"'==`""' local msymbol msymbol(p)
if `"`xline'"'==`""' local xline xline(0, lcol(black))
if `"`yline'"'==`""' local yline yline(-`zcrit' 0 `zcrit', lcol(black))
scatter Z psi if psi>=`psimin' & psi<=`psimax', `connect' `msymbol' ///
`xline' `yline' `xlabel' `ylabel' `title' ytitle(`ytitle') `options'
qui count if psi<`psimin' & !mi(Z)
if r(N) di as text "psimin(`psimin') option: " as result r(N) as text " points with psi < `psimin' have not been graphed"
qui count if psi>`psimax' & !mi(Z)
if r(N) di as text "psimax(`psimax') option: " as result r(N) as text " points with psi > `psimax' have not been graphed"
}
*********************** DISPLAY RESULTS ********************
if "`est'" ~= "" {
di as text _new "Results for log acceleration ratio psi"
local fmt as result %9.0g
di `line'
di as text `col2' "psi" `col3' "P-value" `col4' "[`level'% Conf. interval]"
di `line'
di as text "Best" `col2' `fmt' `est' `col3' %5.3f `PITT' `col4' `fmt' `low' `col5' `fmt' `upp'
if "`ipe'"!="ipe" {
di as text "lower bound" `col2' `fmt' `estlo' `col4' `fmt' `lowlo' `col5' `fmt' `upplo'
di as text "upper bound" `col2' `fmt' `esthi' `col4' `fmt' `lowhi' `col5' `fmt' `upphi'
}
di `line'
if "`ipe'"=="ipe" {
noi di as text "Warning: standard error and confidence interval do not allow for switch"
noi di as text "Use bootstrap to get better standard error and confidence interval"
}
}
********************* END OF Z-PSI ANALYSIS: RESTORE DATA ********************
restore
**************** POST-ESTIMATION OPTIONS: GEN, HR, KMGRAPH ***************
*** GENERATE U
if "`estOK'"!="0" {
if ("`gen'"!="" | "`hr'"!="" | `"`kmgraph'"'!=`""') {
tempvar u du
qui gen double `u'=.
qui gen byte `du'=.
if "`endstudy'"!="." {
tempvar cu
qui gen double `cu'=.
label var `cu' "Recensoring time for `gen'"
local ct `endstudy'
}
local psiround=string(`est',"%9.0g")
label var `u' "Counterfactual untreated outcome for psi=`psiround'"
if "`ipe'"=="ipe" label var `u' "Counterfactual ideal outcome for psi=`psiround'"
label var `du' "Event indicator for `gen'"
* recreate data: corrected for v1.8.3
tempvar toffvar tonvar
if `tsyntax'==0 { // xo0 xo1 syntax
maketonoff, t(_t) treat(`treat') toff(`toffvar') ton(`tonvar') xo0(`xo0') xo1(`xo1')
}
else if !mi("`tonexp'") {
gen `tonvar' = `tonexp'
gen `toffvar' = _t - `tonexp'
}
else if !mi("`toffexp'") {
gen `tonvar' = _t - `toffexp'
gen `toffvar' = `toffexp'
}
makeu, u(`u') du(`du') t(_t) dt(_d) treat(`treat') psi(`est') `psimultopt' ///
toff(`toffvar') ton(`tonvar') `recensopts' recens(`cu') ///
`ipe' `ipecens' endstudy(`endstudy') `debug'
}
*** HR AND KM GRAPH
if "`hr'" != "" | `"`kmgraph'"'!=`""' {
* COMMON CODE
preserve
if !mi("`psimult'") {
tempvar psimultvar
gen double `psimultvar' = `psimult'
local timespsimult *`psimultvar'
}
if "`ipe'"=="ipe" qui replace `u' = exp(`est'`timespsimult')*`u' if `treat'==1
keep `u' `du' _t _d `treat' `adjvars' `strata' `psimultvar'
tempvar u1 treat0 type
rename `treat' `treat0'
gen double `u1' = `u'*exp(-`est'`timespsimult')
* correct rounding errors
qui replace `u1'=_t if abs(_t/`u1'-1) <= epsdouble()
qui replace `u' =_t if abs(_t/`u' -1) <= epsdouble()
qui stack `treat0' _t _d `adjvars' `strata' ///
`treat0' `u' `du' `adjvars' `strata' ///
`treat0' `u1' `du' `adjvars' `strata', group(3) clear
rename _stack `type'
qui gen int `treat' = 10*(`type'-1)+`treat0'
local name0 "0 as observed"
local name1 "1 as observed"
local name10 "0 if untreated"
local name11 "1 if untreated"
local name20 "0 if fully treated"
local name21 "1 if fully treated"
foreach val in 0 1 10 11 20 21 {
label def treat2 `val' "`name`val''", add
}
label val `treat' treat2
qui stset _t _d, noshow
if "`debug'"=="debug" stsum, by(`treat')
if "`hr'" != "" {
* Output results: table header
di _new as text "Hazard ratios derived from psi=" %9.0g `est'
di `line'
di `col2' "Haz. ratio" `col3' %5.3f "P-value" `col4' "[`level'% Conf. Interval]"
di `line'
if "`strata'"!="" local strataopt strata(`strata')
* ITT analysis
qui stcox `treat0' `adjvars' if `type'==1, `strataopt'
local hritt = exp(_b[`treat0'])
local hrittlow = exp(_b[`treat0']-`zcrit'*_se[`treat0'])
local hrittupp = exp(_b[`treat0']+`zcrit'*_se[`treat0'])
di as text "ITT" as result `col2' `hritt' `col3' %5.3f `PITT' `col4' `hrittlow' `col5' `hrittupp'
* Adjusted analysis
qui stcox `treat0' `adjvars' if inlist(`treat',10,21), `strataopt'
cap local hradj = exp(_b[`treat0'])
local seadj = _b[`treat0']/`ZITT'
local hradjlow = exp(_b[`treat0']-`zcrit'*`seadj')
local hradjupp = exp(_b[`treat0']+`zcrit'*`seadj')
di as text "Adjusted*" as result `col2' `hradj' `col3' %5.3f `PITT' `col4' `hradjlow' `col5' `hradjupp'
* No-treatment analysis
qui stcox `treat0' `adjvars' if inlist(`treat',10,11), `strataopt'
cap local hrnon = exp(_b[`treat0'])
di as text "No treatment**" as result `col2' `hrnon' `col3' %5.3f . `col4' . `col5' .
* Other requested analysis
if "`hr'" != "run" {
tokenize "`hr'"
if "`1'"==substr("observed",1,length("`1'")) local one 1
if "`1'"==substr("untreated",1,length("`1'")) local one 11
if "`1'"==substr("fully",1,length("`1'")) local one 21
if "`2'"==substr("observed",1,length("`2'")) local two 0
if "`2'"==substr("untreated",1,length("`2'")) local two 10
if "`2'"==substr("fully",1,length("`2'")) local two 20
tempvar uservar
qui gen `uservar' = `treat'==`one' if inlist(`treat',`one',`two')
qui stcox `uservar' `adjvars' if inlist(`treat',`one',`two'), `strataopt'
cap local hruser = exp(_b[`uservar'])
local seuser = _b[`uservar']/`ZITT'
local hruserlow = exp(_b[`uservar']-`zcrit'*`seuser')
local hruserupp = exp(_b[`uservar']+`zcrit'*`seuser')
di as text "User***" as result `col2' `hruser' `col3' %5.3f `PITT' `col4' `hruserlow' `col5' `hruserupp'
local note3 `""*** User-specified comparison: " as result "`treat'=`name`one''" as text " vs. " as result "`treat'=`name`two''""'
}
* Output results: table footer
di `line'
di as text "* test-based confidence interval"
di as text "** check on estimation procedure: hazard ratio should be near 1"
if !mi(`"`note3'"') di as text `note3'
di
return scalar HR_ITT = `hritt'
return scalar HR_ITT_low = `hrittlow'
return scalar HR_ITT_upp = `hrittupp'
return scalar HR_adj = `hradj'
return scalar HR_adj_low = `hradjlow'