*! stcompet 1.0.7 EC 06 NOV 2012 *! Cumulative Incidence in the Presence of Competing Events program define stcompet, sortpreserve version 10 st_is 2 analysis ******* General frame of the command ***************************************** /* Command needs that previous stset has the form: stset survtime, failure(event==# [#]) where # refers to the code of failure of interest. Than stcompet can compute cumulative incidence, standard error and confidence bounds for this event and for other events that you specify in compet# options. Up to 6 competing events can be defined, but there is no problem in adding other compet# options. The initial syntax of the command is like sts gen because similar commands should share similar syntax. */ ****************************************************************************** local rest `"`0'"' gettoken varname rest : rest, parse(" =,") gettoken eqsign rest : rest, parse(" =,") gettoken thing rest : rest, parse(" =,") if `"`eqsign'"' != "=" { error 198 } while `"`eqsign'"' == "=" { confirm new var `varname' local thing = lower(`"`thing'"') if `"`thing'"' == "ci" { local ci `varname' } else if `"`thing'"' == "se" { local se "`varname'" } else if `"`thing'"'=="hi" { local hi "`varname'" } else if `"`thing'"'=="lo" { local lo "`varname'" } else { di in red `"`thing' unknown function"' exit 198 } local 0 `"`rest'"' gettoken varname rest : rest, parse(" =,") gettoken eqsign rest : rest, parse(" =,") gettoken thing rest : rest, parse(" =,") } local what "`se'`hi'`lo'" syntax /* ... */ [if] [in] , COMPET1(numlist missingokay) [ COMPET2(numlist missingokay) /* */ COMPET3(numlist missingokay) COMPET4(numlist missingokay) COMPET5(numlist missingokay) /* */ COMPET6(numlist missingokay) BY(varname) ALLIGNol Level(integer $S_level) ] marksample touse qui replace `touse' = 0 if _st==0 local type : set type set type double preserve *** 1 - Further checks /* previous stset statement must be as required by stcompet */ if "`_dta[st_bd]'"=="" | "`_dta[st_ev]'"=="" { di as err "failure variable must have been specified as failure(varname==numlist) " /* */ _n "on the stset command prior to using this command" exit 198 } local if_st "`_dta[st_ifexp]'" /* it will need to restore existing if expression */ /* enter cannot be specified as varname numlist - 08nov2008 */ if "`_dta[st_enter]'" !="" { local enter `"`_dta[st_enter]'"' gettoken enter : enter, parse(" =") if "`enter'"!= "time" { di as err "enter option cannot be specified as enter(varname==numlist) " /* */ _n "on the stset command prior to using this command" exit 198 } } /* byvar cannot contain non-integer values */ capture confirm numeric variable `by' local isnum = _rc != 7 if `isnum' { capture assert `by' == int(`by') if `touse' if _rc { di as err "`by' contains non-integer values" exit 459 } } /* competing and interest events must be different */ local compet0 "`_dta[st_ev]'" /* main event */ local i = 0 while "`compet`i''" != "" { local c = `i' + 1 while "`compet`c''" != "" { local ulist : list compet`c' & compet`i' if "`ulist'" != "" { di as err "You specify `ulist' as codes for two competing events" exit 459 } local c = `c' + 1 } local i = `i' + 1 } /* level as usual */ if `level'<10 | `level'>99 { di in red "level() invalid" exit 198 } local lvl = invnorm(1-(1-`level'/100)/2) *** 2 - storing all events togheter and number of events. local all_event = "`compet0' `compet1' " local nc = 2 while "`compet`nc''" != "" { local all_event = "`all_event' `compet`nc''" local nc = `nc' + 1 } local nc = `nc' - 1 *** 3 - Cumulative Incidence estimate: sorting bygroups, if defined, and then starting from *** main event and in competing options order quietly { /* initialize as tempvar */ tempvar byuse s_all cuminc gen `cuminc' = . if "`what'" != "" { tempvar n_all d_all n_com d_com sevar gen `sevar' = . } /* If Byvar */ if "`by'" != "" { tempvar byvar egen `byvar' = group(`by') levelsof `byvar' if `touse' /* byvar values */ local byvalue "`r(levels)'" gen byte `byuse' = 0 foreach X of local byvalue { replace `byuse' = `touse' & `byvar'==`X' streset if `byuse', f(`_dta[st_bd]'==`all_event') sts gen `s_all' = s if "`what'" != "" { sts gen `n_all' = n `d_all' = d } local i = 0 while "`compet`i''" != "" { tempvar I`X'`i' S`X'`i' gen `I`X'`i'' = . /* compet compute cumulative incidence - It resets the code of the events in failvar */ compet `byuse' `s_all' `I`X'`i'', fail(`compet`i'') if "`what'" != "" { /* some n and some d must be missing because of ties */ bysort `byuse' _t (`I`X'`i''): /* */ gen `n_com' = `n_all' if _n==1 sts gen `d_com' = d replace `d_com' = . if `n_com'==. gen `S`X'`i'' = . /* secomp compute standard errors and store them in sevar */ secomp `byuse' `I`X'`i'' `n_com' `d_all' /* */ `d_com' `s_all' `S`X'`i'' drop `d_com' `n_com' replace `sevar' = `S`X'`i'' if `S`X'`i''!=. } /* accumulate estimate in cuminc */ replace `cuminc' = `I`X'`i'' if `I`X'`i''!=. local i = `i' + 1 } drop `s_all' `n_all' `d_all' } } /* No byvar */ else { streset if `touse', f(`_dta[st_bd]'==`all_event') sts gen `s_all' = s if "`what'" != "" { sts gen `n_all' = n `d_all' = d } forvalues i = 0(1)`nc' { tempvar I`i' S`i' gen `I`i'' = . compet `touse' `s_all' `I`i'' , fail(`compet`i'') if "`what'" != ""{ bysort `touse' _t (`I`i''): /* */ gen `n_com' = `n_all' if _n==1 sts gen `d_com' = d replace `d_com' = . if `n_com'==. gen `S`i'' = . secomp `touse' `I`i'' `n_com' `d_all' `d_com' `s_all' `S`i'' drop `d_com' `n_com' replace `sevar' = `S`i'' if `S`i''!=. } replace `cuminc' = `I`i'' if `I`i''!=. } } ***4 - Saving variables if "`ci'" != "" { gen `ci' = `cuminc' label var `ci' "Cumulative Incidence" } if "`se'" != "" { gen `se' = `sevar' label var `se' "Stand Err Cum Incidence" } if "`hi'" != "" { /* log(-log) confidence bouds -> JB Choudhury, Stat in Med 2002; 21: 1129 */ if "`allignol'" == "" { gen `hi' = `cuminc' ^ exp(`lvl'*`sevar'/(`cuminc' *log(`cuminc' ))) label var `hi' "CI() `level'% upper bound" } /* log(-log) confidence bouds -> Beyersmann J et al. Competing Risks and Multistate Models with R, p. 62 */ else { gen `hi' = 1 - (1-`cuminc')^exp(-`lvl'*`sevar'/((1-`cuminc')*log((1-`cuminc' )))) label var `hi' "CI() `level'% upper bound (Allignol)" } } if "`lo'" != "" { if "`allignol'" == "" { gen `lo' = `cuminc' ^ exp(-`lvl'*`sevar'/(`cuminc' *log(`cuminc' ))) label var `lo' "CI() `level'% lower bound" } else { gen `lo' = 1 - (1-`cuminc')^exp(`lvl'*`sevar'/((1-`cuminc')*log((1-`cuminc' )))) label var `lo' "CI() `level'% lower bound (Allignol)" } } set type `type' char _dta[st_ifexp] `if_st' streset, f(`_dta[st_bd]'==`compet0') } restore,not end program define compet /* Compute Crude Cumulative Incidence */ version 7.0 syntax varlist, fail(numlist missingokay) gettoken byuse varlist : varlist gettoken all_surv varlist : varlist tempvar h_comp is_even all_comp gen byte `is_even' = 0 foreach num of local fail { replace `is_even' = 1 if `_dta[st_bd]'==`num' & `byuse' if "`_dta[st_exit]'" != "" { replace `is_even'=0 if `_dta[st_bt]' > `_dta[st_exexp]' & `byuse' } } count if `is_even' if `r(N)'== 0 { exit } streset if `byuse', f(`_dta[st_bd]'==`fail') sts gen `h_comp' = h bysort `all_surv' (`is_even') : gen double `all_comp' = `all_surv' if _n==_N replace `h_comp' = . if `all_comp' == . gsort -`all_comp' replace `varlist' = cond(_n!=1,`h_comp' * `all_comp'[_n-1],`h_comp') gsort -`byuse' _t `varlist' replace `varlist' = sum(`varlist') if `is_even' end program define secomp /* Compute Standard Error - Marubini & Valsecchi p.341 1995 */ version 7 args by Cum n d_all d_com s sevar tempvar newby CI first second third tempname t_CI gen byte `newby' = `by' & `n'!=. /* newby is 1 if byuse e n not missing */ gsort -`newby' _t gen `CI' = cond(_n==1 & `Cum'==.,0,`Cum') replace `CI' = `CI'[_n-1] if `CI'==. & `newby' count if `newby' forvalues i = 1 / `r(N)' { sca `t_CI' = `CI'[`i'] gen `first' = (`t_CI' - `CI')^2*`d_all'/(`n'*(`n'-`d_all')) gen `second' = cond(_n>1,`s'[_n-1]^2*`d_com'*(`n'-`d_com')/`n'^3, /* */ `d_com' * (`n'-`d_com') / `n'^3) if _n<=`i' gen `third' = cond(_n>1,(`t_CI' - `CI')*`s'[_n-1]*`d_com'/`n'^2, /* */ (`t_CI' - `CI') * `d_com' / `n'^2) if _n<=`i' replace `first' = sum(`first') if _n<=`i' replace `second' = sum(`second') if _n<=`i' replace `third' = sum(`third') if _n<=`i' replace `sevar' = sqrt(`first' + `second' - 2*`third') in `i' drop `first' `second' `third' } replace `sevar' = . if `Cum' ==. /* filling missing values with previous estimates */ sort `Cum' `sevar' replace `sevar' = `sevar'[_n-1] if `sevar'==. & `Cum'!=. end