*! v.1.0.0 - 23sept2014
*! N.Orsini, A.Discacciati

capture program drop stqkm
program stqkm, eclass
version 10

local cmdline : copy local 0

syntax [varlist(default=none)] [if] [in]  [ ,  ///
   Reps(integer -1)   ///
   MINReps(integer -1) /// undocumented
   Level(integer $S_level) ///
   SEED(string) ///
   Quantiles(numlist) ///
   WGTv(varlist min=1 max=1) /// undocumented
   NODOTS ///
   ]

	marksample touse, strok novarlist

	if "`log'"!="" | "`dots'"!="" {
			local log "*"
	}

	// Check if any covariate. If none

		if "`varlist'" == "" local overall 1
		else local overall 0

	// Check data st

		st_is 2 analysis

	// Check if data is already weighted

		if "`_dta[st_wv]'" != "" {
			di in red "weighted data not allowed"
			exit 198
		}

	// Check bootstrap

		if `reps' != -1 & `minreps' != -1 {
			di in red "reps() and minreps() can not be used at the same time"
			exit 198
		}

		if `reps' == -1 & `minreps' == -1 local reps 20

		if (`reps' < 2 & `minreps' == -1) | (`reps' == -1 & `minreps' < 2) {
			di in red "insufficient bootstrap replications"
			exit 198
		}

	// Check number of observations

		quietly count if `touse'
		if r(N) < 4 {
			di in red "insufficient observations"
			exit 2001
		}

	// Get exposure variable name

		local firstname : word 1 of `varlist'
	    local xname "`varlist'"
		local ixlist "`varlist'"


		if regexm("`firstname'","_I") == 1 {
				local varlab: variable label `firstname'
				local poseq =  strpos("`varlab'","=")
				local xname = substr("`varlab'",1,`poseq'-1)
				local ixlist "`varlist'"
		}
		else {
			if `: word count `varlist'' > 1 {
					di as err "specify just one variable name"
					exit 198
			}
			if `: word count `varlist'' == 1 {
				qui tab `varlist'
				if r(r) > 2 {
						di as err "`varlist' has more than two levels, use the xi: notation"
						exit 198
				}
				local xname "`varlist'"
			}

		}

		if `overall' == 1  {
						tempvar const1
						gen `const1' =  1
						local xname "`const1'"
						local ixlist "`xname'"
		}

	// Recode exposure as 1,2,3...

		tempvar expr
		qui gen `expr' = .
		local j = 0
		qui levelsof `xname', local(exposurel)
		capture confirm numeric variable `xname'
		if !_rc {
			foreach name of local exposurel  {
				qui replace `expr' = `j' if `xname' == `name'
				local j = `j' + 1
			}
		}
		else {
			foreach name of local exposurel  {
				qui replace `expr' = `j' if `xname' == "`name'"
				local j = `j' + 1
			}
		}

	// Identify the reference level

		qui su `xname'
		local bl = r(min)

	// More than 2 levels

		if "`ixlist'" != "" {
			if "`: char `xname'[omit]'" != "" local bl : char `xname'[omit]
		}

	// Check specified quantiles

		tempname quantlist
		SetQ `quantiles'
		local listquant "`r(quants)'"
		tokenize "`r(quants)'"
		local nq 1
		while "``nq''" != "" {
			local q`nq' ``nq''
			local nq = `nq' + 1
		}
		local nq = `nq' - 1

	// Check quantiles and get a list

		local pow = 10
		local done = 0
		while !`done' {
			local pow = `pow'*10
			local done = 1
			forvalues k=1/`nq' {
				local q = (`q`k'')*`pow'
				local done = `done' & /*
					*/ (reldif(`q',trunc(`q'))<c(epsfloat))
			}
		}
		local f = trunc(log10(`pow'))
		forvalues k=1/`nq' {
			local myeq = "q" + /*
				*/ string(trunc((`q`k'')*`pow'),"%0`f'.0f")
			local eqnames `eqnames' `myeq'
		}
		local eqnames : list uniq eqnames
		local k : word count `eqnames'
		if `k' < `nq' {
			di as err "only `k' out of `nq' quantiles are " /*
			 */ "unique within a relative precision of c(epsfloat)"
			exit 498
		}

	// Define labels for the regression coefficients for each quantile

		tokenize  "`varlist'"

		local eclist ""
		local k 1
		while `k' <= `nq' {
			local myeq : word `k' of `eqnames'
			local eclist "`eclist' `myeq'"
			local i 1
			while "``i''" != "" {
				local eqnams "`eqnams' `myeq'"
				local conams "`conams' ``i''"
				tempvar v
				local vl "`vl' `v'"
				local i = `i' + 1
			}
 			local eqnams "`eqnams' `myeq'"
			local conams "`conams' _cons"

			tempvar v
			local vl "`vl' `v'"
			local k = `k' + 1
		}


	tempname coefb

	local rs : word count `eqnams'
	mat `coefb' = J(1, `rs', .)

	forv s = 1/`rs' {
		local result "`result' `coefb'[1,`s']"
	}

	// Get a clean list of percentiles (25, 50, ...) rather than quantiles (.25, .50, ...)

	local listquant2 : subinstr local eqnames "q" "", all

	// Get name of surival time and failure

	local failure : char _dta[st_bd]
	local depname : char _dta[st_bt]
    local idname : char _dta[st_id]
	preserve

	quietly keep if `touse'
	quietly keep if _st != 0
	quietly keep _* `xname'  `wgtv' `depname'  `failure' `expr' `idname'

	local nobs = `c(N)'
	local tdf = `c(N)' - (`:word count `varlist'' + 1)

	tempfile BOOTRES BOOTJ
	tempname coefs VCE  handle

	qui save `BOOTJ', replace // To be used in the bootstrap loop

	// Point estimates on original sample
	stqkm_f , expvar(`expr' `xname') quantilesf(`listquant2') wgtvv(`wgtv') base(`bl')  nexp(`ixlist') boot(0)
 	mat `coefs' = e(_beta)
	mat colnames `coefs' = `conams'
	mat coleq `coefs' = `eqnams'

	// Confidence intervals on bootstrap samples

	*`log' di in gr "(bootstrapping " _c // Dots by Nicola
	 qui postfile `handle' `vl' using "`BOOTRES'", double

	// set the seed
	if "`seed'" != "" {
		set seed `seed'
	}
	local seed `c(seed)'

	local j 1

	// Reps
	if `reps' != -1 {

	di ""
	if "`nodots'" == "" nois _dots 0, title(Bootstrap replications (`reps')) // reps(`reps')

		while `j'<=`reps' {

			qui use `BOOTJ', clear // Necessary to sample every time from the original dataset

			bsample

			local listquj : subinstr local eqnames "q" "", all

			stqkm_f , expvar(`expr' `xname') quantilesf(`listquj') wgtvv(`wgtv') base(`bl') nexp(`ixlist')
			mat `coefb' = e(_beta)
			version 6: post `handle' `result'
			*if e(mp) == 0 `log' di in gr "." _c // Dots by Nicola
			*else `log' di in red "x" _c // Dots by Nicola

			if "`nodots'" == "" {
				if e(mp) == 0 nois _dots `j' 0
				else nois _dots `j' 1
			}

			local j = `j' + 1

		}

	if "`nodots'" == "" di _n
	}

	// MINReps
	local nsuccess 0

	if `minreps' != -1 {

	di ""
	nois _dots 0, title(Bootstrapping...)

		while `nsuccess' < `minreps' {

			qui use `BOOTJ', clear

			bsample

			local listquj : subinstr local eqnames "q" "", all

			stqkm_f , expvar(`expr' `xname') quantilesf(`listquj') wgtvv(`wgtv') base(`bl') nexp(`ixlist')
			mat `coefb' = e(_beta)
			version 6: post `handle' `result'

			if e(mp) == 0 local fail 0
			else local fail 1

			local nsuccess = `nsuccess' + (`fail' == 0)

			nois _dots `j' `fail'

			local j = `j' + 1
		}

	local reps = `j' - 1

	}

	postclose `handle'

	qui use "`BOOTRES'", clear
	*qui bstat `vl' , stat(`coefs') n(`nobs')
	*estat bootstrap,  all
	*mat list e(ci_normal)
	*mat list e(ci_percentile)

	tempvar missingpercentiles
	qui egen `missingpercentiles' = rowmiss(`vl')
	qui count if `missingpercentiles' != 0
	local missperc = `r(N)'
	local reps = `reps' - `missperc'
	quietly mat accum `VCE' = `vl', dev nocons
	mat rownames `VCE' = `conams'
	mat roweq `VCE' = `eqnams'
	mat colnames `VCE' = `conams'
	mat coleq `VCE' = `eqnams'
	mat `VCE'=`VCE'*(1/(`reps'-1))
	version 6: est post `coefs' `VCE', obs(`c(N)')  /*dof(`tdf')*/  depn(_t)

	*`log' noi di in gr ")" // Dots by Nicola
	restore
	version 6: est repost, esample(`touse')
	capture erase "`BOOTRES'"
	version 11

	ereturn scalar reps = `reps'
	ereturn scalar N = `nobs'

	// Display header

	if `nq' == 1 di in gr "Quantiles from Kaplan-Meier" _col(54) _c
	else di in gr "Simultaneous Quantiles from Kaplan-Meier" _col(54) _c
	di in gr "Number of obs =" in ye %10.0g e(N)


	// Save results

	ereturn scalar reps = `reps'
	ereturn scalar N = `nobs'
	ereturn scalar n_q = `nq'
	ereturn local eqnames "`eqnames'"
	ereturn local vce  "bootstrap"
	ereturn local vcetype "Bootstrap"
	ereturn local predict "sqreg_p"
	ereturn local cmdline `"stqkm `cmdline'"'
	ereturn local cmd "stqkm"

	// Display results

	ereturn display, level(`level')

	// Warning message

	if `missperc' != 0 {
	di as text "Note: one or more quantile estimate could not be estimated in `missperc' bootstrap"
	di as text "replicates; standard errors estimates include only complete replications."
	}

	// Restore original stset
	qui streset [pw = 1]

