*! version 1.5.0 PR 2Feb98 STB-43 sg81
program define mfracpol
version 5.0
if "`1'" == "" | "`1'"=="," {
	if "$S_E_fp"!="mfracpol" { error 301 }
	_fracrep "fractional polynomial" "  df  " "Powers"
	exit
}
local cmd `1'
mac shift
capture _fracchk `cmd'
if _rc {
	di in red "Requires -fracpoly- version 1.5.1 or above (released with STB-41)"
	exit 198
}
if $S_1 {
	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_2
local glm $S_3
local qreg $S_4
local xtgee $S_5
local normal $S_6

local varlist "req ex min(2)"
local if "opt"
local in "opt"
local weight "aweight fweight iweight pweight"
#delimit ;
local options "ALpha(str) ADDpowers(str) CATzero(str) CYCles(int 5) DEAD(str)
 FIXpowers(str) noSCAling POwers(str) SELect(str) DF(str) XOrder(str)
 XPowers(str) ZERo(str) *" ;
#delimit cr
parse "`*'"
_fraccox "`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')" }
/*
	Rearrange order of variables in varlist
*/
if "`xorder'"=="" { local xorder "+" }
_fracord `varlist' `if' `in' [`weight' `exp'], order(`xorder') `regopt' /*
 */ cmd(`cmd')
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 n`nx' "`1'"	/* names of xvars 	*/
	local h`nx' "`1'"	/* names of H(xvars) 	*/
	local df`nx' 1		/* default df for xvar 	*/
	local po`nx' 1 		/* to be final power	*/
	local alp`nx' .05	/* default FP selection level */
	local sel`nx' 1		/* default var selection level */
/*
	Store original order of each RHS variable
*/
	local i=`nx'+1
	local o`nx' ${S_`i'}
/*
	Check for non-unique first 6 characters of var names
*/
	if length("`1'")>=6 {
		local l=`l'+1
		local l`l' "`1'"
	}
	mac shift
}
if `l'>1 {
	local longn
	local i 0
	while `i'<`l' {
		local i=`i'+1
		local li = substr("`l`i''",1,6)
		local j `i'
		while `j'<`l' {
			local j=`j'+1
			if "`li'"==substr("`l`j''",1,6) {
				local longn "`longn' `l`i'',`l`j''"
			}
		}
	}
}
if "`longn'"!="" {
di in red _n "ambiguous variable names" _n
di in bl "The first 6 characters of these pairs of variable names are the same:"
di in gr _n "`longn'" _n
di in bl "You must ensure that each variable is uniquely identified by the"
di in bl "first 6 characters of its name."
di in bl _n "Please rename any variables requiring it and rerun the analysis."
	exit 111
}
/*
	Set up degrees of freedom for each variable
*/
if "`df'"!="" {
	_fracdis "`df'" df 1 4 "`rhs'"
	local j 0
	while `j'<`nx' {
		local j=`j'+1
		if "${S_`j'}"!="" { local df`j' ${S_`j'} }
	}
}
/*
	Set up FP selection level (alpha) for each variable
*/
if "`alpha'"!="" {
	_fracdis "`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'"!="" {
	_fracdis "`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'"!="" {
			_fracin `1' "`rhs'"
			if `df${S_1}'<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'"!="" {
		_fracin `1' "`rhs'"
		local zero$S_1 "zero"
		mac shift
	}
}
/*
	Vars with catzero option
*/
if "`catzero'"!="" {
	unabbrev "`catzero'"
	parse "$S_1", parse(" ")
	while "`1'"!="" {
		_fracin `1' "`rhs'"
		local catz$S_1 "catzero"
		local zero$S_1 "zero"		/* catzero implies zero */
		mac shift
	}
}
/*
	Check for missing values in lhs, rhs and model vars.
*/
tempvar touse
quietly {
	mark `touse' [`weight' `exp'] `if' `in'
	markout `touse' `lhs' `rhs' `dead'
	_fracwgt "`exp'" `touse' "`weight'"
	local wgt $S_2				/* [`weight'`exp'] */
	count if `touse'
	local nobs = _result(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
			if "`rhs2'"=="" { local rhs2 "`h`j''" }
			else local rhs2 "`rhs2' `h`j''"
		}
		if `initial' {
			if "`rhs2'"!="" {
				local fixed "base(`rhs2')"
			}
			else local fixed
			qui fracsel `cmd' `lhs' `ni' `wgt' if `touse', /*
			 */ df(1) `fixed' select(1) `regopt'
			di in gr _n /*
			 */ "Deviance for model with all terms " /*
			 */ "linear = " in ye %9.3f $S_3 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 {
			if "`xpow`i''"!="" { local pw "powers(`xpow`i'')" }
			else local pw "powers(`powers')"
			fracsel `cmd' `lhs' `ni' `wgt' if `touse', `dfi' /*
			 */ `pw' `zero`i'' `catz`i'' `fixed' `h' /*
			 */ `scaling' `fpopt' `regopt' `pvalopt'
		}
		else {
			fracsel `cmd' `lhs' `ni' `wgt' if `touse', df(1) /*
			 */ `fixed' `h' `regopt' `pvalopt'
		}
		local h`i' "$S_1"
		local dev $S_3
		local p "$S_2"
		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 "$S_1"!="" {
			local rhs1 "`rhs1' $S_1"
		}
		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 _n "After cycle " in ye `it' /*
		 */ in gr ", deviance = " in ye %9.3f `dev'
	}
}
if `nx'>1 {
	local s
	if `it'!=1 { local s "s" }
	if `stable' { di _n in gr /*
	 */ "Fractional polynomial fitting algorithm converged" _cont }
	else di _n in gr "No convergence" _cont
	di in gr " after " in ye `it' in gr " cycle`s'."
}
/*
	Store results
*/
local finalvl	/* list of variable names constituting 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''
/*
	Remove variables left behind by fracpoly
*/
	local name=substr("`x'",1,6)
	cap drop `name'_0
	cap drop `name'_1
	cap drop `name'_2
/*
	Create FP vars as necessary
*/
	if trim("`p'")!="." {
		if trim("`p'")=="1" & "`c'"=="" {
			global S_1 `x'
		}
		else qui fracgen `x' `p', `scaling' `z' `c'
		local finalvl "`finalvl' $S_1"
	}
