*! version 2.0.1  24feb2000 
*! syntax [if] [in] [, Quant(integer 50) BY(varname) Level(integer $S_level) ]
program define stquant, rclass
	version 6
        st_is 2 analysis
        /* By option added */
        syntax [if] [in] [, Quant(integer 50) BY(varname) Level(integer $S_level) ]
        tempvar touse
        st_smpl `touse' `"`if'"' `"`in'"' `"`by'"' `""'
        if `level'<10 | `level'>99 {
		di in red "level() invalid"
		exit 198
	}
        if `quant'<5 | `quant'>94 {
                di in red "quant() invalid"
		exit 198
        }
        /* DispHdr adapted from stsum */
        DispHdr `level' `quant' `by'
        quietly {
                tempname fplus fminus
                sca `fplus' = (`quant' + 5)/100
                sca `fminus' = (`quant' - 5)/100
                local quant = `quant'/100
        }
        if `"`by'"' != `""' {
                        tempvar grp subuse
			sort `touse' `by' 
                        qui{
                            by `touse' `by': gen long `grp'=1 if _n==1 & `touse'
                            replace `grp'=sum(`grp') if `touse'
                        }
                        local ng = `grp'[_N]
			local i 1
			while `i' <= `ng' {
                                qui gen byte `subuse'=`touse' & `grp'==`i'
                                  /* DoStats wraps old function */
                                DoStats `subuse' `level' `quant' `fplus' `fminus' `"`by'"'
                                drop `subuse'
				local i = `i' + 1
			}
                        di in gr _dup(9) `"-"' `"+"' _dup(69) `"-"' 
        }
        DoStats `touse' `level' `quant' `fplus' `fminus' `""'
        if $ties > 0 {
            di _n in bl _col(1) "Brookmeyer-Crowley limits could be inaccurate" /*
            */ " in presence of tied survival" _n "times."
        }

        ret add
end


program define DispHdr /* by */
        args level quant ttl1
        di _n in gr _col(10) `"|  no. of"' _col(20) " no. of"  /*
        */ _col(30) "  `quant'" _col(39) /*
        */ `"|-------- `level'% Confidence limits --------|"'
        di in gr `"`ttl1'"' _col(10) `"|  events"' _col(22) `"ties"' /*
        */ _col(28) "percentile" _col(40) "Large Sample Appr. " /*
        */ " Brookmeyer-Crowley" _n /*
        */ _dup(9) `"-"' `"+"' _dup(69) `"-"' 
end

