*! version 1.1.0  01may2001
program define weibullhet2, eclass byable(recall) sort
	version 7
	if replay() {
		if `"`e(cmd)'"' != "weibullhet2" { 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 MLMethod(string) /*
		*/ noCONStant Level(integer $S_level) noLOG /*
		*/ noLRtest HAzard/*
		*/ Robust SCore(string) T0(varname numeric) SKIP TR HR *]
	if _by() {
		_byoptnotallowed score() `"`score'"'
	}

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

	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 "`shared'" != "" {
		local mlmethod = "d0"
	}
	mlopts options, `options'
		
	GetFrailty `frailty' `hazard'
	local mlprog `r(prog)'
	local frtitle `r(frt)'
	local pred `r(pred)'

	if "`shared'" != "" {
		local mlprog "`mlprog'_sh2"
	}
	if "`hazard'"!="" {
                if "`tr'"!="" {
                        di as err "tr invalid with hazard option"
                        exit 198
                }
        }
        else if "`hr'"!="" {
                local hazard "hazard"
	}

	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 weibull model:"
		}
		else {
			`nlog' di as txt _n "Fitting comparison weibull model:"
		}

		qui weibull_s `t' `rhs' `w' if `touse', nocoef `hazard' /*
		*/ dead(`dead') `topt'  `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 weibull model:"
		qui weibull_s `t' `rhs' `w' if `touse', nocoef `hazard' /*
		*/ dead(`dead') `topt' `constant' `mlopts' 
		tempname bp
		mat `bp' = e(b) 
	}
	
	if "`constant'"!="" {
		local skip = "skip"
		`nlog' di as txt _n "Fitting full model:"
	}
	
	global EREGd `dead'
	global EREGt `t'
	global EREGt0 `t0'
	if "`shared'" != "" {
		global EREG_by "`touse' `shared'"
	}
	
	if "`rhs'" != "" & "`skip'"=="" { 
		`nlog' di ""
		`nlog' di as txt "Fitting constant-only model:"
		ml model `mlmethod' `mlprog' (`t': `t'=, ) /*
			*/ /ln_p /ln_the /* 
			*/ `w' if `touse', /*
			*/ 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' ,  `constant') /*
		*/ /ln_p /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'"'
	}
	if "`hazard'"=="" {
		est local title "accelerated failure-time form"
		est local frm2 "time"
	}
	else {
		est local title "log-relative hazard form"
		est local frm2 "hazard"
	}
	est local fr_title `frtitle'
	est local predict `pred'
	est local cmd weibullhet2
	est local t0 "`t0'"
       	est local dead `sdead'
	est local stcurve="stcurve"
	est scalar aux_p = exp([ln_p]_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 weibullhet2
	macro drop EREGd EREGt EREGt0 
	
	Display, level(`level') `coef' `header' `tr' `hr'
	error `e(rc)'
end

program define GetFrailty, rclass
	local frailty `"`1'"'
	local hazard `"`2'"'
	if "`frailty'"=="" {
		di as err "must specify -frailty(gamma | invgauss)-"
		exit 198
	}
	
	local l = length("`frailty'")
	if substr("gamma",1,max(1,`l')) == "`frailty'" {
		if "`hazard'"=="" {
		        return local prog weibhet_glfa
		}
		else {
			return local prog weibhet_glf
		}
		return local frt "Gamma frailty"
		return local pred weibhet_gp 
	}
	else if substr("invgaussian",1,max(1,`l')) == "`frailty'" {
		if "`hazard'"=="" {
		        return local prog weibhet_ilfa
		}
		else {
			return local prog weibhet_ilf
		}
		return local frt "Inverse-Gaussian frailty"
		return local pred weibhet_ip
	}
	else {
		di as err "unknown distribution frailty(`frailty')"
		exit 198
	}
end
		
program define Display 
	if "`e(frm2)'"=="hazard" {
                local options "HR"
        }
        else    local options "TR"

	syntax [, Level(integer $S_level) noCOEF noHEADer `options']
	if "`coef'"=="" {
	        if "`hr'"!="" {
                        local eform "eform(Haz. Ratio)"
	        }
        	else if "`tr'"!="" {
                        local eform "eform(Tm. Ratio)"
                }

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

		_diparm ln_p, level(`level') 
		_diparm ln_the, level(`level') 
		di as txt "{hline 13}{c +}{hline 64}"
		_diparm ln_p, level(`level') noprob exp label("p")
		_diparm ln_p, level(`level') noprob function(exp(-@)) /*
			*/ derivative(-exp(-@)) label("1/p")
		_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