*! version 1.0.0  10oct2007
program scoretest_cox, rclass
	version 10

	// check for unsupported estimation results
	if `"`e(cmd2)'"' != "stcox" {
		error 301
	}
	if `"`e(cmd)'"' == "stcox_fr" {
		di as err "score test not allowed with shared frailty"
		exit 322
	}
	if `"`e(vce)'"' != "oim" {
		di as err "score test not allowed with vce(`e(vce)')"
		exit 322
	}
	if `"`e(texp)'"' != "" {
		di as err "score test not allowed with the option tvc()"
		exit 322
	}
	if `"`e(prefix)'"' != "" {
		di as err "score test not allowed with `e(prefix)' results"
		exit 322
	}

	// these variable identify which coefficients are going to be tested
	syntax varlist [, showcox ]

	tempvar touse
	mark `touse' if e(sample)

	// building the options to be passed directly to -cox-
	if `"`e(wtype)'"' != "" {
		local wt [`e(wtype)'`e(wexp)']
	}
	if `"`e(strata)'"' != "" {
		local strata strata(`e(strata)')
	}
	if `"`e(offset)'"' != "" {
		local offset offset(`e(offset)')
	}
	if `"`e(method)'"' == "partial" {
		local ties exactp
	}
	else if `"`e(method)'"' == "marginal" {
		local ties exactm
	}
	else	local ties `e(method)'

	tempname matb results b1 b0 b g h
	mat `b1' = e(b)
	local cmdline `"`e(cmdline)'"'

	// hold the current results and auto-restore when we exit
	_est hold `results', restore nullok

	* full contains variables in the full model
	* tst contains the variables to test
	* rest contains variables in the restricted model
	local full : colnames `b1'
	local tst : list uniq varlist
	local rest : list full - tst

	* check that tested variables are in the model
	local invalid : list tst - full
	if `"`invalid'"' != "" {
		gettoken invalid : invalid
		di as err "`invalid' not found in the model"
		exit 111
	}

	* fit cox model on the rest, and construct the matrix
	* of coefficients b for the restricted model
	* parse options first
	ParseOpt `cmdline'
	local options `"`r(options)'"'
	qui stcox `rest' if `touse', `options'
	mat `b0' = e(b)
	local npar = colsof(`b1')
	local rnames : colfullnames `b0'

	// fill in the parameter matrix for the null Hypothesis
	mat `b' = J(1, `npar', 0)
	foreach nm in `rnames'{
		local j0 = colnumb(`b0', "`nm'")
		local j1 = colnumb(`b1', "`nm'")
		matrix `b'[1,`j1'] = `b0'[1,`j0']
	}

	* get the gradient and the Hessian for the Cox model at H0
	local coxcmd cox _t `full' `wt' if `touse',	///
				t0(_t0)			///
				dead(_d)		///
				matfrom(`b')		///
				gradient(`g')		///
				hessian(`h')		///
				`ties'			///
				`strata' 
	if "`showcox'" != "" {
		di as txt "{p 0 0 2}calling cox -> "	///
		   in white `"`coxcmd'"'		///
		   "{p_end}"
	}
	quietly `coxcmd'

	// compute the score test statistic
	tempname val pval
	mata: sctest_util()
	local  df = `npar' - colsof(`b0')
	scalar `val'  = r(val)
	scalar `pval' = chi2tail(`df', r(val))
	
	* display results
	local i 0
	di
	foreach var of local tst {
		di as txt " (" %2.0f `++i' ")" as res _col(8) "`var' = 0"
	}
	di
	di as txt _col(12) "chi2(" %3.0f `df' ") ="	as res %8.2f `val'
	di as txt _col(10) "Prob > chi2 =  "		as res %8.4f `pval'

	* return results
	return scalar df = `df'
	return scalar chi2 = `val'
	return scalar p = `pval'
end

program ParseOpt, rclass
	syntax anything [,*]
	return local options `"`options'"'
end

mata:
void sctest_util()
{
	real matrix u1
	real matrix H
	real matrix V

	u1 = st_matrix(st_local("g"))
	H = st_matrix(st_local("h"))
	V = cholinv(H)
	m = u1*V*u1'
	a  = m[1,1]
	st_numscalar("r(val)", a)
}
end

exit