*! version 1.0.0  24oct2007
program survci, sortpreserve
	version 10.0
	if "`e(cmd2)'" != "stcox" {
		error 301
	}
	syntax [if] [in], [ BASESurv(varname)  AT(numlist) CLtype(string) ///
			    Level(cilevel) OUTfile(string) SURVONLY HL 	  ///
			    LEGend(passthru) * ]
	if `"`basesurv'"' == "" {
		if "`e(bases)'" == "" {
			di as err "basesurv() must be specified with stcox"
			exit 198
		}
		local basesurv `e(bases)'
		di as txt `"Note: variable `e(bases)' is used as a baseline survivor function"'
	}
	if `"`hl'"' != "" {
		local cltype loglog
	}
	marksample touse
	qui replace `touse' = 0 if !e(sample)
	local xvars: colnames e(b)
	local nxvars: word count `xvars'
	if `"`at'"' == "" {
		tokenize `xvars'
		forvalues i=1/`nxvars' {
			qui summ ``i'' if `touse', meanonly
			if r(min)==0 & r(max)==1 {
				local at `at' 0
			}
			else {
				local at `at' `=r(mean)'
			}				
		}
	}
	else {
		if `: word count `at''<`nxvars' {
			di as err `"at() must contain `nxvars' values or nothing"'
			exit 198
		} 
	}
	local holdvar `xvars'
	tokenize `at'
	forvalues i=1/`nxvars' {
		gettoken var holdvar: holdvar	
		local ex0b `ex0b' `plus' _b[`var']*``i''
		local plus +
		local note `note'`semic' `var'=``i''
		local semic ;
		tempvar xbar`i' c`i'
		local strvars `strvars' `xbar`i''
	}
	if `nxvars'!=0 {
		local ex0b = exp(`ex0b')
	}
	else {
		local ex0b=1
	}
	tempvar exb lns0 time chat nfail a ahat s se lb ub
	qui {
		gen double `lns0' = ln(`basesurv') if `touse'
		predict double `exb' if `touse', hr
		gen `time' = -_t
		sort `time' _d `touse'
		gen double `chat' = sum(`exb') if `touse'
		tokenize `xvars'
		forvalues i = 1/`nxvars' {
			gen double `c`i'' = sum(``i''*`exb') if `touse'
		}
		by `time' _d: gen long `nfail' = _N if _n==_N & _d==1 & `touse'
		sort _t `nfail' `touse'
		tokenize `at'
		forvalues i = 1/`nxvars' {
			if `"`hl'"' != "" {
gen double `xbar`i'' = sum(`nfail'*(`c`i''/`chat')/(`chat')) ///
							if `touse' & `nfail'<.
replace `xbar`i'' = `xbar`i''/abs(`lns0') - ``i'' if `touse' & `nfail'<.
			}
			else {
gen double `xbar`i'' = sum(`nfail'*(`c`i''/`chat'-``i'')/(`chat')) ///
							if `touse' & `nfail'<.
			}
		}
		gen double `ahat' = sum(`nfail'/(`chat'*(`chat'))) ///
							if `touse' & `nfail'<.
		if `"`hl'"' != "" {
			replace `ahat' = `ahat'/(`lns0')^2 ///
							if `touse' & `nfail'<.
		}
		tempvar touse2
		gen byte `touse2' = `touse'*(`nfail'<.)
		mata: _workSE(`"`strvars'"')
		if `"`hl'"' != "" {
			gen double `se' = sqrt(`ahat'+`bhat') if `touse'
		}
		else {
			gen double `se' = `ex0b'*sqrt(`ahat'+`bhat') if `touse'
		}
		// must have only one nonmis obs per failure time
		by _t: replace `se' = sum(`se') if `touse' & _d==1
		replace `se' = `se'[_n-1] if `touse' & `se'>=. & _d==0
		gen double `s' = `basesurv'^`ex0b'
		local z = invnormal(1-(1-`level'/100)/2)
		if `"`cltype'"'=="" | `"`cltype'"'=="loglog" {
			local cltype loglog
			if `"`hl'"' == "" {
				gen double `lb'=`s'^exp(`z'*`se'/abs(ln(`s')))
				gen double `ub'=`s'^exp(-`z'*`se'/abs(ln(`s')))
			}
			else {
				gen double `lb' = `s'^exp(`z'*`se')
				gen double `ub' = `s'^exp(-`z'*`se')
			}
			local cinote loglog based
		}
		else if `"`cltype'"' == "normal" {
			gen double `lb' = `s'-`z'*`se'*`s'
			gen double `ub' = `s'+`z'*`se'*`s'
			local cinote normal based
		}
		else if `"`cltype'"' == "log" {
			gen double `lb' = `s'*exp(-`z'*`se')
			gen double `ub' = `s'*exp(`z'*`se')
			local cinote log based
		}
		else {
			di as err `"cltype(): `cltype' is not allowed"'
			exit 198
		}
		if `"`hl'"' != "" {
			replace `se' = `se'*abs(ln(`s'))
		}
		replace `se' = `se'*`s'
	}
	if `"`outfile'"' != "" {
		_prefix_saving `outfile'
		local saving `"`s(filename)'"'
		local replace `"`s(replace)'"'
		preserve
		qui keep if `touse'
		keep _t _d `s' `se' `lb' `ub'
		order _t _d `s' `se' `lb' `ub'
		rename `s' _surv
		label variable _surv "S_hat(t,x0,b)"
		rename `lb' _lb
		label variable _lb "`level'% lower bounds (pointwise, `cltype' based)"
		rename `ub' _ub
		label variable _ub "`level'% upper bounds (pointwise, `cltype' based)"
		rename `se' _se
		label variable _se "SE of S_hat(t,x0,b)"
		sort _t _d
		label data "Survival estimates adjusted for `note'"
		compress
		save `saving', `replace'
		restore
	}
	local title Covariate-adjusted survival estimate
	local xtitle analysis time
	if `"`legend'"' == "" {
		local legend legend(label(1 "`level'% CI (`cinote')") label(2 "Survivor function"))
	}
	if `"`survonly'"'=="" {
		local cigraph `cigraph' (rline `lb' `ub' _t if `touse', ///
								sort connect(J))
	}
	local plots `cigraph' (line `s' _t if `touse', sort connect(J))
	twoway `plots' , title(`title') xtitle(`xtitle') `legend' ///
				note(Covariates: `note') `options'
end

version 10
mata:
mata set matastrict on

void _workSE(string scalar vars)
{
	real rowvector indvar, touse2
	real colvector r, Bhat
	real matrix Xc, V
	string scalar bhat

	indvar = st_varindex(tokens(vars))
	st_view(touse2, ., st_local("touse2"))
	r = select(range(1,rows(touse2),1), touse2)
	st_view(Xc, ., indvar, st_local("touse2"))
	V = st_matrix("e(V)")
	if (V == J(0,0,.)) Bhat = J(rows(Xc),1,0) /* no covariates */
	else Bhat = rowsum(Xc*V:*Xc)
	(void) st_addvar("double", bhat = st_tempname())
	st_store(r,bhat, Bhat)
	st_local("bhat", bhat)
}
end