*! version 6.0.0  19apr1999
program define stgrby
	version 6
	st_is 2 analysis
	syntax varlist [if] [in] ,  [ /*
		*/ B2title(string) MEAN MEDian  P P25 P75 Rmean /*
		*/ Symbol(string) T1title(string) T2title(string) /*
		*/ XSCale(string) XLAbel YLine(string) * ]
	local summari="_t"
	if "`xlabel'"!="" { error 198 }		/* xlabel not allowed */

	if "`p'"~="" {
		di in red "must specify either p25 or p75"
		exit 198	
	}

	if "`mean'`median'`p25'`p75'`rmean'"=="" {
		local median="median"
	}

	if "`mean'"~="" & "`median'"~="" {
		di in red "mean & median may not be specified together"
		exit 198	
	}

	if "`median'"~="" {
		local sumopts
		local sumres "r(p50)"   /* where result found after -stsum- */
		local ttlhead "Median survival time"
		local median 1
		local myflag 1
	}
	else if "`p25'"~="" {
		local sumopts
		local sumres "r(p25)"   /* where result found after -stsum- */
		local ttlhead "25th percentile of survival time"
		local median 1
		local myflag 1
	}
	else if "`p75'"~="" {
		local sumopts
		local sumres "r(p75)"   /* where result found after -stsum- */
		local ttlhead "75th percentile of survival time"
		local median 1
		local myflag 1
	}
	else if "`rmean'"~="" {
		local sumopts
		local sumres "mean"  /* where result found after -stsum- */
		local ttlhead "Restricted mean survival time"
		local median 0
		local myflag 0
	}
	else {
		local sumopts
		local sumres "mean"  /* where result found after -stsum- */
		local ttlhead "Exponetally extended mean survival time"
		local median 0
		local myflag 0
	}
	
	tempvar touse x mean mean2 lbl grp strv id mine hi lo sim
	mark `touse' `if' `in'
	qui replace `touse'=0 if _st~=1
	preserve
	qui drop if `touse'==0

	quietly { 
		gen `x' = . 
		gen str1 `sim' ="" 
		gen `mean' = . 
		gen `mean2'= .
		gen str16 `lbl' = ""
		gen long `id'=_n
		tokenize `varlist'
		local ttl "`1'"
		local o 1
		local i 1
		while "``i''"!="" {
			if `i'!=1 {
				local ttl "`ttl', ``i''"
			}
			egen `grp'=group(``i'') if `touse'
			local t : type ``i''
			if substr("`t'",1,3)=="str" {
				gen str8 `strv' = ``i''
			}
			else {
				local t : value label ``i''
				if "`t'"=="" {
					gen str10 `strv'=string(``i'')
				}
				else	decode ``i'', gen(`strv') maxlength(8)
			}
			if `median'== 1 {
				summ `grp'
				local maxg=r(max)
				local j 1
				while `j'<=`maxg' {
					if `o'+1>_N { error 1001 }
					stsum if `grp'==`j'
					if `sumres'== . {
						sum _t if `grp'== `j'
						replace `mean' = /*
						*/ r(max) if `id'==`o'
						replace `mean2'= /*
						*/ r(max) if `id'==`o'+1
						replace `sim'= "(*)" if `id'==`o'
					}
					else {
						replace `mean' = /*
						*/ `sumres' if `id'==`o'
						replace `mean2'= /*
						*/ `sumres' if `id'==`o'+1
						replace `sim'= "O" if `id'==`o'
					}
					replace `x'=`i' if `id'==`o' | /*
					*/ `id'==`o'+1
					gen `mine'=1 if `grp'==`j'
					sort `mine'
					replace `lbl'=`strv'[1] if `id'==`o'+1
					drop `mine' 
					local o=`o'+2
					local j=`j'+1
				}
				drop `strv' `grp'
				capture label drop `strv'
			}
                        else {                 /* mean calculation */

				summ `grp'
				local maxg=r(max)
                                local j 1
                                while `j'<=`maxg' {
                                        if `o'+1>_N { error 1001 }
					GETMean `grp' `j'
					if "`rmean'"~="" {
						local area=`s(a1)'
					}
					else {
						local area=`s(a1)' + `s(a2)'
					}
					if `s(a2)'~=0 {
						local myflag 1
					}
					replace `mean' =  `area' if `id'==`o'
					replace `mean2'=  `area' if `id'==`o'+1
					replace `sim'= "O" if `id'==`o'
					replace `x'=`i' if `id'==`o' | /*
                                        */ `id'==`o'+1
                                        gen `mine'=1 if `grp'==`j'
                                        sort `mine'
                                        replace `lbl'=`strv'[1] if `id'==`o'+1
                                        drop `mine'
                                        local o=`o'+2
                                        local j=`j'+1
                                }
                                drop `strv' `grp'
                                capture label drop `strv'
                        }
                        local i=`i'+1
		}	
		if `median'==1 {
			stsum if `touse'==1 
			local yli=`sumres'
		}
		else {
			GETGMean
			if "`rmean'"~="" {
				local yli=`s(a1)'
			}
			else {
				local yli=`s(a1)' + `s(a2)'
			}
		}
		local adj=.18*(`i'-1)/3
		replace `x'=`x'-`adj' if `lbl'!=""
		sort `x' `mean'
		by `x': gen `lo' = `mean'[1] if `x'!=.
		replace `mean'=-`mean'
		sort `x' `mean'
		replace `mean'=-`mean'
		by `x': gen `hi' = `mean'[1] if `x'!=.
	}
	if `myflag'==0 {
		local ttlhead /*
		*/ "Mean survival time"
	}
	if `"`t1title'"'=="" & `"`t2title'"'=="" {
		local t1 " "
		local t2 `"`ttlhead'"' 
	}
	else {
		if `"`t2title'"'=="" {
			local t1 " "
			local t2 `"`t1title'"'
		}
		else {
			local t1 `"`t1title'"'
			local t2 `"`t2title'"'
		}
	}
	if `"`b2title'"'=="" {
		local b2title `"`ttl'"'
	}

	if "`yline'"=="" {
		local yline "`yli'"
		local opyline="yline(`yline')"
	}
	if "`yline'"=="." {
		local opyline=""
	}

	if "`xscale'"=="" {
	local losca=`i'-.5
		local xscale ".5,`losca'"
	}

	if "`symbol'"=="" {
		local symbol "s([`sim']..[`lbl'])"
	}
	else {
		local symbol = "s(" + substr("`symbol'",1,1)+"..[`lbl'])"
	}
	graph `mean' `hi' `lo' `mean2' `x', `symbol' c(.||.) trim(10) /*
		*/ border xsca(`xscale') t2(`"`t2'"') t1(`"`t1'"') /*
		*/ b2(`"`b2title'"') `opyline' `options'
end

program def GETGMean, sclass   /*Grand Mean */
	tempvar s
	sts gen `s'=s
	tempvar prior gap area sarea ns
	gsort -`s' _t
	gen double `prior'=_t0 in 1
	gen double `ns'=1 in 1
	replace `prior'=_t[_n-1] if `prior'==.
	replace `ns'=`s'[_n-1] if `ns'==.
	gen `gap'=_t-`prior'
	gen `area'=`gap'*`ns'
	gen `sarea'=sum( `area')
	sum `sarea'
	local area1=r(max)
	sreturn local a1=r(max)
	sum _t
	local mt=r(max)
	tempvar ms
	gen `ms'=`s' if _t==`mt'
	sum `ms'
	tempname msurv lamb area2
	scalar `msurv'=r(max) 	
	scalar `lamb'=-log(`msurv')/`mt'
	scalar `area2'=(1/`lamb')*exp(-`lamb'*`mt')
	if `lamb'==. {
		local area2 0
	}
	sreturn local a2=`area2'
*noi di in red `lamb' "  " `mt' "  " `area2'+`area1'
end

program def GETMean, sclass
	args grp j
	tempvar s
	sts gen `s'=s if `grp'==`j'
	tempvar prior gap area sarea ns
	gsort -`s' _t
	gen double `prior'=_t0 in 1
	gen double `ns'=1 in 1
	replace `prior'=_t[_n-1] if `prior'==.
	replace `ns'=`s'[_n-1] if `ns'==.
	gen `gap'=_t-`prior'
	gen `area'=`gap'*`ns'
	gen `sarea'=sum(`area')
	sum `sarea'
        sreturn local a1=r(max)
        sum _t if `grp'==`j'
        local mt=r(max)
        tempvar ms
        gen `ms'=`s' if _t==`mt'
        sum `ms' if `grp'==`j'
	tempname msurv lamb area2
        scalar `msurv'=r(max)
        scalar `lamb'=-log(`msurv')/`mt'
        scalar `area2'=(1/`lamb')*exp(-`lamb'*`mt')
	if `lamb'==. {
		local area2 0
	}
        sreturn local a2=`area2'

end