/*
	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 _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 */
	global S_E_k`i' `po`o'' /* "powers" for ith predictor */
}
global S_E_dist `dist'
global S_E_wgt `weight'
global S_E_exp "`exp'"
global S_E_depv `lhs'
global S_E_dev `dev'
global S_E_rhs	/* deliberately blank for consistency with fracpoly */
global S_E_opts `regopt'
global S_E_fvl `finalvl'
global S_E_nx `nx'
global S_E_t1t "Fractional Polynomial"
_fracrep "fractional polynomial" "  df  " "Powers"
global S_E_fp "mfracpol"
end

* Constructed according to pseudocode analysis in fpsel2.alg
*?? Program fracsel needs tidying up.
*?? Weight handling, is full parse necessary, is it stable and safe,
*?? compatibility with similar routine in mrsnb, general standards.
program define fracsel
local cmd `1'
mac shift
local varlist "req ex min(2)"
local if "opt"
local in "opt"
local weight "fweight aweight"
#delimit ;
local options "ALpha(real .05) SELect(real .05) noHEad DF(int -1)
 POwers(string) BAse(string) noTIDy CATzero *" ;
#delimit cr
parse "`*'"
local cz="`catzero'"!=""
local omit=(`select'<1)
if `df'>4 | `df'==0 | `df'<-1 {
	di in red "invalid df"
	exit 198
}
local meq2 = `df'>2 | `df'==-1
if `df'==4 { local df -1 }
if `df'==1 {					/* linear */
	if "`powers'"!="" {
		di in red "powers() invalid with df(1)"
		exit 198
	}
	local fixpowe 1
}
else {
	if "`powers'"=="" { local powers "-2,-1,-.5,0,.5,1,2,3" }
	local powers "powers(`powers')"
	local degree=1+`meq2'
	local degree "degree(`degree')"
}

if "`weight'"!="" { 
	local weight "[`weight'`exp']"
}
parse "`varlist'", parse(" ")
local lhs `1'
mac shift
local rhs
if `meq2' {
	local two "Powers   DDif(m=2)    P"
}
else {
	local two "                       "
	}
