*! version 1.0.1 29Jul2008
*! Nair (EP) and Hall-Wellner confidence band for survival and cumulative hazard function
program define stcband, rclass sortpreserve
	version 9.0
	syntax [if] [in] , [ GENHIgh(name) GENLOw(name) ///
		  NAIR TLOwer(real -1) TUPper(real -1) NA ///
		  TRANsform(name) Level(integer 95) noGraph * ] // graphic options
	st_is 2 analysis
	marksample touse
	qui replace `touse' = 0 if _st==0
	if "`genhigh'"!="" {
		cap confirm new var `genhigh'
		if _rc {
			di in smcl as err "{p}`genhigh' already defined.{p_end}"
			exit 110
		}
	}
	if "`genlow'"!="" {
		cap confirm new var `genlow'
		if _rc {
			di in smcl as err "{p}`genlow' already defined.{p_end}"
			exit 110
		}
	}
	if "`nair'"!=""	local method "nair"
	else		local method "hw"
	if `tlower' < 0 local tlower
	if `tupper' < 0 local tupper
	su _t if `touse' & _d , meanonly
	if "`tupper'"!="" {
		cap assert `tupper' <= `r(max)'
		if _rc {
			di in smcl as err "{p}Upper time limit is greater than the maximum observed risk time.{p_end}"
			exit 198
		}
	}
	else local tupper = `r(max)'
	if "`tlower'"!="" {
		su _t if `touse' & _d & _t<=`tlower', meanonly
		if "`r(max)'" != ""	local tlower = `r(max)'
		else {
			if "`method'" == "nair" {
				di in smcl as err "{p}Equal precision (Nair) bands does not allow lower time limit less than the minimum observed risk time.{p_end}"
				exit 198
			}
			else local tlower = 0
		}
	}
	else	local tlower = `r(min)'
	qui count if _t>=`tlower' & _t<=`tupper' & `touse' & _d
	if r(N)<1 {
		di in smcl as err "{p}No observations in the range of time of confidence band{p_end}"
		exit 2000
	}
	if "`transform'"!="" {
		if "`transform'" != "log" & "`transform'" != "linear" & "`transform'" != "arcsine" {
			di in smcl as err "{p}transform() option incorrectly specified. Choose: linear, log or arcsine.{p_end}"
			exit 198
		}
	}
	else local transform "log"
	if `level' != 90  & `level' != 95 & `level' != 99 {
			di in smcl as err "{p}Confidence level for simultaneous bands must be either 90, 95 or 99.{p_end}"
			exit 198
	}

	_get_gropts , graphopts(`options') getallowed(LColor LPattern LWidth YLAbel SAving legend scheme xsize ysize plot addplot)
	local options `"`s(graphopts)'"'
	local lpattern `"`s(lpattern)'"'
	local lcolor `"lc(`s(lcolor)')"'
	local lwidth `"lw(`s(lwidth)')"'
	local ylabel `"`s(ylabel)'"'
	local saving `"`s(saving)'"'
	_check4gropts saving, opt(`saving')
	local scheme `"`s(scheme)'"'
	local legend `"`s(legend)'"'
	_check4gropts legend, opt(`legend')
	local plot `"`s(plot)'"'
	local addplot `"`s(addplot)'"'
	local xsize `"`s(xsize)'"'
	local ysize `"`s(ysize)'"'
	if "`saving'" != "" local saving `"saving(`saving')"'
	if "`scheme'" != "" local scheme `"scheme(`scheme')"'
	if "`xsize'" != ""   local xsize `"xsize(`xsize')"'
	if "`ysize'" != ""   local ysize `"ysize(`ysize')"'

	qui count if `touse'
	local n = r(N)
	tempvar s se fa hi_s lo_s time
	tempname al au vperc hperc
	if "`na'" ==""{
		sts gen `s'=s `se'=se(s) if `touse'
		qui replace `se'=`se'^2/`s'^2
	}
	else {
		sts gen `s'=na `se'=se(na) if `touse'
		qui replace `se'=`se'^2
	}
	* Computing al
	g byte `fa' = float(_t)<=float(`tlower') & _d & `touse'
	sort `fa' _t
	scalar `al' = `n' * `se'[_N] / (1+`n'*`se'[_N])
	* Computing au
	qui replace `fa' = float(_t)<=float(`tupper') & _d & `touse'
	sort `fa' _t
	scalar `au' = `n' * `se'[_N] / (1+`n'*`se'[_N])
	sca `vperc' = 1 - (ceil(`au'/2*100)-`au'/2*100)
	sca `hperc' = 1 - (ceil(`al'/2*100)-`al'/2*100)
	qui findcrit `level' `al' `au' `vperc' `hperc' `method'