end



*******************
**				 **
**	SUBPROGRAMS  **
**				 **
*******************

capture program drop stqkm_f
program stqkm_f, eclass
version 11
syntax , expvar(string) quantilesf(numlist) base(string) nexp(string) [ wgtvv(string) boot(string) ]

tokenize `expvar'
tempvar exposure origexp
gen `exposure' = `1'
gen `origexp'  = `2'

	qui tab `exposure'
	local levx = r(r)


	if "`wgtvv'" != "" {
		qui streset  [pw = `wgtvv']
	}

	* Get the percentiles of the Kaplan-Meier

	tempvar atrisk
	qui gen double `atrisk' = _t - _t0

	tempvar flag
	qui bysort `exposure':  gen `flag' = _n == 1

	tempvar grp subuse
	qui bysort  `exposure': gen long `grp'=1 if _n==1
	qui replace `grp'=sum(`grp')
	local ng = `grp'[_N]

    foreach k of local quantilesf {
			tempvar out`k'
			qui gen `out`k'' = .
			label var `out`k'' "q`k'"
			char `out`k'' [varname] "q`k'"
			local listvarquant "`listvarquant' `out`k''"
	}


	local i 1
	while `i' <= `ng' {
			    capture drop `subuse'
				qui gen byte `subuse'= `grp'==`i'

				DoStats `subuse' `atrisk' `exposure' , percentiles(`quantilesf')

				foreach k of local quantilesf {

					if "`boot'" == "0" { // Check if all quantiles are estimable in the original dataset; Andrea 12sep2011
							if r(p`k') == . {
								di as err "Quantile `k'% does not exist in at least one of the exposure groups"
								exit 198
							}


					}

					qui replace `out`k'' = r(p`k') if `grp'==`i'

				}

				local i = `i' + 1
	}

	* list `origexp'  `listvarquant' if `flag' == 1, clean noobs subvarname

	// Save differences among percentiles as vectors

	tempname _bacc
	char `origexp'[omit] `base'
	local s = 1
	foreach v of local listvarquant {
			if `levx' > 1  qui reg `v' `nexp' // basta passare i nomi e non usare i.
			// Check missing values
			local names : colfullnames e(b)

			foreach x of local names {
				if substr("`x'",1,2)=="o." local missingp = substr("`x'",3,length("`x'")-2)
			}
			else qui reg `v'

			tempname matclean
			mat `matclean' = e(b)

			if "`missingp'" != "" mat `matclean'[1, colnumb(`matclean',"`missingp'")]  = .

			if `s' == 1 {
					mat `_bacc' = `matclean'
					}
			else  {
					mat `_bacc' = `_bacc' , `matclean'
			}

		local `++s'
	}

	if "`missingp'" != "" ereturn scalar mp = 1
	else ereturn scalar  mp   = 0
	ereturn matrix _beta = `_bacc'

end

capture program drop DoStats
program define DoStats, rclass /* touse atrisk by */
	syntax varlist  , percentiles(string)

	tempvar touse atrisk by
	tokenize `varlist'
    qui gen `touse' = `1'
	qui gen `atrisk' = `2'
	qui gen `by' = `3'

	local id : char _dta[st_id]
	local wv  : char _dta[st_wv]

	tempname tatr m s j0
	quietly {
		if `"`wv'"'==`""' {
			summ `atrisk' if `touse'
			scalar `tatr' = r(sum)
			summ _d if `touse'
			local fail = r(sum)
			if `"`id'"'==`""' {
				count if `touse'
				local nsubj = r(N)
			}
			else {
				sort `touse' `id'
				by `touse' `id': /*
				*/ gen byte `m'=1 if _n==1 & `touse'
				summ `m'
				local nsubj = r(sum)
			}
		}
		else {
			tempvar z
			gen double `z' = `wv'*`atrisk'
			summ `z' if `touse'
			scalar `tatr' = r(sum)
			replace `z' = `wv'*_d
			summ `z' if `touse'
			local fail = r(sum)
			drop `z'
			if `"`id'"'==`""' {
				summ `wv' if `touse'
				local nsubj = r(sum)
			}
			else {
				sort `touse' `id'
				by `touse' `id': /*
				*/ gen double `z'=`wv' if _n==1 & `touse'
				summ `z' if `touse'
				local nsubj = r(sum)
			}
		}

		local ttl `"$S_1"'
		capture GetS `touse' `s'
		if _rc==0 {
			replace `touse' = 0 if `s'>=.
			sort `touse' _t
			gen long `j0' = _n if `touse'==0
			summ `j0'
			local j0 = cond(r(max)>=.,1,r(max)+1)

			foreach k of local percentiles {
				  Findptl `=(100-`k')/100' `s' `j0'
				  local p`k' $S_1
			}
		}
		else {

				foreach k of local percentiles {
				local p`k' .
			}
		}
	}


    noisily foreach k of local percentiles {

	ret scalar p`k' = `p`k''
	}


end


capture program drop SetQ
program define SetQ /* <nothing> | # [,] # ... */ , rclass
	version 6.0
		if "`*'"=="" {
		ret local quants ".5"
		exit
	}
	local orig "`*'"
	tokenize "`*'", parse(" ,")

	while "`1'" != "" {
		FixNumb "`orig'" `1'
		ret local quants "`return(quants)' `r(q)'"
		mac shift
		if "`1'"=="," {
			mac shift
		}
	}
