*! 1.02 Ariel Linden 05August2017 // fixed touse statements, made btime() required
*! 1.01 Chuck Huber 26July2017 // fixed mean table output
*! 1.00 Ariel Linden 06June2017

program define stbrier, rclass

version 13

	syntax [varlist(default=none fv)] [if] [in] ,  	///
		BTime(numlist) 								/// desired timepoint for the estimation
		[ Distribution(string)						/// distribution for -streg-
		COMPete(string)								/// competing risk for -stcreg-
		Ipcw(varlist fv)							/// varlist for IPCW
		Gen(string)  								/// generate brier score
		*											/// model options
		]

********************************************* points of interest *****************************************
// 1. St_i is the adjusted hazard (risk) probability at btime (can be estimated using any model)
// 2. Gt_i is the adjusted censoring probability at time _t. (always estimated using Cox for consistency)
// 3. Gt is the pooled censoring probability at btime. (always estimated using Cox for consistency)
**********************************************************************************************************

	// extract failure and failure value from stset
	local event "`_dta[st_bd]'"
	local event_val "`_dta[st_ev]'"

	marksample touse
	quietly {
	count if `touse'
	if r(N) == 0 error 2000
	local N = r(N)
	replace `touse' = -`touse'
	} // end quietly

	// Check on model options
    if ("`distribution'" !="") & ("`compete'" != "") {
		 di as err "Only one model type can be specified"
    exit 198
    }

	tempvar Sh t_old xbi hazi hazti St_i xbC hazC Gt1 lasthazCi Gt Gt_i