*	scalar list
	ret scalar cr = `r(cr)'
	if "`na'" == "" {
		if "`method'"=="hw" {
			if "`transform'"=="log" {
				tempvar theta
				qui g double `theta' = exp( (`r(cr)'*(1+`n'*`se')) / (sqrt(`n')*ln(`s')) ) ///
					if `touse' & float(_t)<=float(`tupper') & float(_t)>=float(`tlower')
				qui g `hi_s' = `s'^`theta'
				qui g `lo_s' = `s'^(1/`theta')
			}
			else if "`transform'"=="arcsine" {
				qui g `lo_s' = sin( max(0,asin(sqrt(`s'))     - 0.5*(`r(cr)'*(1+`n'*`se'))/sqrt(`n') * sqrt(`s'/(1-`s'))) )^2 if `touse' ///
					& float(_t)<=float(`tupper') & float(_t)>=float(`tlower')
				qui g `hi_s' = sin( min(_pi/2,asin(sqrt(`s')) + 0.5*(`r(cr)'*(1+`n'*`se'))/sqrt(`n') * sqrt(`s'/(1-`s'))) )^2 if `touse' ///
					& float(_t)<=float(`tupper') & float(_t)>=float(`tlower')
			}
			else if "`transform'"=="linear" {
				qui g `lo_s' = `s' - (`r(cr)'*(1+`n'*`se')/sqrt(`n')) *`s' if `touse' ///
					& float(_t)<=float(`tupper') & float(_t)>=float(`tlower')
				qui g `hi_s' = `s' + (`r(cr)'*(1+`n'*`se')/sqrt(`n')) *`s' if `touse' ///
					& float(_t)<=float(`tupper') & float(_t)>=float(`tlower')
				qui count if `lo_s' < 0
				if `r(N)'>0 {
					qui replace `lo_s' = 0 if `lo_s'<0 & `touse' ///
						& float(_t)<=float(`tupper') & float(_t)>=float(`tlower')
					di _n in smcl in gr `"{p}Some value of low confidence band lower than 0 has been changed to 0.{p_end}"'
				}
				qui count if `hi_s' > 1
				if `r(N)'>0 {
					qui replace `hi_s' = 1 if `hi_s'>1 & `touse' ///
						& float(_t)<=float(`tupper') & float(_t)>=float(`tlower')
					di _n in smcl in gr `"{p}Some value of high confidence band greater than 1 has been reduced to 1.{p_end}"'
				}
			}
		}
		if "`method'"=="nair" {
			if "`transform'"=="log" {
				tempvar theta
				qui g double `theta' = exp( `r(cr)'*sqrt(`se')/ln(`s') ) if `touse' ///
					& float(_t)<=float(`tupper') & float(_t)>=float(`tlower')
				qui g `hi_s' = `s'^`theta'
				qui g `lo_s' = `s'^(1/`theta')
			}
			else if "`transform'"=="arcsine" {
				qui g `lo_s' = sin( max(0,asin(sqrt(`s'))     - 0.5*(`r(cr)'*sqrt(`se')) * sqrt(`s'/(1-`s'))) )^2 if `touse' ///
						& float(_t)<=float(`tupper') & float(_t)>=float(`tlower')
				qui g `hi_s' = sin( min(_pi/2,asin(sqrt(`s')) + 0.5*(`r(cr)'*sqrt(`se')) * sqrt(`s'/(1-`s'))) )^2 if `touse' ///
					& float(_t)<=float(`tupper') & float(_t)>=float(`tlower')
			}
			else if "`transform'"=="linear" {
				qui g `lo_s' = `s' - `r(cr)'*sqrt(`se')*`s' if `touse' ///
					& float(_t)<=float(`tupper') & float(_t)>=float(`tlower')
				qui g `hi_s' = `s' + `r(cr)'*sqrt(`se')*`s' if `touse' ///
					& float(_t)<=float(`tupper') & float(_t)>=float(`tlower')
				qui count if `lo_s' < 0
				if `r(N)'>0 {
					qui replace `lo_s' = 0 if `lo_s'<0 & `touse' ///
						& float(_t)<=float(`tupper') & float(_t)>=float(`tlower')
					di _n in smcl in gr `"{p}Some value of low confidence band lower than 0 has been changed to 0.{p_end}"'
				}
				qui count if `hi_s' > 1
				if `r(N)'>0 {
					qui replace `hi_s' = 1 if `hi_s'>1 & `touse' ///
						& float(_t)<=float(`tupper') & float(_t)>=float(`tlower')
					di _n in smcl in gr `"{p}Some value of high confidence band greater than 1 has been reduced to 1.{p_end}"'
				}
			}
		}
	}

	if "`na'" != "" {
		if "`method'"=="hw" {
			if "`transform'"=="log" {
				tempvar theta
				qui g double `theta' = exp( (`r(cr)'*(1+`n'*`se')) / (sqrt(`n')*`s')  ) ///
					if `touse' & float(_t)<=float(`tupper') & float(_t)>=float(`tlower')
				qui g `hi_s' = `s'*`theta'
				qui g `lo_s' = `s'/`theta'
			}
			else if "`transform'"=="arcsine" {
				qui g `lo_s' = -2*ln( sin( min( _pi/2,asin(exp(-`s'/2)) + 0.5*(`r(cr)'*(1+`n'*`se')/sqrt(`n'))*(exp(`s')-1)^-0.5 ))) if `touse' ///
					& float(_t)<=float(`tupper') & float(_t)>=float(`tlower')
				qui g `hi_s' = -2*ln( sin( max( 0,asin(exp(-`s'/2)) - 0.5*(`r(cr)'*(1+`n'*`se')/sqrt(`n'))*(exp(`s')-1)^-0.5 ))) if `touse' ///
					& float(_t)<=float(`tupper') & float(_t)>=float(`tlower')
			}
			else if "`transform'"=="linear" {
				qui g `lo_s' = `s' - (`r(cr)'*(1+`n'*`se')/sqrt(`n')) if `touse' ///
					& float(_t)<=float(`tupper') & float(_t)>=float(`tlower')
				qui g `hi_s' = `s' + (`r(cr)'*(1+`n'*`se')/sqrt(`n')) if `touse' ///
					& float(_t)<=float(`tupper') & float(_t)>=float(`tlower')
				qui count if `lo_s' < 0
				if `r(N)'>0 {
					qui replace `lo_s' = 0 if `lo_s'<0 & `touse' ///
						& float(_t)<=float(`tupper') & float(_t)>=float(`tlower')
					di _n in smcl in gr `"{p}Some value of low confidence band lower than 0 has been changed to 0.{p_end}"'
				}
			}
		}
		if "`method'"=="nair" {
			if "`transform'"=="log" {
				tempvar theta
				qui g double `theta' = exp( `r(cr)'*sqrt(`se')/`s' ) if `touse' ///
					& float(_t)<=float(`tupper') & float(_t)>=float(`tlower')
				qui g `hi_s' = `s'*`theta'
				qui g `lo_s' = `s'/`theta'
			}
			else if "`transform'"=="arcsine" {
				qui g `lo_s' = -2*ln( sin( min( _pi/2,asin(exp(-`s'/2)) + 0.5*(`r(cr)'*sqrt(`se'))*(exp(`s')-1)^-0.5 ))) if `touse' ///
					& float(_t)<=float(`tupper') & float(_t)>=float(`tlower')
				qui g `hi_s' = -2*ln( sin( max( 0,asin(exp(-`s'/2)) - 0.5*(`r(cr)'*sqrt(`se'))*(exp(`s')-1)^-0.5 ))) if `touse' ///
					& float(_t)<=float(`tupper') & float(_t)>=float(`tlower')
			}
			else if "`transform'"=="linear" {
				qui g `lo_s' = `s' - `r(cr)'*sqrt(`se') if `touse' ///
					& float(_t)<=float(`tupper') & float(_t)>=float(`tlower')
				qui g `hi_s' = `s' + `r(cr)'*sqrt(`se') if `touse' ///
					& float(_t)<=float(`tupper') & float(_t)>=float(`tlower')
				qui count if `lo_s' < 0
				if `r(N)'>0 {
					qui replace `lo_s' = 0 if `lo_s'<0 & `touse' ///
						& float(_t)<=float(`tupper') & float(_t)>=float(`tlower')
					di _n in smcl in gr `"{p}Some value of low confidence band lower than 0 has been changed to 0.{p_end}"'
				}
			}
		}
	}

**To generate varlist
	if "`method'" == "hw" local v "H-W"
	if "`method'" == "nair" local v "EP"
	if "`genhigh'"!="" {
		qui g double `genhigh'  = `hi_s'
		label var `genhigh' `"`v' `level' Upper Band"'
	}
	if "`genlow'"!="" {
		qui g double `genlow'  = `lo_s'
		label var `genlow' `"`v' `level' Lower Band"'
	}

**To graph
	if "`graph'"==""{
		label var `hi_s' `"`v' `level' Upper Band"'
		label var `lo_s' `"`v' `level' Lower Band"'
		if "`scheme'" == "" {
			cap findfile scheme-lean2.scheme
			if !_rc 	local scheme `"scheme(lean2)"'
		}
		if "`lpattern'"=="" {
			local lpattern `"lpattern(l - -)"'
		}
		else local lpattern lp(`lpattern')"'
		if "`xtitle'"=="" {
			local xtitle `"xtitle("Time")"'
		}
		if "`ytitle'"=="" {
			if "`na'" == "" local ytitle `"ytitle("Survival Probability")"'
			else		local ytitle `"ytitle("Cumulative Hazard")"'
		}
		if `"`legend'"'==`""' {
			if "`na'"==""	local legend `"legend(pos(2) col(1) ring(0))"'
			else		local legend `"legend(pos(11) col(1) ring(0))"'
		}
		else local legend `"legend(`legend')"'
		if "`ylabel'"=="" {
			local ylabel `"ylabel(,nogrid)"'
		}
		else local ylabel `"ylabel(`ylabel')"'

		qui g `time' = _t if `touse'
		qui replace `time' = `tlower' if _t<`tlower' & `touse'
		qui replace `time' = `tupper' if _t>`tupper' & `touse'
		qui {
			replace `s' = . if _t > `tupper' & `touse'
			replace `hi_s' = . if _t > `tupper' & `touse'
			replace `lo_s' = . if _t > `tupper' & `touse'
			sort `touse' _t
			replace `s' = `s'[_n-1] if `s'==. & `touse'
			replace `hi_s' = `hi_s'[_n-1] if `hi_s'==. & `touse'
			replace `lo_s' = `lo_s'[_n-1] if `lo_s'==. & `touse'
			replace `s' = . if _t < `tlower' & `touse'
			replace `hi_s' = . if _t < `tlower' & `touse'
			replace `lo_s' = . if _t < `tlower' & `touse'
			if "`na'" ==""	sort `touse' `s' _t
			else	       gsort `touse' -`s' _t
			replace `s' = `s'[_n-1] if `s'==. & `touse'
			replace `hi_s' = `hi_s'[_n-1] if `hi_s'==. & `touse'
			replace `lo_s' = `lo_s'[_n-1] if `lo_s'==. & `touse'
		}
		twoway	(line `s' `hi_s' `lo_s' `time' , sort c(J J J) `lpattern' `lcolor' `lwidth' ///
			`title' `ytitle' `xtitle' `legend' `saving' `scheme' `ylabel' `ysize' `xsize'  `options')  ///
				|| `plot' || `addplot'
	}
end


program define findcrit, rclass
	* Computing critical values from the Hall-Wellner and Nair tables
	args level al au vperc hperc method
	preserve
	tempname au_c al_c sc1 sc2 sc3 sc4 int1 int2 int3 int4 crit
	tempvar c1 c2 c3 c4
	if "`method'"=="hw"	cap findfile HallWellnerTables.dta
	if _rc {
		di _n as err `"HallWellnerTables.dta not found. This file must be in one of the {help adopath:ado-path} directories"'
		exit 601
	}
	if "`method'"=="nair"	cap findfile NairTables.dta
	if _rc {
		di _n as err `"NairTables.dta not found. This file must be in one of the {help adopath:ado-path} directories"'
		exit 601
	}
	local fileful `"`r(fn)'"'
	use `"`fileful'"',clear
	scalar `au_c' = round(`au',.02)
	if  `au_c' > `au'  sca `au_c' = scalar(`au_c') - .02
	sca `al_c' = round(`al',.02)
	if  `al_c' > `al'  sca `al_c' = scalar(`al_c') - .02
	if "`method'"=="hw" {
		g byte `c1' = ( float(lower)==float(scalar(`al_c')) )       & (float(upper)==float(scalar(`au_c')) )
		g byte `c2' = ( float(lower)==float(scalar(`al_c') + .02) ) & (float(upper)==float(scalar(`au_c')) )
		g byte `c3' = ( float(lower)==float(scalar(`al_c')) )       & (float(upper)==float(scalar(`au_c') + .02) )
		g byte `c4' = ( float(lower)==float(scalar(`al_c') + .02) ) & (float(upper)==float(scalar(`au_c') + .02) )
	}
	else {
		tempname al_nc au_nc
		sca `al_nc' = max(`al_c',.02)
		sca `au_nc' = min(`au_c'+.2,.98)
		g byte `c1' = ( float(lower)==float(scalar(`al_nc')) )       & (float(upper)==float(scalar(`au_c'))  )
		g byte `c2' = ( float(lower)==float(scalar(`al_c') + .02) )  & (float(upper)==float(scalar(`au_c'))  )
		g byte `c3' = ( float(lower)==float(scalar(`al_nc')) )       & (float(upper)==float(scalar(`au_nc')) )
		g byte `c4' = ( float(lower)==float(scalar(`al_c') + .02) )  & (float(upper)==float(scalar(`au_nc')) )
	}

	sort `c1'
	scalar `sc1' = k`level'[_N]
	sort `c2'
	scalar `sc2' = k`level'[_N]
	sort `c3'
	scalar `sc3' = k`level'[_N]
	sort `c4'
	scalar `sc4' = k`level'[_N]
	if `sc2' == . scalar `sc2' = (`sc1'+`sc4')/2
*	scalar list `sc1' `sc2' `sc3' `sc4'
	*interpolations
	scalar `int1' = scalar(`sc1') - abs(scalar(`sc1')-scalar(`sc2'))* `hperc'
	scalar `int2' = scalar(`sc4') - abs(scalar(`sc4')-scalar(`sc2'))* `vperc'
	scalar `int3' = scalar(`sc3') - abs(scalar(`sc3')-scalar(`sc4'))* `hperc'
	scalar `int4' = scalar(`sc3') - abs(scalar(`sc3')-scalar(`sc1'))* `vperc'
	scalar `crit' = ( scalar(`int1')*(1-`vperc') + scalar(`int4')*(1-`hperc') + ///
		scalar(`int2')*`hperc' + scalar(`int3')*`vperc' ) / 2
	ret scalar cr = scalar(`crit')
end