program define DoStats, rclass 
        args touse level quant fplus fminus by
        local d : char _dta[st_d]
        local t : char _dta[st_t] 
        preserve
        qui {
             keep if `touse'
             count if `d'==1
        }
        local N_deads = r(N)
        st_ct "" -> t pop die
        gen die_cum=sum((die>1)*die)
        local N_ties=die_cum[_N]
        global ties `N_ties'
        restore, preserve
        qui keep if `touse'
             /* Estimate Large Sample Approximation */
        tempname s gren lb ub median Var_med t_u t_l /*
           */ up lo s_u s_l br_up br_lo
        capture {
*                sts gen `s'=s
                /* bisogna rivedere per i casi estremi con stcox,estimate */
                stcox,estimate bases(`s')
                sts gen `gren'=se(s) 
                sts gen `lb'=lb(s), l(`level')
                sts gen `ub'=ub(s), l(`level')
        }
        if _rc!=0 {
                 local median .
                 local lo .
                 local up .
                 local br_lo .
                 local br_up .
                 Settitle `touse' `"`by'"'
                 local ttl `"$S_4"'
                 Displine `"`ttl'"' `N_deads' `N_ties' `median' /*
                  */ `lo' `up' `br_lo' `br_up'
                 exit
        }
        sort `t'
        /* Findptl as well as Settitle and Displine are adapted from stsum */
        Findptl `quant' `s' `gren'
        scalar `median' = $S_1
        if `median'==. {
                 local lo .
                 local up .
                 local br_lo .
                 local br_up .
                 Settitle `touse' `"`by'"'
                 local ttl `"$S_4"'
                 Displine `"`ttl'"' `N_deads' `N_ties' `median' /*
                  */ `lo' `up' `br_lo' `br_up'
                 exit
        }
        scalar `Var_med' = $S_2^2
        Findptl2 `fplus' `s'
        scalar `t_u' = $S_1
        scalar `s_u' = $S_2
        Findptl `fminus' `s' 
        scalar `t_l' = $S_1
        scalar `s_l' = $S_3
        tempname f_t Var_t up lo
        /* Hosmer & Lemeshow (2.12) */
        sca `f_t'=(`s_u'-`s_l')/(`t_l'-`t_u')
        /* (2.13) */
        sca `Var_t'=`Var_med'/`f_t'^2
        local z = invnorm(1-(1-`level'/100)/2)
        /* (2.14) */
        sca `up' = `median' + `z'*sqrt(`Var_t')
        sca `lo' = `median' - `z'*sqrt(`Var_t')
        sca `up' = round(`up',.0001)
        sca `lo' = round(`lo',.0001)

        /* Brookmeyer-Crowley limits are estimated according to p. 51
           as t_min(percentile,1] and not according to p. 154 as
           t_min[percentile,1]. 154 */
        Findptl2 `quant' `ub'
        scalar `br_up' = $S_1
        Findptl `quant' `lb' 
        scalar `br_lo' = $S_1
        /* saving results */
        ret scalar Med = `median'
        ret scalar L_up = `up'
        ret scalar L_lo = `lo'
        ret scalar BC_up = `br_up'
        ret scalar BC_lo = `br_lo'
        ret scalar f_t50 = `f_t'
        ret scalar Var_t50 = `Var_t'
        Settitle `touse' `"`by'"'
        local ttl `"$S_4"'
        Displine `"`ttl'"' `N_deads' `N_ties' `median' `lo' `up' `br_lo' `br_up'
end

program define Findptl /* percentile s greenwood-s.e.*/
        args p s gren
        tempvar j
        /* j0 is useless if `touse' keeped */
        quietly {
                gen long  `j' = _n if `s'<=`p'
                summ `j'
		local j=r(min)
	}
	global S_1 = _t[`j']
        global S_2 = `gren'[`j']
        global S_3 = `s'[`j']
end


program define Findptl2 /* percentile s */
        args p s 
	tempvar j
	quietly {
                gen long  `j' = _n if `s' >= `p'
                summ `j'
                local j=r(max)
	}
	global S_1 = _t[`j']
        global S_2 = `s'[`j']
end


program define Settitle /* touse byvars */
        args touse byvars
	if `"`byvars'"'==`""' {
                global S_4 `"total"'
		exit
	}
	quietly {
		tempvar j
		gen long `j' = _n if `touse'
		summ `j'
		local j = r(min)
	}
        mac shift
        local ty : type `1'
        if substr(`"`ty'"',1,3)==`"str"' {
                 local v = trim(substr(trim(`1'[`j']),1,8))
        }
        else {
                 local v = `1'[`j']
                 local v : label (`1') `v' 8
        }
        global S_4 `"`v'"'
end


program define Displine
        args ttl N_deads N_ties median lo up br_lo br_up
	local ttl = trim(`"`ttl'"')
	local skip = 8-length(`"`ttl'"')
	di in gr _skip(`skip') `"`ttl' | "' in ye /*
        */ %6.0f `N_deads' `" "' %6.0f `N_ties' /*
        */ _col(24) %9.0g `median' `"     "' /*
        */ %9.0g `lo' `" "' %9.0g `up' `" "' /*
        */ %8.0g `br_lo' `" "' %8.0g `br_up'

end

exit


stquant, by(hercoc)

         |  no. of  no. of     50     |-------- 95% Confidence limits --------|
hercoc   |  events   ties  percentile  Large Sample Appr.  Brookmeyer-Crowley
---------+---------------------------------------------------------------------
       1 |     92     20      150      102.3492  197.6508      106      190
       2 |    100     25      151      110.0249  191.9751      110      177
       3 |    136     46      186      143.8734  228.1266      147      224
       4 |    165     50      181      150.7872  211.2128      153      211
---------+---------------------------------------------------------------------
   total |    493    349      168      149.6521  186.3479      153      184

Brookmeyer-Crowley limits could be inaccurate in presence of tied survival
times.