*! version 1.0.1 21Nov2011

/*
History
MJC 21Nov2011 version 1.0.1 - added nolog and keepcons options. Additional user defined constraints can now be used.
MJC 07Oct2011 version 1.0.0
*/

program stmix, eclass byable(onecall)
	version 11.2
	if _by() {
		local by "by `_byvars'`_byrc0':"
	}
	if replay() {
		syntax  , [DISTribution(string) *]
		if "`distribution'" != "" {
			`by' Estimate `0'
		}
		else {
			if "`e(cmd)'" != "stmix" {
				error 301
			}
	                if _by() {
        	                error 190
	                }
			Replay `0'
		}
		exit
	}
	`by' Estimate `0'
end


program Estimate, eclass byable(recall)
	st_is 2 analysis
		syntax [varlist(default=empty)] [if] [in] 	[,										///
														DISTribution(string)				///
														LAMBDA1(varlist)					///
														GAMMA1(varlist)						///
														LAMBDA2(varlist)					///
														GAMMA2(varlist)						///
														PMIX(varlist)						///
																							///
														NOHR								///
														SHOWCons							///
														KEEPCons							///
														SHOWINIT							///
																							///
														PMIXConstraint(real 0)				///
														Level(cilevel)						///
														NOLOG								///
														*									///
													]

																	
		/*  Weights */
        if "`weight'" != "" {
			display as err "weights not allowed"
			exit 198
        }
        local wt: char _dta[st_w]       
        if "`wt'" != "" {
			display as err "weights not allowed"
			exit 198
		}
																	
		marksample touse
		markout `touse' `lambda1' `gamma1' `lambda2' `gamma2' `pmix'
		qui replace `touse' = 0 if _st==0
		
		mlopts mlopts , `options'	/* Extract any ml options to pass to ml model */
		local extra_constraints `s(constraints)'

		/*********************************************************************************************************************************************/
		/* Preliminaries */
		
			/* Factor variables not allowed for tvc varables */
			fvexpand `varlist' `lambda1' `gamma1' `lambda2' `gamma2' `pmix'
			if "`r(fvops)'" != "" {
				display as error "Factor variables not allowed. Create your own dummy varibles."
				exit 198
			}

			if "`distribution'"=="" {
				local dist "ww"
			}
			else {
				local l = length("`distribution'")
				if substr("weibexp",1,max(5,`l')) == "`distribution'" | "we" == "`distribution'" {
					local dist "we"
				}
				else if substr("weibweib",1,max(5,`l')) == "`distribution'" | "ww" == "`distribution'" {
					local dist "ww"
				}
				else {
					di as error "Unknown distribution"
					exit 198
				}
			}			
			
			if "`dist'"=="ww" {
				local geqn "(ln_gamma2: =`gamma2')"
				global gg "gg"
				local text "Weibull-Weibull"
			}
			else {
				global gg
				local text "Weibull-exponential"
			}
		
			local showin "quietly"
			if "`showinit'"!="" {
				local showin
			}
					
			if "`dist'"=="we" & "`gamma2'"!="" {
				di as error "gamma2 invalid when distribution is weibexp"
				exit 198
			}
			
		/* Check Time Origin */
		qui summ _t0 if `touse', meanonly
		if r(max)>0 {
			display in green  "note: delayed entry models are being fitted"
			global del_entry = 1
		}
		else {
			global del_entry = 0
		}
		
		if wordcount("`lambda1' `gamma1' `lambda2' `gamma2' `pmix'")>0 {
			global nhr "no"
		}
		else {
			global nhr "yes"
		}
		
		/*********************************************************************************************************************************************/
		/* Obtain initial values */
		
			constraint free
			constraint `r(free)' [logit_p_mix][_cons] = `pmixconstraint'
			local cons "`r(free)'"
			local conslist "`r(free)'"
			constraint free
			constraint `r(free)' [xb][_cons] = 0
			local conslist "`conslist' `r(free)'"
			local conslist2 "`r(free)'"
			
			local dropconslist `conslist'
			/* If further constraints are listed stpm2 then remove this from mlopts and add to conslist */
				if "`extra_constraints'" != "" {
					local mlopts : subinstr local mlopts "constraints(`extra_constraints')" "",word
					local conslist `conslist' `extra_constraints'
					local conslist2 `conslist2' `extra_constraints'
				}
				
			local constopts "constraints(`conslist')"
			di in green "Obtaining initial values:"
			di ""
			`showin' ml model lf stmix_lf_`dist'							///
							(xb: _t _d _t0= `varlist')						///
							(logit_p_mix: =`pmix')	 						///
							(ln_lambda1: =`lambda1')						/// 
							(ln_gamma1: =`gamma1') 							///
							(ln_lambda2: =`lambda2')						///
							`geqn'	 										///
							if `touse',										///
							`constopts'										///
							`mlopts'										///
							`nolog'											///
							iters(100)										///
							maximize
			
			`showin' ml display
			
			/* Extract initial values for final estimation */
			mat initmat = e(b)
			
		/*********************************************************************************************************************************************/
		/* Final estimation */

			local constopts "constraints(`conslist2')"
			di in green "Fitting full model:"
	
			ml model lf stmix_lf_`dist'										///
							(xb: _t _d _t0= `varlist')						///
							(logit_p_mix: =`pmix')	 						///
							(ln_lambda1: =`lambda1')						/// 
							(ln_gamma1: =`gamma1') 							///
							(ln_lambda2: =`lambda2')						///
							`geqn'	 										///
							if `touse',										///
							`constopts'										///
							`mlopts'										///
							search(off)										///
							waldtest(0)										///
							init(initmat, skip)								///
							`nolog'											///
							maximize

			ereturn local title "Mixture `text' proportional hazards regression"
			ereturn local cmd stmix
			ereturn local predict stmix_pred
			ereturn local distribution "`dist'"
			ereturn local varlist "`varlist'"
			ereturn local lambda1varlist "`lambda1'"
			ereturn local gamma1varlist "`gamma1'"
			ereturn local lambda2varlist "`lambda2'"
			ereturn local gamma2varlist "`gamma2'"
			ereturn local pmixvarlist "`pmix'"
			
			
			if "`keepcons'" == "" {
				constraint drop `dropconslist'
			}
			else {
				ereturn local sp_constraints `dropconslist'
			}
			
			Replay, level(`level') `showcons' `nohr' `gg'
		
end

program Replay
		syntax [, Level(cilevel) SHOWCons NOHR GG]
			
			if "`showcons'"!="" {
				local constdi
			}
			else {
				local constdi "nocnsreport"
			}
			
			if "`nohr'"=="" & "$nhr"=="yes"{
				local hr "hr"
			}
			
			ml display, `constdi' /*first plus*/ `hr' level(`level')
			
		/*	_diparm ln_lambda1, exp label("lambda_1")
			_diparm ln_gamma1, 	exp label("gamma_1")
			_diparm __sep__
			_diparm ln_lambda2, exp label("lambda_2")
			if "$gg"!="" {
				_diparm ln_gamma2, 	exp label("gamma_2")
			}
			_diparm __sep__
			_diparm logit_p_mix, invlogit label("p_mix")
			_diparm __bot__		*/
end