*! stpepemori 1.0.3 EC 02 OCT 2010
*! Testing Cumulative Incidence and Conditional Probability

program define stpepemori, sortpreserve rclass
	version 9.0
	st_is 2 analysis
	syntax varname [if] [in] , COMPET(numlist missingokay) [ CONDitional ]
	marksample touse
	qui replace `touse' = 0 if _st==0
		/* it will need to restore existing if expression */
	local if_st "`_dta[st_ifexp]'" /* it will need to restore existing if expression */
	preserve
	qui keep if `touse'
	
*** 1 - Checks
		/* previous stset statement must be as required by stcompet */
	if "`_dta[st_bd]'"=="" | "`_dta[st_ev]'"=="" {   
	    di as err  "{p}failure variable must have been specified as failure(varname==numlist) " /*
        */ _n "on the stset command prior to using this command.{p_end}"
		exit 198
        }
	if `"`_dta[st_id]'"' != "" {
		cap bysort `_dta[st_id]' : assert _N==1
		if _rc {
			di as err  "stpepemori requires that the data are stset with only one observation per individual"
			exit 198
		}   
	}
		/* Test compares just two groups */
	tempvar bybin
	qui bysort `varlist' : g byte `bybin' = _n==1
	qui replace `bybin' = sum(`bybin') - 1
	if `bybin'[_N] != 1 {
		di as err "Pepe and Mori test compares the cumulative incidence of just two groups."
		exit 459
	}
		/* Competing and interest events must be different */
	local compet0 "`_dta[st_ev]'"   /* main event */ 
	local intlist : list compet0 & compet
	if "`intlist'" != "" {
		di as err "You specify `intlist' as codes for two competing events"
                  	exit 459
        }

	local orev "`_dta[st_bd]'"
	tempvar cens w_1 w_2 S_1 S_2 S_3 sigma1 sigma2 n0 n1 byuse 
	quietly {
*** 2 - Create -the appropriate censoring variable : 0 = Alive, 1 = Main event, 2 = Competing event -n_risk
*		g byte `cens' = (`_dta[st_bd]'==`compet0') + 2*(`_dta[st_bd]'==`compet')
		recode `_dta[st_bd]' (`compet0' = 1) (`compet' = 2), gen(`cens')
		streset, f(`cens' == 0 1 2)
		sts gen `n0' = n , by(`bybin')
		g long `n1' = `n0' if `bybin'==1
		replace `n0' = .   if `bybin'==1

*** 3 - Cumulative Incidence estimate and censoring distribution for event 1 and 2
		g byte `byuse' = 0 
		forval X = 0/1 {
			replace `byuse' = (`bybin'==`X')
			tempvar s_all`X' d`X' I`X'_1 I`X'_2 C`X'_1 C`X'_2 d`X'_1 d`X'_2
			streset if `byuse', f(`cens'== 1 2)
			sts gen `s_all`X'' = s  
			g byte  `d`X''   = _d
			replace `d`X''   = 0 if `d`X''==.
			g  double `I`X'_1' = .
			g  double `I`X'_2' = .
			g byte `d`X'_1' = .
			g byte `d`X'_2' = .
			/* compet compute cumulative incidence - It resets the code 
			   of the events in failvar */
			compet `byuse' `s_all`X'' `d`X'_1' `I`X'_1' , fail(1) 
			replace `d`X'_1'   = 0 if `d`X'_1'==.
			compet `byuse' `s_all`X'' `d`X'_2' `I`X'_2' , fail(2) 
			replace `d`X'_2'   = 0 if `d`X'_2'==.
			streset if `byuse', f(`cens'== 0 2)
			sts gen `C`X'_1' = s
			streset if `byuse', f(`cens'== 0 1)
			sts gen `C`X'_2' = s
			count if `byuse'
			local N`X' = `r(N)'
		}
*** 4 - Adapting basic quantities : Cum. Inc. , Cond. Prob. and Censoring Distrbution
		char _dta[st_ifexp] `if_st'
		streset, f(`cens'== 1 2)
		foreach var of varlist `I0_1' `I1_1' `I0_2' `I1_2' {
			sort _t `var'
			replace `var' = 0 if _n == 1  & `var'== .
			replace `var' = `var'[_n-1] if `var'== . & _n > 1 
		}	
		foreach var of varlist `C0_1' `C1_1' `C0_2' `C1_2' `s_all0' `s_all1' {
			sort _t `var'
			replace `var' = 1 if _n == 1  & `var'== .
			replace `var' = `var'[_n-1] if `var'== . & _n > 1 
		}
		if "`conditional'" != "" {
			tempvar CP0_1 CP0_2 CP1_1 CP1_2
			g  double `CP0_1' = `I0_1' / (1 - `I0_2') 
			g  double `CP0_2' = `I0_2' / (1 - `I0_1') 
			g  double `CP1_1' = `I1_1' / (1 - `I1_2') 
			g  double `CP1_2' = `I1_2' / (1 - `I1_1') 
		}

		/* Max time at the minimun last time between groups */
		su _t if !`bybin', meanonly
		local tau = `r(max)'
		su _t if `bybin', meanonly
		local tau = min(`r(max)',`tau')
		keep if _t <= `tau'
					
		collapse `I0_1' `I1_1' `I0_2' `I1_2' `C0_1' `C1_1' `C0_2' `C1_2' `CP0_1' `CP0_2' `CP1_1' `CP1_2' `s_all0' `s_all1'  ///
			`n0' `n1' (sum) `d0' `d0_1' `d0_2' `d1' `d1_1' `d1_2', by(_t) fast

