*! version 1.8.9 PR 21Dec98.
program define mfracpol, eclass
version 6
if "`1'" == "" | "`1'"=="," {
	if e(fp_cmd2)!="mfracpol" { error 301 }
	_fracrep "fractional polynomial" "  df  " "Powers"
	exit
}
gettoken cmd 0 : 0
frac_chk `cmd'
if `s(bad)' {
	di in red "invalid or unrecognised command, `cmd'"
	exit 198
}
*!! st commands disabled
/*
if $S_7 {
	_fracst m `cmd' `*'
	exit
}
*/
/*
	dist=0 (normal), 1 (binomial), 2 (poisson), 3 (cox), 4 (glm),
	5 (xtgee), 6(ereg/weibull).
*/
local dist `s(dist)'
local glm `s(isglm)'
local qreg `s(isqreg)'
local xtgee `s(isxtgee)'
local normal `s(isnorm)'

syntax varlist(min=2) [if] [in] [aw fw pw iw] [, /*
 */ ADDpowers(str) ADJust(str) ALpha(str) ALL CATzero(str) DF(str) /*
 */ DFMax(int 4) CYCles(int 5) DEAD(str) FIXpowers(str) FP01 noSCAling /*
 */ POwers(str) SELect(str) SEQuential XOrder(str) XPowers(str) ZERo(str) * ]
frac_cox "`dead'" `dist'
/*
	Process options
*/
local regopt `options'
if "`dead'"!="" { local regopt "`regopt' dead(`dead')" }
if "`powers'"==""  { local powers "-2,-1,-.5,0,.5,1,2,3" }
if "`addpowe'"!="" { local fpopt "`fpopt' addp(`addpowe')" }
if "`fixpowe'"!="" { local fpopt "`fpopt' fixp(`fixpowe')" }
/*
	Check for missing values in lhs, rhs and model vars.
*/
tempvar touse
quietly {
	mark `touse' [`weight' `exp'] `if' `in'
	markout `touse' `varlist' `dead'
	frac_wgt "`exp'" `touse' "`weight'"
	local wgt `r(wgt)'				/* [`weight'`exp'] */
	count if `touse'
	local nobs = _result(1)
}
/*
	Rearrange order of variables in varlist
*/
if "`xorder'"=="" { local xorder "+" }
_fracord `varlist' `wgt' if `touse', order(`xorder') `regopt' /*
 */ cmd(`cmd') `fp01'
local varlist "$S_1"
parse "`varlist'", parse(" ")
local lhs `1'
mac shift
local rhs "`*'"
/*
	Initialisation.
	Count RHS vars, store var names in n1, n2, ...,
	and store default df and powers.
*/
local nx 0
local l 0
while "`1'"!="" {
	local nx = `nx'+1
	local alp`nx' .05	/* default FP selection level */
	local h`nx' "`1'"	/* names of H(xvars) 	*/
	local n`nx' "`1'"	/* names of xvars 	*/
	local po`nx' 1 		/* to be final power	*/
	local sel`nx' 1		/* default var selection level */
/*
	Remove old I* variables for xvar
*/
	frac_mun `n`nx'' purge
/*
	Store original order of each RHS variable
*/
	local i=`nx'+1
	local o`nx' ${S_`i'}

	mac shift
}
/*
	Adjustment
*/
if "`adjust'"!="" {
	frac_dis "`adjust'" adjust "`rhs'"
	local j 0
	while `j'<`nx' {
		local j=`j'+1
		local adj`j' ${S_`j'}
	}
}
/*
	Set up degrees of freedom for each variable
*/
if "`df'"!="" {
	frac_dis "`df'" df 1 . "`rhs'"
	local j 0
	while `j'<`nx' {
		local j=`j'+1
		if "${S_`j'}"!="" { local df`j' ${S_`j'} }
	}
}
/*
	Assign default df for vars not so far accounted for.
	Award 1 df if 2-3 distinct values, 2 df for 4-5 values,
	dfmax df for >=6 values.

	Adjustment: adjust by default if >2 distinct values, not otherwise.
*/
local j 1
while `j'<=`nx' {
	qui inspect `n`j'' if `touse'
	local uniq=_result(7)
	if "`df`j''"=="" {
		local n `n`j''
		if `uniq'==1 {
			noi di in red "`n`j'' does not vary."
			exit 198
		}
		else if `uniq'<=3 { local df`j' 1 }
		else if `uniq'<=5 { local df`j' 2 }
		else local df`j' `dfmax'
	}
	if "`adj`j''"!="" {
		if "`adj`j''"=="no" {
			local adj`j'
		}
		else local adj`j' adjust(`adj')
	}
	else if `uniq'>2 {
		local adj`j' adjust(mean)
	}
	local j=`j'+1
}
/*
	Set up FP selection level (alpha) for each variable
*/
if "`alpha'"!="" {
	frac_dis "`alpha'" alpha 0 1 "`rhs'"
	local j 0
	while `j'<`nx' {
		local j=`j'+1
		if "${S_`j'}"!="" { local alp`j' ${S_`j'} }
	}
}
/*
	Set up selection level for each variable
*/
if "`select'"!="" {
	frac_dis "`select'" select 0 1 "`rhs'"
	local j 0
	while `j'<`nx' {
		local j=`j'+1
		if "${S_`j'}"!="" { local sel`j' ${S_`j'} }
	}
}
/*
	Powers for individual FP terms.
*/
if "`xpowers'"!="" {
	parse "`xpowers'", parse(",")
	local ncl 0 /* # of comma-delimited clusters */
	while "`1'"!="" {
		if "`1'"=="," { mac shift }
		local ncl=`ncl'+1
		local clust`ncl' "`1'"
		mac shift
	}
	if "`clust`ncl''"=="" { local ncl=`ncl'-1 }
	if `ncl'>`nx' {
		di in red "too many sets of powers() values specified"
		exit 198
	}
/*
	Disentangle each varlist:string cluster
*/
	local i 1
	while `i'<=`ncl' {
		parse "`clust`i''", parse("=:")
		local powk `3'
		unabbrev "`1'"
		parse "$S_1", parse(" ")
		while "`1'"!="" {
			frac_in `1' "`rhs'"
			if `df`s(k)''<2 {
				di in red "xpowers invalid for `1', it has 1 df"
				exit 198
			}
			local xpow$S_1 `powk'
			mac shift
		}
		local i = `i'+1
	}
}
/*
	Vars with zero option
*/
if "`zero'"!="" {
	unabbrev "`zero'"
	parse "$S_1", parse(" ")
	while "`1'"!="" {
		frac_in `1' "`rhs'"
		local zero`s(k)' "zero"
		mac shift
	}
}
/*
	Vars with catzero option
*/
if "`catzero'"!="" {
	unabbrev "`catzero'"
	parse "$S_1", parse(" ")
	while "`1'"!="" {
		frac_in `1' "`rhs'"
		local catz`s(k)' "catzero"
		local zero`s(k)' "zero"		/* catzero implies zero */
		mac shift
	}
}
/*
	Reserve names for H(predictors) by creating a dummy variable
	for each predictor which potentially needs transformation.
*/
local i 1
while `i'<=`nx' {
	if `df`i''>1 | "`zero`i''`catz`i''"!="" {
		frac_mun `n`i''
		local stub`i' `s(name)'
		qui gen byte `stub`i''_1=.
	}
	local i=`i'+1
}
/*
	Build FP model.
*/
local it 0
local initial 1
local stable 0 /* convergence flag */
while !`stable' & `it'<=`cycles' {
	local it = `it'+1
	local pwrs
	local rhs1
	local stable 1 /* later changed to 0 if any power or status changes */
	local lastch 0 /* becomes index of last var which changed status */
	local i 1
	while `i'<=`nx' {
		local ni "`n`i''"
		local dfi "df(`df`i'')"
/*
	Build up RHS2 from the i+1th var to the end
*/
		local rhs2
		local j `i'
		while `j'<`nx' {
			local j = `j'+1
			local rhs2 `rhs2' `h`j''
		}
		if `initial' {
			if "`rhs2'"!="" {
				local fixed "base(`rhs2')"
			}
			else local fixed
			qui frac_sel `cmd' `lhs' `ni' `wgt' if `touse', /*
			 */ df(1) `fixed' select(1) `regopt'
			di in gr _n /*
			 */ "Deviance for model with all terms " /*
			 */ "untransformed = " in ye %9.3f r(dev) in gr ", " /*
			 */ in ye `nobs' in gr " observations"
		}
		if "`rhs1'`rhs2'"!="" {
			local fixed "base(`rhs1' `rhs2')"
		}
		else local fixed
/*
	Vars with df(1) are straight-line
*/
		local pvalopt "alpha(`alp`i'') select(`sel`i'')"
		if `i'==1 { di }
		if `df`i''==1 {
			local pw
			local fpo
		}
		else {
			if "`xpow`i''"!="" { local pw "powers(`xpow`i'')" }
			else local pw "powers(`powers')"
			local fpo `fpopt'
		}
		if `df`i''==1 & `sel`i''==1 {	/* var is included anyway */
			local rhs1 `rhs1' `h`i''
			di in gr "[`ni' included with 1 df in model]" _n
		}
		else {
			if "`stub`i''"!="" { local n name(`stub`i'') }
			else local n
			frac_sel `cmd' `lhs' `ni' `wgt' if `touse', `dfi' /*
			 */ `pw' `zero`i'' `catz`i'' `fixed' `h' /*
			 */ `fpo' `regopt' `pvalopt' `n' `sequent'
			local h`i' `r(rhs)'
			local dev=r(dev)
			local p `r(pwrs)'
			if "`p'"!="`po`i''" {
				if `nx'>1 { local stable 0 }
				local po`i' "`p'"
				local lastch `i'
			}
			if "`pwrs'"=="" { local pwrs "`p'" }
			else local pwrs "`pwrs',`p'"
			if "`h`i''"!="" {
				local rhs1 `rhs1' `h`i''
			}
		}
		if `initial' {
			local h "nohead"
			local initial 0
		}
		local i = `i'+1
	}
	if `lastch'==1 { local stable 1 } /* 1 change only, at i=1 */
	if !`stable' {
		di in gr _dup(70) "-" _n "Cycle " in ye `it' in gr /*
		 */ ": deviance =   " in ye %9.3f `dev' _n in gr _dup(70) "-"
	}
}
if `nx'>1 {
	local s
	if `it'!=1 { local s "s" }
	if !`stable' {
		di _n in gr "No convergence" _cont
	}
	else di _n in gr /*
	 */ "Fractional polynomial fitting algorithm converged" _cont
	di in gr " after " in ye `it' in gr " cycle`s'."
}
if `stable' {
	di _n in gr "Transformations of covariates:" _n
}
/*
	Remove variables left behind by _frac154
*/
local i 1
while `i'<=`nx' {
	if "`stub`i''"!="" { cap drop `stub`i''* }
	local i=`i'+1
}
/*
	Store results
*/
local finalvl	/* predictors in final model */
local i 1
while `i'<=`nx' {
	local o `o`i''
	local p `po`o''
	local x `n`o''
	local z `zero`o''
	local c `catz`o''
	local a `adj`o''
/*
	Create FP vars as necessary, with new unique names
*/
	if trim("`p'")!="." {
		if trim("`p'")=="1" & "`z'`c'`a'"=="" {
			local namex `x'
		}
		else {
			frac_mun `x'
			local vn `s(name)'
			fracgen `x' `p' if e(sample), `all' `scaling' /*
			 */ `z' `c' `a' name(`vn')
			local namex `r(names)'
		}
		local finalvl `finalvl' `namex'
	}
/*
	Save stuff for _fracrep to use for display purposes
*/
	local npows: word count `p'
	local fdf=2*`npows'-(`npows'==1&"`p'"=="1")+("`c'"!="")
	local idf=`df`o''+("`c'"!="")
	global MFPfd`i' `fdf'		/* final degrees of freedom */
	global MFPdf`i' `idf'		/* initial degrees of freedom */
	global MFPal`i' `alp`o''	/* FP selection level */
	global MFPse`i' `sel`o''	/* var selection levl */
	local i=`i'+1
}
/*
	Estimate final (conditionally linear) model.
*/
qui `cmd' `lhs' `finalvl' `wgt' if e(sample), `regopt'
/*
	Store results
*/
global S_1 `finalvl'
global S_2 `dev'
local i 0
while `i'<`nx' {
	local i=`i'+1
	local o `o`i''
	global S_E_x`i' `n`o''	/* name of ith predictor in user order */
	est local fp_x`i' `n`o''
	global S_E_k`i' `po`o'' 	/* "powers" for ith predictor */
	est local fp_k`i' `po`o''
	if "`catz`o''"!="" {
		global S_E_c`i' 1
		est local fp_c`i' 1
	}
}
global S_E_dist `dist'
est scalar fp_dist=`dist'
global S_E_wgt `weight'
est local fp_wgt `weight'
global S_E_exp "`exp'"
est local fp_exp `exp'
global S_E_depv `lhs'
est local fp_depv `lhs'
global S_E_dev `dev'
est scalar fp_dev=`dev'
global S_E_rhs	/* deliberately blank for consistency with fracpoly */
est local fp_rhs
global S_E_opts `regopt'
est local fp_opts `regopt'
global S_E_fvl `finalvl'
est local fp_fvl `finalvl'
global S_E_nx `nx'
est scalar fp_nx=`nx'
global S_E_t1t "Fractional Polynomial"
est local fp_t1t "Fractional Polynomial"
_fracrep "fractional polynomial" "  df  " "Powers"
global S_E_fp "fracpoly"
est local fp_cmd "fracpoly"
global S_E_fp2 "mfracpol"
est local fp_cmd2 "mfracpol"

end