*! version 1.9.0 PR 27Feb1999.  STB-49 sg81.1
program define mfracpol, eclass
/*
1.9.0/PR:- _fracord replaced with frac_ord, _fracrep with frac_rep.
	 - Handling of adjustment and #unique values for each predictor
	   updated to use new frac_adj ado-file.
1.8.9/PR:- translated to Stata 6
*/
version 6
if "`1'" == "" | "`1'"=="," {
	if e(fp_cmd2)!="mfracpol" { error 301 }
	frac_rep "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
}
/*
	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) /*
 */ DFDefault(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 {
	marksample touse
	markout `touse' `varlist' `dead'
	frac_wgt "`exp'" `touse' "`weight'"
	local wgt `r(wgt)'				/* [`weight'`exp'] */
	count if `touse'
	local nobs = r(N)
}
/*
	Rearrange order of variables in varlist
*/
if "`xorder'"=="" { local xorder "+" }
frac_ord `varlist' `wgt' if `touse', order(`xorder') `regopt' cmd(`cmd') `fp01'
local varlist `s(names)'
gettoken lhs rhs : varlist
local nx : word count `rhs'
local i 1
while `i'<=`nx' {
/*
	Store original order of each RHS variable
*/
	local o`i' `s(pos`i')'
	local i = `i'+1
}
/*
	Initialisation.
*/
tokenize `rhs'
local i 1
while "``i''"!="" {
	local alp`i' .05	/* default FP selection level */
	local h`i' ``i''	/* names of H(xvars) 	*/
	local n`i' ``i''	/* names of xvars 	*/
	local po`i' 1 		/* to be final power	*/
	local sel`i' 1		/* default var selection level */
/*
	Remove old I* variables
*/
	frac_mun `n`i'' purge
	local i=`i'+1
}
/*
	Adjustment
*/
frac_adj "`adjust'" "`rhs'" `touse'
local i 1
while `i'<=`nx' {
	if "`r(adj`i')'"!="" {
		local adj`i' adjust(`r(adj`i')')
	}
	local uniq`i'=r(uniq`i')
	local i=`i'+1
}
/*
	Set up degrees of freedom for each variable
*/
if "`df'"!="" {
	frac_dis "`df'" df 1 . "`rhs'"
	local i 1
	while `i'<=`nx' {
		if "${S_`i'}"!="" { local df`i' ${S_`i'} }
		local i=`i'+1
	}
}
/*
	Assign default df for vars not so far accounted for.
	Give 1 df if 2-3 distinct values, 2 df for 4-5 values,
	dfdefault df for >=6 values.
*/
local i 1
while `i'<=`nx' {
	if "`df`i''"=="" {
		if `uniq`i''<=3 	{ local df`i' 1 }
		else if `uniq`i''<=5 	{ local df`i'=min(2,`dfdefau') }
		else 			{ local df`i' `dfdefau' }
	}
	local i=`i'+1
}
/*
	Set up FP selection level (alpha) for each variable
*/
if "`alpha'"!="" {
	frac_dis "`alpha'" alpha 0 1 "`rhs'"
	local i 1
	while `i'<=`nx' {
		if "${S_`i'}"!="" { local alp`i' ${S_`i'} }
		local i=`i'+1
	}
}
/*
	Set up selection level for each variable
*/
if "`select'"!="" {
	frac_dis "`select'" select 0 1 "`rhs'"
	local i 1
	while `i'<=`nx' {
		if "${S_`i'}"!="" { local sel`i' ${S_`i'} }
		local i=`i'+1
	}
}
/*
	Individual FP powers for variables.
*/
if "`xpowers'"!="" {
	frac_dis "`xpowers'" xpowers "`rhs'"
	local i 1
	while `i'<=`nx' {
		if "${S_`i'}"!="" { local xpow`i' ${S_`i'} }
		local i=`i'+1
	}
}
/*
	Vars with zero option
*/
if "`zero'"!="" {
	tokenize `zero'
	while "`1'"!="" {
		frac_in `1' "`rhs'"
		local zero`s(k)' "zero"
		mac shift
	}
}
/*
	Vars with catzero option
*/
if "`catzero'"!="" {
	tokenize `catzero'
	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 frac_154
*/
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 frac_rep to display
*/
	local npows: word count `p'
	local fd`i'=2*`npows'-(`npows'==1 & "`p'"=="1")+("`c'"!="")
	local id`i'=`df`o''+("`c'"!="")
	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 1
while `i'<=`nx' {
	local o `o`i''

	est scalar fp_fd`i'=`fd`i''	/* final degrees of freedom */
	est scalar fp_id`i'=`id`i''	/* initial degrees of freedom */
	est scalar fp_al`i'=`alp`o''	/* FP selection level */
	est scalar fp_se`i'=`sel`o''	/* var selection levl */

	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
	}
	local i=`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"
frac_rep "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