end

capture program drop FixNumb
program define FixNumb /* # */ , rclass
	version 6.0
	local orig "`1'"
	mac shift
	capture confirm number `1'
	if _rc {
		Invalid "`orig'" "`1' not a number"
	}
	if `1' >= 1 {
		ret local q = `1'/100
	}
	else 	ret local q `1'
	if `return(q)'<=0 | `return(q)'>=1 {
		Invalid "`orig'" "`return(q)' out of range"
	}
end

capture program drop Invalid
program define Invalid /* "<orig>" "<extra>" */
	version 6.0
	di in red "quantiles(`1') invalid"
	if "`2'" != "" {
		di in red "`2'"
	}
	exit 198
end

capture program drop Findptl
program define Findptl /* percentile s `j0' */
	args p s j0

	if `j0'>=. {
		global S_1 .
		exit
	}

	tempvar j
	quietly {
		gen long  `j' = _n if float(`s')<=float(`p') in `j0'/l
		summ `j' in `j0'/l
		local j=r(min)
	}
	local t : char _dta[st_t]
	global S_1 = _t[`j']
end

capture program drop GetS
program define GetS /* touse s */
	args touse s

	tempname prior
	capture {
		capture estimate hold `prior'
		stcox if `touse', estimate bases(`s')
	}
	local rc=_rc
	capture estimate unhold `prior'
	exit `rc'
end

exit


webuse catheter, clear
xi: stqkm i.female, q(25 50 75) reps(20)