*! version 1.1.0 Stephen P. Jenkins, May 2001
*! Split-population survival ('cure') model, discrete-time

program define spsurv, eclass

	version 6
	if replay() {
		if "`e(cmd)'" ~= "spsurv" {
			error 301
		}
		Replay `0'
	}
	else	Estimate `0'
end


program define Estimate, eclass

	syntax varlist  [if] [in], Id(varname) Seq(varname) /*
		*/ [CLOGlog CPR0(real -4) Level(passthru) NOCONS EForm * ]


	tokenize `varlist'
	local dead "`1'"
	mac shift
	local rhs "`*'"

			/* identify estimation subsample */

	marksample touse
	markout `touse' `id' `seq' `dead' `varlist'

			/* produce initial values */

*	quietly cloglog `dead' `rhs'  if `touse', `nocons'

	if "`cloglog'" == "" {
		quietly cloglog `dead' `rhs'  if `touse', `nocons'
	}
	else { cloglog `dead' `rhs'  if `touse', `nocons'  }

	tempname b0 cpr b1 V0
	matrix `b0' = e(b)
	matrix `V0' = e(V)

	matrix coleq `b0' = hazard

	matrix `cpr' = (`cpr0')
	matrix colnames `cpr' = cure_p:_cons
	matrix `b1' = `b0',`cpr'


	if "`nocons'" == "" {

		local ll0 = e(ll_0)
		local ll0  lf0(1 `ll0')  /* logL from constant-only model */
	}
					/* if NOCONS, then this omitted */
					

	local ll_noc = e(ll)


				/* obtain the full estimation results */


	global S_MLE_id "`id'"

	sort $S_MLE_id `touse' `seq'

	ml model d0 spsur_ll (hazard: `dead'=`rhs', `nocons') (cure_p:) /*
		*/ if `touse', maximize missing /*
		*/ `ll0' search(off) init(`b1', skip) /*
		*/ nopreserve title("Split population survival model")  /*
		*/ `options'



			/* fill in e()  */

	est scalar ll_noc = `ll_noc'

	est matrix b0  `b0'
	est matrix V0  `V0'

	est scalar cpr0 = `cpr0'

	est scalar curep = 1/(1+exp(-_b[cure_p: _cons]))
	est scalar securep = e(curep)*(1 - e(curep))*_se[cure_p: _cons]

			/* set chi2 stata for LR test =0 if c v.v. small,
			   and also then set p-value = 1 
	            	   Code adapted from -xtreg_ml-		*/

        est scalar chi2_c = cond(e(curep)<1e-5, 0, /*
                        */ abs(-2*(e(ll_noc)-e(ll))))


	estimate local cmd "spsurv"


	Replay, `level' `eform'

end

program define Replay

	syntax [, Level(int $S_level) eform ]

	if "`eform'" == "" {
		ml display, level(`level')  
	}
	else ml display, level(`level')  eform("exp(b)")
	di in gr "c = Pr(never fail) = " in ye e(curep)  _c
	di in gr "; Std.Err. = " in ye e(securep) _c 
	di in gr "; z = " in ye e(curep)/e(securep) 

			/* Boundary-value LR test */

	tempname pval
        scalar `pval' =  chiprob(1, e(chi2_c))*0.5
	if e(chi2_c)==0 { scalar `pval'= 1 }
	if ((e(chi2_c) > 0.005) & (e(chi2_c)<1e4)) | (e(chi2_c)==0) {
		local fmt "%8.2f"
       	}
        else    local fmt "%8.2e"
        di in green "Likelihood ratio test of c=0: " _c
        di in green "chibar2(01)= " in ye `fmt' e(chi2_c) _c
        di in green " Prob>=chibar2 = " in ye %5.3f `pval'
	

end