*** 5 - Weights and cumulative weighted differences
		g `w_1' = cond(_n>1,((`N0'+`N1')*`C0_1'[_n-1]*`C1_1'[_n-1]) / (`N0'*`C0_1'[_n-1] + `N1'*`C1_1'[_n-1]),1) 
		g `w_2' = cond(_n>1,((`N0'+`N1')*`C0_2'[_n-1]*`C1_2'[_n-1]) / (`N0'*`C0_2'[_n-1] + `N1'*`C1_2'[_n-1]),1) 
		if "`conditional'" == "" {
			g `S_1' = sum((`I0_1'-`I1_1') * `w_1' * (_t[_n+1] - _t[_n]) ) 
			g `S_2' = sum((`I0_2'-`I1_2') * `w_2' * (_t[_n+1] - _t[_n]) ) 
		}
		else {
			g `S_1' = sum((`CP0_1'-`CP1_1') * `w_1' * (_t[_n+1] - _t[_n]) ) 
			g `S_2' = sum((`CP0_2'-`CP1_2') * `w_2' * (_t[_n+1] - _t[_n]) ) 
		}
		local s1 = (sqrt(`N0'*`N1'/(`N0'+`N1')) * `S_1'[_N] )^2
		local s2 = (sqrt(`N0'*`N1'/(`N0'+`N1')) * `S_2'[_N] )^2
		sort _t 
		replace `n0' = `n0'[_n-1] - `d0'[_n-1] if `n0'==. 
		replace `n1' = `n1'[_n-1] - `d1'[_n-1] if `n1'==. 
*** 6 - Variance
		forval i = 1/2 {
* Group 1
			if `i' == 1 local cmp = 2
			else        local cmp = 1
			if "`conditional'" == "" {
				replace `S_1' = sum(`w_`i''*(_t[_n+1] - _t[_n])*(1-`I0_`i''))  
				su `S_1' , meanonly
				replace `S_1' = `r(max)' - `S_1' +  `w_`i''*(_t[_n+1] - _t[_n])*(1-`I0_`i'')	// nu_1
				replace `S_2' = sum(`w_`i''*(_t[_n+1] - _t[_n]))  
				su `S_2' , meanonly
				replace `S_2' = `r(max)' - `S_2' +  `w_`i''*(_t[_n+1] - _t[_n])			// nu_2
				g `S_3' = sum(`w_`i''*(_t[_n+1] - _t[_n])*`I0_`i'')
				su `S_3' , meanonly
				replace `S_3' = `r(max)' - `S_3' +  `w_`i''*(_t[_n+1] - _t[_n])*`I0_`i''	
				g `sigma1'    = ((`S_1' - `I0_`cmp''*`S_2')^2 * `d0_`i'' + `S_3'^2*(`d0' - `d0_`i'')) / (`n0'*(`n0'-1))
				su `sigma1',meanonly
			}
			else {
				replace `S_1' = sum(`w_`i''*(_t[_n+1]-_t[_n])*`s_all0'/(1-`I0_`cmp'')^2)  
				su `S_1' , meanonly
				replace `S_1' = `r(max)' - `S_1' +  `w_`i''*(_t[_n+1]-_t[_n])*`s_all0'/(1-`I0_`cmp'')^2	
				g `sigma1' = `S_1'^2*((1-`I0_`cmp'')^2 *`d0_`i''+(`d0'-`d0_`i'')*`I0_`i''^2) / (`n0'*(`n0'-1))
				su `sigma1',meanonly
			}
			local var1_1 = `r(sum)'
* Group 2
			if "`conditional'" == "" {
				replace `S_1' = sum(`w_`i''*(_t[_n+1] - _t[_n])*(1-`I1_`i''))     		
				su `S_1' , meanonly
				replace `S_1' = `r(max)' - `S_1' +  `w_`i''*(_t[_n+1] - _t[_n])*(1-`I1_`i'')	
				replace `S_2' = sum(`w_`i''*(_t[_n+1] - _t[_n]))	 	
				su `S_2' , meanonly
				replace `S_2' = `r(max)' - `S_2' +  `w_`i''*(_t[_n+1] - _t[_n])			
				replace `S_3' = sum(`w_`i''*(_t[_n+1] - _t[_n])*`I1_`i'') 
				su `S_3' , meanonly
				replace `S_3' = `r(max)' - `S_3' +  `w_`i''*(_t[_n+1] - _t[_n])*`I1_`i''	
				g `sigma2'    = ((`S_1' - `I1_`cmp''*`S_2')^2 * `d1_`i'' + `S_3'^2*(`d1' - `d1_`i'')) / (`n1'*(`n1'-1))    
				su `sigma2',meanonly
				drop `S_3' 
			}
			else {
				replace `S_1' = sum(`w_`i''*(_t[_n+1]-_t[_n])*`s_all1'/(1-`I1_`cmp'')^2)  
				su `S_1' , meanonly
				replace `S_1' = `r(max)' - `S_1' +  `w_`i''*(_t[_n+1]-_t[_n])*`s_all1'/(1-`I1_`cmp'')^2	
				g `sigma2' = `S_1'^2*((1-`I1_`cmp'')^2 *`d1_`i''+(`d1'-`d1_`i'')*`I1_`i''^2) / (`n1'*(`n1'-1))  
				su `sigma2',meanonly
			}
			local var1_2 = `r(sum)'
			local var1 = (`N0'*`N1'*(`var1_1' + `var1_2')) / (`N0' + `N1')
			local z`i' = `s`i'' / `var1'
			local p`i' = chiprob(1,`z`i'')
			drop `sigma1' `sigma2'
		}
	}
*** 7 - Show and return results for two types of event
	local tst "cumulative incidence"
	if "`conditional'" != "" local tst "conditional probability"
	ret scalar r1    = `z1'
	ret scalar r2    = `z2'
	ret scalar p1    = `p1'
	ret scalar p2    = `p2'
	local rnd1 "0000"
	local rnd2 "0000"

	foreach i of numlist 1 10 100 1000 10000{
		if `z1' > `i' local rnd1 = substr("`rnd1'",2,.)
		if `z2' > `i' local rnd2 = substr("`rnd2'",2,.)
	}
	foreach i in p1 p2 {
		local eq`i' "= "
		if ``i'' < 0.00001 {
			local `i' ".00001"
			local eq`i' "< "
		}
		else	local `i' = round(``i'',.00001)
	}
	di _n(2) in smcl in gr `"{title:Pepe and Mori test comparing the `tst' of two groups of `varlist'}"'
*	di as res "Pepe and Mori test comparing the `tst' of two groups of `varlist' "
	di _n in gr "     Main event failure:  " in ye "`orev' == `compet0'"
	di in gr "Chi2(1) = " in ye round(`z1',.`rnd1'1)  in gr "  -  p `eqp1' "  in ye "0" `p1'
	di
	di in gr "Competing event failure:  " in ye "`orev' == `compet'"
	di in gr "Chi2(1) = " in ye round(`z2',.`rnd2'1)  in gr "  -  p `eqp2' "  in ye "0" `p2'
end


program define compet    /* Compute Crude Cumulative Incidence */ 
	version 9
	syntax varlist, fail(numlist missingokay) 
	gettoken byuse varlist : varlist
	gettoken all_surv varlist : varlist
	gettoken d varlist : varlist
	tempvar h_comp is_even all_comp
        gen byte `is_even' = 0
        replace `is_even' = 1 if `_dta[st_bd]'==`fail' 
	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')
	replace `d' = _d 
	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' & `byuse'
end