*! version 1.1.0  31aug2001
program define gompertzhet2, eclass byable(recall) sort
	version 7
	if replay() {
		if `"`e(cmd)'"' != "gompertzhet2" { error 301 }
		if _by() { error 190 }
		Display `0'
		error `e(rc)'
		exit
	}

	syntax [varlist] [if] [in] [fweight pweight iweight] /*
		*/ [, noCOEF FRailty(string) SHared(varname) /*
		*/ CLuster(varname) Dead(varname numeric)/*
		*/ DEBUG FROM(string) noHEADer HR MLMethod(string) /*
		*/ noCONStant Level(integer $S_level) noLOG /*
		*/ OFFset(varname numeric) noLRtest/*
		*/ Robust SCore(string) T0(varname numeric) SKIP *]
	if _by() {
		_byoptnotallowed score() `"`score'"'
	}

	tokenize `varlist'
	local t `1'
	mac shift
	local rhs `*'

	if "`shared'"!="" {
		local mlmethod = "d0"
	}

	if "`cluster'"!="" {
		local cluopt cluster(`cluster')
		if "`shared'" != "" {
			di as err "cluster() not allowed with shared()"
			exit 198
		}
	}
	if "`from'" != "" { local iniopt init(`from') }
	if "`mlmethod'" == "" { local mlmethod = "d2" }
	if "`offset'" !="" { local offopt = "offset(`offset')" }
	mlopts options, `options'

	GetFrailty `frailty'
	local mlprog `r(prog)'
	local frtitle `r(frt)'
	local pred `r(pred)'

	if "`shared'"!="" {
		local mlprog "`mlprog'_sh2"
	}

	if "`score'" != "" {
		local n : word count `score'
		if `n' != 3 {
			di as err /*
		*/ "score() invalid:  three new variable names required"
			exit 198
		}
		confirm new var `score'
		local scopt "score(`score')"
	}

	if "`weight'" != "" {
		tempvar wv
		qui gen double `wv' `exp'
		local w [`weight'=`wv']
	}

	tempvar touse
	mark `touse' `w' `if' `in'
	markout `touse' `t' `rhs' `dead' `t0'
	markout `touse' `cluster' `shared', strok

	tempvar nn
	gen long `nn' = _n
	sort `touse' `shared' `nn'

	if "`dead'" != "" {
		local sdead "`dead'"
		capture assert `dead'==0 | `dead'==1 if `touse'
		if _rc {
			tempvar mydead
			qui gen byte `mydead' = `dead'!=0 if `touse'
			local dead "`mydead'"
		}
	}
	else {
		tempvar dead
		qui gen byte `dead'=1
		local sdead 1
	}
	if "`t0'" == "" {
		local t0 0
	}
	else { local topt t0(`t0') }

	capture assert `t0' < `t' if `touse'
	if _rc {
		di as err "`t0' >= `t' in some obs."
		exit 498
	}

	_rmcoll `rhs' `w' if `touse', `constant'
	local rhs "`r(varlist)'"
	global S_1

	if "`log'"!="" { local nlog="*" }

	if `"`from'"'=="" & "`lrtest'"=="" {
		if "`robust'"!="" | "`cluster'"!="" | /*
		*/ "`weight'"=="pweight" {
			local lrtest "nolrtest"
			`nlog' di as txt _n /*
			*/ "Getting starting values from gompertz model:"
		}
		else {
			`nlog' di as txt _n "Fitting comparison gompertz model:"
		}

		qui gompertz `t' `rhs' `w' if `touse', nocoef /*
		*/ dead(`dead') `topt' `offopt' `constant' `mlopts'

		tempname bp
		mat `bp' = e(b)

		if "`lrtest'"=="" {
			tempname llc
			scalar `llc' = e(ll)
		}
	}

	else if `"`from'"'=="" { /*nolrtest*/
		`nlog' di as txt _n /*
		*/ "Getting starting values from gompertz model:"
		qui gompertz `t' `rhs' `w' if `touse', nocoef /*
		*/ dead(`dead') `topt' `offopt' `constant' `mlopts'
		tempname bp
		mat `bp' = e(b)
	}

	global EREGd `dead'
	global EREGt `t'
	global EREGt0 `t0'
	if "`shared'"!="" {
		global EREG_by "`touse' `shared'"
	}

	tempvar num den
        tempname bc0

	quietly {
                if "`weight'"=="aweight" | "`weight'"=="pweight" {
                        tempvar wvn
                        summ `wv' if `touse', meanonly
                        gen double `wvn' = `wv'/r(mean)
                        local wvngen 1
                }
                else if "`weight'"!="" {
                        local wvn `wv'
                }
                else {
                        local wvn 1
                }
                gen double `num' = sum(`wvn'*`dead'*`touse')
                gen double `den' = sum(`wvn'*(`t'-$EREGt0)*`touse')
                local cons = ln(`num'[_N]/`den'[_N])
                drop `num' `den'
                if "`wvngen'"!=""  { drop `wvn' }
                mat `bc0' = ( `cons', 0 , 0 )    /* GAMMA, LNTHETA = 0 */
                mat colnames `bc0' = $EREGt:_cons gamma:_cons ln_the:_cons
        }


	if "`constant'"!="" {
		local skip = "skip"
		`nlog' di as txt _n "Fitting full model:"
	}


	if "`rhs'" != "" & "`skip'"=="" {
		`nlog' di ""
		`nlog' di as txt "Fitting constant-only model:"
		ml model `mlmethod' `mlprog' (`t': `t'=, `offopt') /*
			*/ /gamma /ln_the /*
			*/ `w' if `touse', init(`bc0') /*
			*/ missing collin nopreserve wald(0) `mlopts' /*
			*/ max search(quietly) noout `log' `options'
		local cont continue
		`nlog' di ""
		`nlog' di as txt "Fitting full model:"
	}
	else {
		 local cont wald(1)
	}

	if `"`from'"'=="" {
		tempname b0
		mat `b0' = (0)
		mat colnames `b0' = ln_the:_cons
		mat `bp' = (`bp',`b0')
		local iniopt "init(`bp')"
	}
	else { local iniopt `"init(`from')"' }

	ml model `mlmethod' `mlprog' /*
		*/ (`t': `t'=`rhs' , `offopt' `constant') /*
		*/ /gamma /ln_the `w' if `touse', `cont' noout /*
		*/ `robust' `cluopt' `scopt' `iniopt' `mlopts' /*
		*/ missing collin nopreserve /*
		*/ max search(quietly) `log' `options'

	if "`e(wtype)'" != "" {
		est local wexp `"`exp'"'
	}
	est local title "log relative-hazard form"
	est local fr_title `frtitle'
	est local predict `pred'
	est local cmd gompertzhet2
	est local t0 "`t0'"
       	est local dead `sdead'
	est local stcurve="stcurve"
	est scalar gamma = [gamma]_b[_cons]
	est scalar theta = exp([ln_the]_b[_cons])

	if `"`lrtest'"'=="" & `"`from'"'=="" {
		est scalar ll_c = `llc'
		if (_b[/ln_the] < -20) | (e(ll)<e(ll_c)) {
			est scalar chi2_c = 0
			est scalar p_c = 1
		}
		else {
			est scalar chi2_c = 2*(e(ll)-e(ll_c))
			est scalar p_c = chi2tail(1, e(chi2_c))*0.5
		}
	}

	global S_E_cmd gompertzhet2
	macro drop EREGd EREGt EREGt0 EREG_by

	Display, level(`level') `coef' `header' `hr'
	error `e(rc)'
end

program define GetFrailty, rclass
	local frailty `"`1'"'
	if "`frailty'"=="" {
		di as err "must specify -frailty(gamma | invgauss)-"
		exit 198
	}

	local l = length("`frailty'")
	if substr("gamma",1,max(1,`l')) == "`frailty'" {
		return local prog gomphet_glf
		return local frt "Gamma frailty"
		return local pred gomphet_gp
	}
	else if substr("invgaussian",1,max(1,`l')) == "`frailty'" {
		return local prog gomphet_ilf
		return local frt "Inverse-Gaussian frailty"
		return local pred gomphet_ip
	}
	else {
		di as err "unknown distribution frailty(`frailty')"
		exit 198
	}
end

program define Display
	syntax [, Level(integer $S_level) noCOEF noHEADer HR ]
	if "`coef'"=="" {
		if `"`hr'"'!=`""' {
			local hr `"eform(Haz. Ratio)"'
                }

		else {local hr }

		if "`header'" == "" {
			di _n as txt /*
			*/ "Gompertz regression, " /*
			*/ "`e(fr_title)'" " -- entry time `e(t0)'"
			ml di, `hr' level(`level') first plus
		}
		else	est di, `hr' level(`level') first plus

		_diparm gamma, level(`level') label("gamma")
		_diparm ln_the, level(`level')
		di as txt "{hline 13}{c +}{hline 64}"
		_diparm ln_the, level(`level') noprob exp label("theta")
		di as txt "{hline 13}{c BT}{hline 64}"
		if "`e(ll_c)'"!="" {
			if ((e(chi2_c) > 0.005) & (e(chi2_c)<1e4)) /*
                             */ | (e(chi2_c)==0) {
                		local fmt "%8.2f"
			}
			else local fmt "%8.2e"
			di as txt "Likelihood ratio test of theta=0: " /*
	        */ as txt "{help j_chibar:chibar2(01) = }" as res `fmt' /*
        	*/ e(chi2_c) as txt " Prob>=chibar2 = " as res %5.3f /*
        	*/ e(p_c)
		}
	}
end
exit