if "`head'"=="" {
	di in gr _col(11) "--m = 1 vs linear--" /*
	 */ _col(34) "--- m = 2 vs m = 1 ---  Final     P(FP models)"
	di in gr "Variable" _col(11) "Power  Dev.diff.  P" /*
	 */ _col(34) "Powers   Dev.diff.  P   Powers    m = 1  m = 2"
	di in gr _dup(79) "-"
}
local i 1
while "``i''"!="" { 
	local n ``i''
	local pwrs 1
	local name = substr("`n'", 1, 6)
	qui fracpoly `cmd' `lhs' `n' `fixpowe' `base' `if' `in' /*
	 */ `weight', `degree' `powers' `options' `catzero'
	local normal=($S_E_dist==0)
	local nobs $S_E_nobs
	local dev0 $S_E_d0
	local devlin $S_E_dlin
	local dev `devlin'
	local devm1 $S_E_d1
	local pm1 $S_E_p1
	local rdf $S_E_rdf

	local n1=1+`cz'
	local n2=`rdf'-`n1'
	local d=`dev0'-`devlin'
	_fracpv `normal' "`weight'" `nobs' `d' `n1' `n2'
	local P0lin $S_1

	if `df'==1 {
		local pm1 1
		local P0m1 `P0lin'
	}
	else {
		local n1 1
		local n2=`rdf'-2-`cz'
		local d=`devlin'-`devm1'
		_fracpv `normal' "`weight'" `nobs' `d' `n1' `n2'
		local Plinm1 $S_1
		local n1=2+`cz'
		local n2=`rdf'-`n1'
		local d=`dev0'-`devm1'
		_fracpv `normal' "`weight'" `nobs' `d' `n1' `n2'
		local P0m1 $S_1
		if `meq2' {
			local n1 2
			local n2=`rdf'-4-`cz'
			local pm21: word 1 of $S_E_p2
			local pm22: word 2 of $S_E_p2
			local devm2 $S_E_d2
			local d=`devm1'-`devm2'
			_fracpv `normal' "`weight'" `nobs' `d' `n1' `n2'
			local Pm1m2 $S_1
			local n1=4+`cz'
			local n2=`rdf'-`n1'
			local d=`dev0'-`devm2'
			_fracpv `normal' "`weight'" `nobs' `d' `n1' `n2'
			local P0m2 $S_1
		}
	}
	di in gr "``i''" _col(12) in ye `pm1' _col(16) _cont
	if `pm1'==1 {
		local Plinm1 1
		if `df'!=1 {
			di in ye %8.3f 0 %7.3f 1 _cont
		}
		else di _skip(14) " " _cont
	}
	else {
		di in ye %8.3f `dev'-`devm1' /*
		 */ %7.3f `Plinm1' _cont
	}
	di "   " _cont

	local done 1
	if `meq2' {
		local pm = substr("`pm21',`pm22'     ",1,8)
		di in ye "`pm'" %8.3f `devm1'-`devm2' /*
		 */ %7.3f `Pm1m2' _cont

		if `omit' & (`P0m2'>=`select' & `P0m2'!=.) {
/*
	`Dropping' RHS variable since 4+cz df test of m=2 is non-sig.
        Note: using variable selection level here, not FP sel level.
*/
			local dev `dev0' 
			local pwrs "."
			local n
		}
		else {
			local done 0
			local omit 0
		}
		if !`done' {
/*
        Test m=2 vs m=1 at FP sel level
*/
			if `Pm1m2'<`alpha' & `Pm1m2'!=. {
				local done 1
				local dev `devm2'
				qui fracgen `n' `pm21' `pm22', /*
				 */ `options' `catzero' replace
				local n $S_1
				local pwrs "`pm21' `pm22'"
			}
		}
	}
	if `df'==2 | !`done' {
		local done 0
		if `omit' {
/*
        Test for `dropping' RHS variable at variable selection level
        using 2 df test of m=1.
*/
                        if `P0m1'>`select' & `P0m1'!=. {
				local done 1
				local dev `dev0' 
				local pwrs "."
				local n
			}
		}
		if !`done' {
/*
        Test m=1 vs linear using FP sel level
*/
                        if `Plinm1'<=`alpha' {
				local dev `devm1'
				qui fracgen `n' `pm1', /*
				 */ `options' `catzero' replace
				local pwrs `pm1'
				local n $S_1
			}
			else {
				local dev `devlin'
			}
		}
	}
	if `df'==1 {
		local done 0
		if `omit' {
/*
        Test for `dropping' RHS variable at variable selection level
        using 1 df test of beta=0.
*/
                        if `P0lin'>`select' & `P0lin'!=. {
				local done 1
				local dev `dev0' 
				local pwrs "."
				local n
			}
		}
		if !`done' {
			local dev `devlin'
		}
	}
	if !`meq2' { di _skip(23) _cont }
	di " " substr("`pwrs'       ",1,8) %7.3f `P0m1' _cont
	if `meq2' { di in ye %7.3f `P0m2' }
	else di
	local rhs "`rhs' `n'"
	local i = `i'+1
}
global S_1 `rhs'
global S_2 "`pwrs'"
global S_3 `dev'
global S_4 `n'
end