*! version 1.0.0  20 July 98, Kim Lyngby Mikkelsen, kimmik@doktor.dk
program define stcumh
/* Syntax: strata_var baseh_var,
		[STRCode(x,y)] [noCurtail] [Keep(n)] [SAVE(outsheet)] */
version 5.0
st_is
local varlist "req ex min(2) max(2)"
local if "opt"
local in "opt"
local options "STRCode(string) Keep(string) noCurtail save(string) * "
local options "`options' T1title(string) T2title(string) * "
local options "`options' B2title(string) L1title(string) * "
parse "`*'"
parse "`varlist'", parse(" ")
local level `1'
local baseh `2'

* Extract st data
local time : char _dta[st_t]
local aw   : char _dta[st_w]

* Check for weights
if "`aw'"!="" {
       di in re "weighted analysis not implemented"
       exit
   }
di "'stcumh': Graphical 'check' of the proportional hazards assumption."
st_show `show'
di in green "   stratified by:" in yellow "  `level'"
di ""

preserve

* 'if() and/or in()'-restriction.
tempvar touse
mark `touse' `if' `in'
local n_all1 = [_N]
qui keep if `touse'
local n_res1 = [_N]
local res = `n_all1' - `n_res1'
if `res' > 0	{
di "`res' observation dropped due to 'if() + in()'-restriction!"
}
quietly tab `level'
local nstrata = _result(2)
/* 				If more than two strata in strata_var*/
if `nstrata'>2	{
		if "`strcode'"==""		{
di in red "When more than two strata in the strata-specifying "
di in red "variable '`level'', the two strata to be plotted must"
di in red "be specified in the 'STRCode(x,y)'-option!"
				exit 101
		}
		else	{
			parse "`strcode'", parse(",")
			capture confirm integer number `1'
			if _rc!=0 {
		di in red "`1'? Strata-indicators must be integers"
					exit 101
			}
			local scx `1'
			capture confirm integer number `3'
			if _rc!=0 {
		di in red "`3'? Strata-indicators must be integers"
					exit 101
			}
			local scy `3'
			if "`4'" !=""	{
di in red "Two, and only two strata are allowed in: , strc(x,y)"
						exit 101
				}
				local gate "="
		 }
}
/*	If two strata in strata_var*/
else if	`nstrata'==2	{
		qui su `level'
		local scx = _result(5)
		local scy = _result(6)
		local gate "="
}
/*	If one stratum in strata_var*/
else if	`nstrata'==1	{
di in red "The strata-specifying variable '`level''"
di in red "must contain at least two strata!"
exit 101
}
/*	Calculate cumulative hazards */
keep `time' `level' `baseh'
local n_all = [_N]
qui keep if `baseh'!=.
local n_endp = [_N]
local n_cens = `n_all' - `n_endp'
#d ;
di "`n_cens' " in red "censored" in yellow "
(of `n_all') observations were dropped." ;
#d cr
di ""
sort `level' `time'
tempvar cumh
quietly by `level': gen cumh=sum(`baseh')
tempvar cumh0 cumh1
gen `cumh0'=0
gen `cumh1'=0
*	lab var `cumh0' "Stratum coded `scx'"
*	lab var `cumh1' "Stratum coded `scy'"
sort `time'
quietly {
	replace `cumh0'=cumh if `level'==`scx'
	replace `cumh1'=cumh if `level'==`scy'
	replace `cumh0'=`cumh0'[_n-1] if `cumh0'==0 & _n>1
	replace `cumh1'=`cumh1'[_n-1] if `cumh1'==0 & _n>1
}
/* Titles		*/
if "`t1title'"==""	{
#d ;
local options "`options' t1(Graphical check of the
proportional hazards assumptions.)" ;
#d cr
}
else	local options "t1(`t1title') `options'"
if "`t2title'"!=""	{
		di "The t2title (`t2title') is ignored!"
}
if "`b2title'"!=""	{
		local options "`options' b2(`b2title')"
}
else	local options "`options' b2(Stratum coded `scx')"
if "`l1title'"!=""	{
		local options "`options' l1(`l1title')"
}
else	local options "`options' l1(Stratum coded `scy')"

if "`keep'"!=""		{
		capture confirm integer number `keep'
		if _rc!=0 {
				di in red "keep(`keep')? `keep' must be an integer"
				exit 101
		}
		if `keep' < 1	{
		di in red "keep(`keep')? The integer in keep() must be > 1"
				exit 101
		}
		if 	"`curtail'"==""	{
#d ;
		di in red "When 'keep(n)' is specified,
so must be the 'nocurtale' option!" ;
#d cr
				exit 101
		}
