clear set mem 10m // Example 1 set obs 5000 set seed 1234 gen t25 = -1/.25 * log(1-uniform()) gen t99 = -1/.99 * log(1-uniform()) gen time = min(t99,t25) gen byte fail = (t25=t99) replace fail = 0 if time>=2 replace time = 2 if time>2 stset time, f(fail==1) noshow sts gen KM = s gen Complement = 1- KM stcompet CumInc=ci, compet1(2) gen CI=CumInc if fail==1 twoway line Complement CI time, sort xla(0(.5)2) yla(0(.1).5) clw(thick thick) /// yti("Probability") legend(off) text(.31 1.2 "1-KM", place(e)) /// text(.15 1.6 "CI", place(e)) scheme(sj) more set obs 10000 gen treat = _n>5000 replace t25 = -1/.25 * log(1-uniform()) replace t99 = -1/.495 * log(1-uniform()) in 5000/l replace time = min(t99,t25) replace fail = (t25=t99) replace fail = 0 if time>=2 replace time = 2 if time>2 stset time, f(fail==1) noshow stcompet Newci=ci,compet1(2) by(treat) gen CI_treat1=Newci if treat==1 gen CI_treat0=Newci if treat==0 label var CI_treat1 "Treated" label var CI_treat0 "Untreated" label define fail 1 "fail 1" 2 "fail 2" label val fail fail twoway line CI_* time if fail, sort yti("Cumulative Incidence") xti("time") xla(0(.5)2) yla(#5) /// by(fail, style(compact) rows(1) yrescale legend(pos(12)) note("")) /// subtitle(, ring(0) pos(11) nobexpand) clp("-##" l) clw(thick thick) scheme(sj) more // Example 2 use "prostatecancer", clear stset time,f(status==1) stcompet CI = ci hilim = hi lowlim = lo, compet1(2) compet2(3) by(treatment) gen CI_tr = CI if treat label var CI_tr "Treated" gen CI_untr = CI if !treat label var CI_untr "Untreated" gen byte cause=status label val cause fail // The following lines aim to extend the last cumulative incidence estimate // of each competing cause of death to the greatest observed time in each treatment arm. // That makes sense when a censoring occurs at the last observed time as in the data at hand. gen long Months=time qui su Months if treatment, meanonly local max = r(max) gsort treatment -cause time // This sorting instruction puts the censored observations of treated patients // to the end of data set. replace Months= `max' in -3/l // In the last three censored treated patients the value of // var "Months" is replaced with the greatest observed time of the // treated patients. Of course we have checked that there are // at least three censored observations for each treatment arm. // If there weren't, we would have to add some record to the data. replace cause=1 in l // The codes relating to each competing event are substituted in these last three records. replace cause=2 in -2 // Negative numbers may be used to specify a record starting from the end of the data. replace cause=3 in -3 // Note that in the case of censored observations the cumulative incidence estimate is missing. // So now we have for each cause of death and treatment arm one record where // Months has its greatest value and the cumulative incidence estimate is missing. qui su Months if !treatment, meanonly // Previous instructions are repeated for the untreated patients local max = r(max) gsort -treatment -cause time replace Months= `max' in -3/l replace cause=1 in l replace cause=2 in -2 replace cause=3 in -3 // Using bysort prefix as in the lines below, in each group of treatment // and cause of death the missing cumulative incidence estimate // of the last record is replaced with the previous one. bysort treatment cause (CI_tr) : replace CI_tr=CI_tr[_n-1] if _n==_N bysort treatment cause (CI_untr) : replace CI_untr=CI_untr[_n-1] if _n==_N twoway line CI_* Months if cause, sort c(J J) by(cause, note("") style(compact) rows(1) /// legend(pos(12))) subtitle(, ring(0) pos(11) nobexpand) clw(thick thick) clp(-# l) /// yla(0(.2).6) xla(0(20)80) yti("Cumulative Incidence") legend(rows(1)) scheme(sj) more gen hi_tr = hilim if treat gen hi_untr = hilim if !treat bysort treatment cause (hi_tr) : replace hi_tr=hi_tr[_n-1] if _n==_N bysort treatment cause (hi_untr) : replace hi_untr=hi_untr[_n-1] if _n==_N gen lo_tr = lowlim if treat gen lo_untr = lowlim if !treat bysort treatment cause (lo_tr) : replace lo_tr=lo_tr[_n-1] if _n==_N bysort treatment cause (lo_untr) : replace lo_untr=lo_untr[_n-1] if _n==_N twoway line CI_* hi_* lo_* time if cause==1, sort c(J J J J J J) /// clp(l -## l -## l -##) clw(medthick medthick medthick medthick medthick medthick) /// clc(black black black black black black) /// yla(0(.2).6) xla(0(20)80) yti("Cancer Cumulative Incidence") ysca(range(0 .6)) /// legend(order(1 2) ring(0) pos(11) cols(1)) scheme(sj) ****** Mary Lunn Analysis stset,clear expand 3 if touse bysort id : gen byte n = _n if touse gen byte deltacvd = n==2 if touse gen byte deltaother = n==3 if touse rename status fail gen byte status=fail==1 & n==1 | fail==2 & n==2 | fail==3 & n==3 if touse list id time fail status delta* if id==1 | id==4 | id==6 | id==32,sepby(id) abbreviate(10) more stset time if touse,f(status) gen byte treat_canc = (1-deltacvd)*(1-deltaother)*treatment gen byte treat_cvd = deltacvd*treatment gen byte treat_other = deltaother*treatment gen byte AG_canc = (1-deltacvd)*(1-deltaother)*AG gen byte AG_cvd = deltacvd*AG gen byte AG_other = deltaother*AG gen byte WT_canc = (1-deltacvd)*(1-deltaother)*WT gen byte WT_cvd = deltacvd*WT gen byte WT_other = deltaother*WT gen byte PF_canc = (1-deltacvd)*(1-deltaother)*PF gen byte PF_cvd = deltacvd*PF gen byte PF_other = deltaother*PF gen HX_canc = (1-deltacvd)*(1-deltaother)*HX gen byte HX_cvd = deltacvd*HX gen byte HX_other = deltaother*HX gen byte HG_canc = (1-deltacvd)*(1-deltaother)*HG gen byte HG_cvd = deltacvd*HG gen byte HG_other = deltaother*HG gen byte SZ_canc=(1-deltacvd)*(1-deltaother)*SZ gen byte SZ_cvd = deltacvd*SZ gen byte SZ_other = deltaother*SZ gen byte SG_canc = (1-deltacvd)*(1-deltaother)*SG gen byte SG_cvd = deltacvd*SG gen byte SG_other = deltaother*SG *****table 7 di _n as text "Table 7 Mary Lunn" stcox delta* treat_* AG_* WT_* PF_* HX_* HG_* SZ_* SG_* , nohr nolog exit