quietly {

	/**************** parametric regression for St_i ************************/

	if "`distribution'" !="" {

		// check if specified btime is in the range
		sum _t if `touse', meanonly
		if `btime' < r(min) | `btime' > r(max) {
			di as err "The specified btime is beyond the range of the data:" ///
			" check timevar and respecify btime accordingly"
			exit 498
		}

		// Generate St_i = adjusted hazard probability at btime  //
		gsort _t -`event'
		streg `varlist' if `touse', dist(`distribution') `options'
		gen `t_old' = _t
		replace _t = `btime' if `touse'
		predict `Sh', surv
		gen `St_i' = 1 - `Sh'
		replace _t = `t_old'

	}	// end St_i using streg

	/**************** Cox regression for St_i ************************/

	else if ("`distribution'" == "") & ("`compete'" == "") {

		// check if btime exists within the data
		sum _t if `touse', meanonly
		if `btime' < r(min) | `btime' > r(max) {
			di as err "The specified btime is beyond the range of the data:" ///
			" check timevar and respecify btime accordingly"
			exit 498
		}
		else if `btime' > r(min) | `btime' < r(max) {
			count if _t == `btime' & `touse'
		}
		// if btime is not in data, reset btime to previous value of _t
		if r(N) == 0 {
			sum _t if _t < `btime' & `touse'
			local btime = r(max)
		}

		// Generate St_i = adjusted hazard probability at btime  //
		gsort _t -`event'
		stcox `varlist' if `touse', estimate `options'
		predict `xbi' if `touse', xb
		predict `hazi' if `touse', basechazard
		gen `hazti' = `hazi' if _t==`btime' & `touse'
		sum `hazti' if `touse', meanonly
		replace `hazti' = r(mean) if `touse'
		gen `St_i' = 1-(exp(-`hazti')^exp(`xbi')) if `touse'

	}	// end St_i using Cox

	/**************** Competing-risk regression for St_i ************************/

	else if "`compete'" != "" {

		// check if btime exists within the data
		sum _t if `touse', meanonly
		if `btime' < r(min) | `btime' > r(max) {
			di as err "The specified btime is beyond the range of the data:" ///
			" check timevar and respecify btime accordingly"
			exit 498
		}
		else if `btime' > r(min) | `btime' < r(max) {
			count if _t == `btime' & `touse'
		}
		// if btime is not in data, reset btime to previous value of _t
		if r(N) == 0 {
			sum _t if _t < `btime' & `touse'
			local btime = r(max)
		}

		// Generate St_i = adjusted hazard probability at btime  //
		gsort _t -`event'
		stcrreg `varlist' if `touse', compete(`compete') `options'
		predict `hazi' if `touse', basecif
		gen `hazti' = `hazi' if _t==`btime' & `touse'
		sum `hazti' if `touse', meanonly
		replace `hazti' = r(mean) if `touse'

		if "`varlist'" != "" {
			predict `xbi' if `touse', xb
			gen `St_i' = 1-(1-`hazti')^exp(`xbi') if `touse'
		}
		else {
			rename `hazti' `St_i'
		}

		// get values used for competing risks
		local crvals "`e(crevent)'"
		local crvals = substr("`crvals'", strpos("`crvals'", "==") + 3, .)
		local crvals : subinstr local crvals " " ",", all

	}	// end St_i using stcrreg

	/********************** Generate IPCW ******************************/

	// Check if btime exists in the data. If not, reset btime to previous value of _t
		count if _t == `btime' & `touse'
		if r(N) == 0 {
			sum _t if _t < `btime' & `touse'
			local btime = r(max)
		}

	// Generate Gt = the adjusted censoring probability at btime //
		stset _t if `touse', failure(`event'=0)
		stcox `ipcw' if `touse', estimate efron
		predict `xbC' if `touse', xb
		predict `hazC' if `touse', basechazard
		gen `Gt1' = `hazC' if _t==`btime' & `touse'
		sum `Gt1' if `touse', meanonly
		replace `Gt1' = r(mean) if `touse'
		gen `Gt' = (exp(-`Gt1')^exp(`xbC')) if `touse'

	// Generate Gt_i = adjusted censoring probability at time _t //
		gen `lasthazCi' = `hazC'[_n-1] if _t > _t[_n-1]
		replace `lasthazCi' = `lasthazCi'[_n-1] if missing(`lasthazCi')
		gen `Gt_i' = (exp(-`lasthazCi')^exp(`xbC'))
		replace `Gt_i' = 1 in 1

	// reset failure in stset
		if "`event_val'" == "" {
			stset _t if `touse', failure(`event' = 1)
		}
		else {
			stset _t if `touse', failure(`event' = `event_val')
		}

	/********************* Generate brier scores ****************************/
		tempvar brier
		gen `brier' = 0										if `touse'
		replace `brier' = (1/`Gt_i') * ((1 - `St_i')^2)		if _t <= `btime' & _d==1 & `touse'
		replace `brier' = (1/`Gt') * ((0 - `St_i')^2)			if _t > `btime' & `touse'

		if "`compete'" != "" {
			replace `brier' = (1/`Gt_i') * ((0 - `St_i')^2)	if _t <= `btime' & inlist(`event',`crvals') & `touse'
		}


} // end quietly

	// get mean brier score and save as scalar
		quietly mean `brier' if `touse'

		// display the results
		tempname meanout n mean se ll ul level
		local `level' = r(level)
		local `n' = e(N)
		matrix `meanout' = r(table)
		matrix colnames `meanout' = "brier"
		scalar `mean' = `meanout'[1,1]
		scalar `se'   = `meanout'[2,1]
		scalar `ll'   = `meanout'[5,1]
		scalar `ul'   = `meanout'[6,1]

		disp ""
		disp "Mean estimation                      Number of obs   = "  %9.0g as result  ``n''
		disp ""
		tempname mytab
		.`mytab' = ._tab.new, col(5) lmargin(0)
		.`mytab'.width    15   |12    12    12    12
		.`mytab'.titlefmt  .     .    . %24s   .
		.`mytab'.pad       .     2     3     3    3
		.`mytab'.numfmt    . %9.0g %9.0g %9.0g %9.0g
		.`mytab'.strcolor result  .  .  . .
		.`mytab'.strfmt    %14s  .  .  . .
		.`mytab'.strcolor   text  .  .  . .
		.`mytab'.sep, top
		.`mytab'.titles " "           				/// 1
				"Mean"								/// 2
				"   Std. Err."                     	/// 3
				"[``level''% Conf. Interval]" ""  	//  4 5
		.`mytab'.sep, middle
		.`mytab'.row    "Brier score" 				///
			`mean'									///
			`se' 									///
			`ll'									///
			`ul'
		.`mytab'.sep, bottom


		matrix br  = e(b)
		return scalar brier = br[1,1]
		return matrix table = `meanout'

	// to save, or not to save, this is the question....
	quietly {
		if "`gen'" != "" {
			gen `gen' = `brier'
			label var `gen' "brier scores"
		}
	} // end gen quietly
end