*			else
			qui su `time' if `level'==`scx'
			local xmin = round(_result(5),0.1)
			local xmax = round(_result(6),0.1)
			di "Time range for stratum coded `scx': `xmin'-`xmax'"
			qui su `time' if `level'==`scy'
			local ymin = round(_result(5),0.1)
			local ymax = round(_result(6),0.1)
			di "Time range for stratum coded `scy': `ymin'-`ymax'"
					local d  = [_N]
				qui su `cumh0'
qui drop if `cumh0'[_n]==`cumh0'[_n-`keep'] & `cumh0'[_n] == `cumh0'[_N]
				local d0 = [_N]
				qui su `cumh0'
qui drop if `cumh1'[_n]==`cumh1'[_n-`keep'] & `cumh1'[_n] == `cumh1'[_N]
				local d1 = [_N]
				local drop = `d' - `d1'
	di "While `keep' observations were kept, `drop' were dropped."
				qui su `time' if `level'==`scx'
			local xmin = round(_result(5),0.1)
			local xmax = round(_result(6),0.1)
			qui su `time' if `level'==`scy'
			local ymin = round(_result(5),0.1)
			local ymax = round(_result(6),0.1)
			local t_min = min(`xmin',`ymin')
			local t_max = max(`xmax',`ymax')
#d ;
di in red "'Un-curtailed - keep'" in yellow
" time range:        `t_min'-`t_max'" ;
#d cr
local options "`options' t2(Plotted time range: `t_min'-`t_max')"
			local restric = 1
		local gate "="
}
if "`curtail'" `gate'=""	{
		qui su `time' if `level'==`scx'
		local xmin = round(_result(5),0.1)
		local xmax = round(_result(6),0.1)
		di "Time range for stratum coded `scx': `xmin'-`xmax'"
		qui su `time' if `level'==`scy'
		local ymin = round(_result(5),0.1)
		local ymax = round(_result(6),0.1)
		di "Time range for stratum coded `scy': `ymin'-`ymax'"
		local t_min = max(`xmin',`ymin')
		local t_max = min(`xmax',`ymax')
di in red "'Curtailed'" in yellow " time range:        `t_min'-`t_max'"
	local options "`options' t2(Plotted time range: `t_min'-`t_max')"
		local N = [_N]
		di "'Un-curtailed' number of endpoints: `N'"
		qui drop if `time'<`t_min' | `time'>`t_max'
		local N = [_N]
		di "'Curtailed' number of endpoints:    `N'"
		qui su `cumh0'
		local cumh0tm=_result(6)
		qui su `cumh1'
		local cumh1tm=_result(6)
		if `cumh0tm' >= `cumh1tm'	{
				local s_max = round(`cumh0tm',.05)
		}
		else	{
				local s_max = round(`cumh1tm',.05)
		}
		local a = .2*`s_max'
		local b = .4*`s_max'
		local c = .6*`s_max'
		local d = .8*`s_max'
		local e = `s_max'
		local lab "0,`a',`b',`c',`d',`e'"
#d ;
		graph `cumh1' `cumh0', xscale(0,`s_max') yscale(0,`s_max')
		xlab(`lab') ylab(`lab') s(.) c(l) `options'	;
#d cr
*			reg `cumh1' `cumh0'
}
else	{
		if "`keep'"==""		{
		qui su `time' if `level'==`scx'
		local xmin = round(_result(5),0.1)
		local xmax = round(_result(6),0.1)
	di "Time range for stratum coded `scx': `xmin'-`xmax'"
		qui su `time' if `level'==`scy'
		local ymin = round(_result(5),0.1)
		local ymax = round(_result(6),0.1)
	di "Time range for stratum coded `scy': `ymin'-`ymax'"
		local t_min = min(`xmin',`ymin')
		local t_max = max(`xmax',`ymax')
di in red "'Plotted'" in yellow " time range:          `t_min'-`t_max'"
	local options "`options' t2(Plotted time range: `t_min'-`t_max')"
		}
		qui su `cumh0'
		local cumh0tm=_result(6)
		qui su `cumh1'
		local cumh1tm=_result(6)
		if `cumh0tm' >= `cumh1tm'	{
				local s_max = round(`cumh0tm',.05)
		}
		else	{
				local s_max = round(`cumh1tm',.05)
		}
		local a = .2*`s_max'
		local b = .4*`s_max'
		local c = .6*`s_max'
		local d = .8*`s_max'
		local e = `s_max'
		local lab "0,`a',`b',`c',`d',`e'"
#d ;
		graph `cumh1' `cumh0', xscale(0,`s_max') yscale(0,`s_max')
		xlab(`lab') ylab(`lab') s(.) c(l) `options' ;
#d cr
}
if "`save'"!=""	{
		local w = length("`save'")
		if `w'>8	{
#d ;
di in red "Data not saved, as filename `save'
has more than 8 characters" ;
#d cr
				exit 101
		}
		else	{
				outs `time' `cumh0' `cumh1' using `save'
				di in red "Data saved in file: `save'"
		}
}
end