*! version 2.0.9 PR/EW 03Apr97.  STB-38 sbe15
program define xrigls
version 4.0
local varlist "req ex min(2) max(2)"
local if "opt"
local in "opt"
#delimit ;
local options "ALpha(real 0.05) COVArs(str) CEntile(str) CYcles(int 2)
 CV DETail FP(str) noGRaph noLEAve noTIdy PARam(str) POwers(str)
 ROpts(str) SAVing(str) SE noSELect ";
#delimit cr
local small 1e-6
parse "`*'"
if "`select'"=="noselect" { local omit "noomit" }
local dist bc
local dupper "N "
if "`centile'"=="" { local centile "3 97" }
parse "`centile'", parse(" ")
local i 0
while "`1'"!="" {
	local nc=`nc'+1
	local cent`nc' `1'
	local C "`C' __C`1'"
	mac shift
}
local ginit 1 /* i.e. normal distribution */
if "`powers'"!="" { local powers "powers(`powers')" }
if "`param'"!="" {
	parse "`param'", parse(" ")
	while "`1'"!="" {
		local fl=lower(substr("`1'",1,2))
		if "`fl'"=="cv" 	{ local cv "cv" } /* !! not working properly yet */
		else if "`fl'"=="ln" 	{ local lns "lns" }
		else {
			di in red "invalid param() option `1'"
			exit 198
		}
		mac shift
	}
}
if "`cv'"!="" { local cv "cv" }
if "`cv'"!="" { local slab "CV"}
else { local slab "SD" }
/*
	Extract regression options
*/
_jprxrpa "`ropts'" ms "regression options"
local ncl $S_1
local i 1
while `i'<=`ncl' {
	local param ${S_p`i'}
	local `param'opt "${S_c`i'}"
	local i = `i'+1
}
/*
	Extract base variables for each parameter
*/
_jprxrpa "`covars'" ms "covariates"
local ncl $S_1
local i 1
while `i'<=`ncl' {
	local param ${S_p`i'}
	local `param'base "${S_c`i'}" /* covariates -- check validity? */
	local i = `i'+1
}
/*
	Most of the next 5 stmts are probably unnecessary
*/
global S_trunc = "`trunc'"!=""
global S_cv = "`cv'"!=""
global S_tau = "`tau'"!=""
global S_lns = "`lns'"!=""
global S_gamma `ginit' /* !! fix this up -- needed for centcalc? */

parse "`varlist'", parse(" ")
tempvar y touse
local lhs "`1'"
local rhs "`2'"
tempvar y touse
mark `touse' `if' `in'
markout `touse' `lhs' `rhs' `mbase' `sbase'
/*
	Extract details of FP model
*/
_jprxrpa "`fp'" ms "FP model terms"
local ncl $S_1 /* could be 0 */
local i 1
while `i'<=`ncl' {
	local curve ${S_p`i'}
	local spec "${S_c`i'}"
	parse "`spec'", parse(" ")
	local first = lower("`1'")
	if "`first'"=="fix" {
		if "`curve'"=="m" | "`curve'"=="s" {
			di in red "invalid `spec'"
			exit 198
		}
		local `curve'fix "`curve'fix"
	}
	else if "`first'"=="df" {
		mac shift
		confirm num `*'
		if `1'<0 {
			di in red "invalid df"
			exit 198
		}
		local `curve'df `1'
	}
	else {
		if "`first'"=="powers" { mac shift }
		qui fracgen `rhs' `*' if `touse',name(X`curve') replace
		local `curve'vars $S_1
		local `curve'fixpow `*'
	}
	local i = `i'+1
}
/*
	Deal with df
*/
local c1 m
local c2 s
local i 1
while `i'<=2 {
	local curve `c`i''
	if "``curve'fixpow'"!="" {
		if "``curve'df'"!="" {
			di in red "invalid df" /* 'cos have fixpowers */
			exit 198
		}
		local `curve'pow ``curve'fixpow'
		local `curve'df -1	/* means have fixpowers */
	}
	else if "``curve'df'"=="" {
		if "`curve'"=="m" {
			local mdf 4	/* default for M */
		}
		else	local sdf 2	/* default for S */
	}
	local `curve'powS = 2+2*``curve'df' /* pos'n of powers in S_ */
	local i = `i'+1
}
quietly {
/*
	Calc y-fit, then sd from regression on abs residuals,
	then refit both using weights based on fitted sd.
*/
	local exp
/*
	!! [Boxcox option deleted]
	If Boxcox transformation is requested,
	find the geometric mean of y and transform y to standardized form.
*/
	gen `y' = `lhs' if `touse'

* ----------------------------------------------------------------------
/*
	Initialisation (cycle 0, if iteration required)
*/
	tempvar M S wt
	tempvar mse sse
	gen `wt' = 1 if `touse'
	local exp "[aweight=`wt']"
	if "`cv'"!="" { local expy "[aweight=`wt'*`M'^-2]" }
	else		local expy `exp'
	if "`mfixpow'"!="" | `mdf'==0 {
		regress `y' `mvars' `mbase', `mopt'
	}
	else {
		if "`mbase'"!="" { local base "base(`mbase')" }
		else local base
		frachoos regress `y' `rhs', `base' `powers' `omit' /*
		 */ `scaling' alpha(`alpha') `mopt' name(Xm) df(`mdf')  
		local mpow $S_4
		local mvars $S_5
		regress `y' `mvars' `mbase', `mopt'
	}
	if _result(5)<1 { error 2001 } /* need at least 1 residual df */
	local obs = _result(1)
	local ssr = _result(4)
	local rmse = _result(9)
	_devian `wt' `obs' `ssr' 1
	local dev $S_1
	local devold `dev'
	pred `M' if `touse'
	if "`se'"!="" { predict `mse' if `touse', stdp }
	if `sdf'==0 & "`cv'`sbase'"=="" {
		local spow " "
		gen `S' = `rmse'
	}
	else {
		local spowold " "
		tempvar absres
		if "`lns'"=="" {
			local f=sqrt(_pi/2)
			if "`cv'"=="" {
				gen `absres' = `f'*abs(`y'-`M')
			}
			else gen `absres' = `f'*abs(`y'/`M'-1)
		}
		else {
			local f 0.63518142
			if "`cv'"=="" {
				gen `absres' = `f'+ln(abs(`y'-`M'))
			}
			else gen `absres' = `f'+ln(abs(`y'/`M'-1))
		}
		if "`sfixpow'"!="" | `sdf'==0 {
			regress `absres' `svars' `sbase', `sopt'
		}
		else {
			if "`sbase'"!="" { local base "base(`sbase')" }
			else local base
			frachoos regress `absres' `rhs', `base' `omit' /*
			 */ `powers' `scaling' alpha(`alpha') /*
			 */ `sopt' name(Xs) df(`sdf') 
			local spow $S_4
			local svars $S_5
			regress `absres' `svars' `sbase', `sopt'
		}
		pred `S' if `touse'
		if "`se'"!="" { predict `sse' if `touse', stdp }
		if "`lns'"!="" { replace `S'=exp(`S') }
		replace `wt' = `S'^-2 if `touse'
	}
/*
	End of `sdf'==0 conditional
*/
	noi di _n in gr _col(11) "--- FP Powers ---" _n /*
	 */ "Cycle" _col(11) "Mean" _col(22) "`slab'" /*
	 */ _col(33) " Deviance" _col(44) "   Change" /*
	 */ _col(55) "Residual SS" _n _dup(65) "-"

        noi di in ye 0 _col(11) in ye "`mpow'" /*
	 */ _col(22) "`spowold'" _col(33) in ye %9.3f `dev' /*
	 */ _col(44) in ye %9.3f 0 _col(55) in ye %9.0g `ssr'
	if !(`sdf'==0 & "`cv'`sbase'"=="") {
/*
	Cycle between mean fit and SD fit
*/
		local it 1
		while `it'<=`cycles' {
			local lastit = (`it'==`cycles')
/*
	Weighted fit for mean
*/
			if "`mfixpow'"!="" | `mdf'==0 {
				regress `y' `mvars' `mbase' `expy', `mopt'
			}
			else {
				if "`mbase'"!="" { local base "base(`mbase')" }
				else local base
				frachoos regress `y' `rhs' `expy', `omit' /*
				 */ `base' `powers' `scaling' /*
				 */ alpha(`alpha') `mopt' name(Xm) df(`mdf')
				local mpow $S_4
				local mvars $S_5
				regress `y' `mvars' `mbase' `expy', `mopt'
			}
			local ssr = _result(4)
			local rmse = _result(9)
			drop `M'
			cap drop `mse'
			pred `M' if `touse'
			if "`se'"!="" { predict `mse' if `touse', stdp }
			if "`cv'"!="" { local mf `M' }
			else		local mf 1
			_devian `wt' `obs' `ssr' `mf'
			local dev $S_1
			local spowold `spow'
			if !`lastit' {
/* ************************************************************************

	Weighted fit for SD

	Note use of `exp' (i.e. weight = `wt') in following regression,
	to allow for the variance of the absolute residuals
	(see Carroll and Ruppert p 80, expression (3.20).
*/
	if "`lns'"=="" {
		if "`cv'"=="" {
			replace `absres' = `f'*abs(`y'-`M')
		}
		else replace `absres' = `f'*abs(`y'/`M'-1)
	}
	else {
		if "`cv'"=="" {
			replace `absres' = `f'+ln(abs(`y'-`M'))
		}
		else replace `absres' = `f'+ln(abs(`y'/`M'-1))
	}
	if "`sfixpow'"!="" | `sdf'==0 {
		regress `absres' `svars' `sbase' `e', `sopt'
	}
	else {
		if "`sbase'"!="" { local base "base(`sbase')" }
		else local base
		frachoos regress `absres' `rhs' `e', `omit' /*
		 */ `base' `powers' `scaling' /*
		 */ alpha(`alpha') `mopt' name(Xs) df(`sdf')
		local spow $S_4
		local svars $S_5
		regress `absres' `svars' `sbase' `e', `sopt'
	}
	drop `S' 
	cap drop `sse'
	pred `S' if `touse'
	if "`se'"!="" { predict `sse' if `touse', stdp }
	if "`lns'"!="" { replace `S'=exp(`S') }
	sum `S'
	if _result(5)<=0 {
		noi di in red "negative fitted standard deviation"
		exit 2002
	}
	replace `wt' = `S'^-2 if `touse'
************************************************************************
			}
			noi di in ye `it' _col(11) "`mpow'" /*
			 */ _col(22) "`spowold'" _col(33) %9.3f `dev' /*
			 */ _col(44) in ye %9.3f `dev'-`devold' /*
			 */ _col(55) in ye %9.0g `ssr'
			local devold `dev'
			local it=`it'+1
		} /* end of cycling loop */
*-----------------------------------------------------------------------------
	}
/*
	Standardized residuals (Z-scores)
*/
	tempvar Z
	if "`cv'"=="" {
		gen `Z' = (`y'-`M')/`S' if `touse'
	}
	else	gen `Z' = (`y'/`M'-1)/`S' if `touse'
/*
	Reference interval
*/
	centcalc `M' `S', centile(`centile') prefix(__C) /*
	 */ dist(`dist') gamma(1) `cv'
}
/*
	SEs of centile values.
*/
	quietly {
	count if `touse'
	local nobs=_result(1)
	if "`se'"!="" {
		loc i 0
		while `i'<`nc' {
		    loc i=`i'+1
		    tempvar c`i'se 
		    local z=invnorm(`cent`i''/100)
		    gen `c`i'se'=sqrt(`mse'^2 + (`z'*`sse')^2)
    	    	    if `sdf'==0 {
			replace `sse'=`S'/sqrt(2*(`nobs'-1))
		        replace `c`i'se'=sqrt(`mse'^2 + (`z'*`sse')^2)
		    }
		    if "`cv'"!="" {
	 	    	if `sdf'==0 {
			   replace `sse'=sqrt(`S'^2/(2*(`nobs'-1)*`M'^2) /*
			    */ +(`mse'*`S'/`M'^2)^2)
			}
 			replace `c`i'se'=sqrt( (`z'*`M'*`sse')^2 /*
			 */ +((1+(`z'*`S'))*`mse')^2 +(`z'*`mse'*`sse')^2)
		    }
		}
	}
	}
/*
	Graph
*/
if "`graph'"!="nograph" {
	if "`saving'"!="" { local saving "saving(`saving')" }
	if "`mbase'`sbase'`gbase'`dbase'"!="" { loc cs "s(.pppppppp)" }
	else loc cs "c(.ssssssss) s(oiiiiiiii)"
	loc title "Normal model, `centile' centiles"
	graph `y' `M' `C' `rhs', title("`title'") /*
	 */ l1("`lhs'") b2("`rhs'") xlab ylab `cs' sort `saving'
}
/*
	Save variables if required
*/
/* need to save weights first */
tempvar wts
if "`cv'"!="" { qui gen `wts'=`wt'*`M'^-2 }
else { qui gen `wts'=`wt' } 
if "`leave'"!="noleave" {

	local prog "_gls"

	local vn Z`prog'
	cap drop `vn'
	rename `Z' `vn'
	lab var `vn' "`dupper' Z-scores"

	if "`cv'"!="" { local cvl "CV" }
	else		local cvl "SD"
	if "`lns'"!="" {
		qui replace `S' = ln(`S')
		local lnsl "log "
	}

	local a1 m
	local a2 s

	local b1 "mean"
	local b2 "`lnsl'`cvl'"

	local i 1
	local p 2
	while `i'<=`p' {
		local a `a`i''
		local A = upper("`a'")
		local b `b`i''
		local vn `A'`prog'
		cap drop `vn'
		rename ``A'' `vn'

		if ``a'df'==0	{ local add "constant" }
		else	 	{ local add "powers ``a'pow'" }
		lab var `vn' "`dupper' `b': `add'"
		local i = `i'+1
	}
	loc i 0
	while `i'<`nc' {
		loc i=`i'+1
		loc cent `cent`i''
		loc vn=substr("C`cent'`prog'",1,8)
		cap drop `vn'
		rename __C`cent' `vn'
		lab var `vn' "`dupper' `cent' centile"
		if "`se'"!="" {
			loc vn=substr("se`vn'",1,8)
			cap drop `vn'
			rename `c`i'se' `vn'
			lab var `vn' "`dupper' se[`cent' centile]"
		}
	}
}
else {
	local i 0
	while `i'<`nc' {
		local i=`i'+1
		local cent `cent`i''
		cap drop __C`cent'
	}
}

di in gr _n "Final deviance = " %9.3f in ye `dev' /*
 */ in gr " (" in ye `obs' in gr " observations)."
di in gr "Power(s) for mean curve = " in ye "`mpow'" in gr "." _cont
if trim("`spow'")!="" {
	di in gr " Power(s) for " "`slab'" " curve = " in ye "`spow'" in gr "." _cont
}
di
if "`detail'"!="" {
    if trim("`mpow'")!="" {
    	di in gr _n "Regression for mean curve" _n _dup(25) "-"
	regress `y' `mvars' `mbase' [aweight=`wts'], `mopt' depname(`lhs')
        describe `mvars'
    }
    if trim("`spow'")!="" {
	di in gr _n "Regression for " "`slab'" " curve" _n _dup(23) "-"
	regress `absres' `svars' `sbase' `e', `sopt' depname("Abs. res")
        describe `svars'
    }
}
if "`tidy'"!="notidy" {
	if "`mvars'"!="" & "`mvars'"!="`rhs'" { drop `mvars' }
	if "`svars'"!="" & "`svars'"!="`rhs'" { drop `svars' }
}
global S_1 `dev'
global S_2 `mpow'
global S_3 `spow'
end

program define _devian /* deviance for xrigls */
	local wt `1'
	local obs `2'
	local ssr `3'
	local M `4'
/*
	Mean log normalized weights - use means.ado
	Note: local mnlnwt = log($S_4/$S_3) --- log (geom mean/arith mean)

	S_1 = deviance
*/
	tempvar w
	qui gen `w'=`wt'*`M'^-2
	means `w'
	global S_1 = `obs'*(1+log(2*_pi*`ssr'/`obs')-log($S_4/$S_3))
end


*! V 1.0.2 PR  12Mar97.
program define frachoos
loc 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) noOMit noHEad DF(int 0) BAse(string)
 Name(string) POwers(string) SCAling noTIDy *" ;
#delimit cr
parse "`*'"
if `df'<0 {
	di in red "invalid df"
	exit 198
}
if `df'==0 { local df 4 }
if `df'==1 {					/* linear */
	if "`powers'"!="" {
		di in red "powers() invalid with df(1)"
		exit 198
	}
}
else {
	if "`powers'"=="" { local powers "-2,-1,-.5,0,.5,1,2,3" }
	local powers "powers(`powers')"
	local degree=int(`df'/2+.5)
	local deg "degree(`degree')"
}
if "`name'"!="" { local name "name(`name')" }
if "`weight'"!="" { 
	local weight "[`weight'`exp']"
}
parse "`varlist'", parse(" ")
local lhs `1'
mac shift
local rhs
if "`head'"=="" {
	di in gr "Variable" _col(11) "df" _col(21) "Deviance" /*
	 */ _col(33) "Deviance" _col(47) "P" /*
	 */ _col(55) "Final"
	di in gr _col(32) "difference" _col(55) "Powers     df"
	di in gr _dup(67) "-" _cont
}
local i 1
while "``i''"!="" { 
	local n ``i''
	qui fracpoly `cmd' `lhs' `n' `base' `if' `in' /*
	 */ `weight', compare `deg' `powers' `options'
	drop $S_E_xp
	cap drop _sample
	local devdeg ${S_E_d`degree'}
	local vname "``i''"
	local degx `degree'
	local done 0
	while `degx'>0 & !`done' {
		local degx1 `degx'
		local degx=`degx'-1
		local df=2*`degx1'
		if `degx'==0 {
			local dev $S_E_dlin
		}
		else {
			local dev ${S_E_d`degx'}
		}
		local d=`dev'-`devdeg'
		local P ${S_E_P`degx1'}
		local pwrs ${S_E_p`degx1'}
		di _n in gr "`vname'" in ye _col(12) `df' /*
		 */ _col(20) %9.3f `devdeg' /*
		 */ _col(32) %9.3f `d' /*
		 */ _col(44) %6.3f `P' _col(55) "`pwrs'" _col(67) _cont
		local vname
		if `P'<`alpha' | "`omit'"!="" {
			local done 1
			local dev `devdeg'
			qui fracgen `n' `pwrs', replace `scaling' `name'
			local n "$S_1"
		}
		local devdeg `dev'
	}
/*
	Deal with linear.
*/
	if !`done' {
		local devdeg $S_E_dlin
		local P $S_E_Plin
		local d=$S_E_d0-`devdeg'
		local df 1
		local pwrs 1
		di _n in gr "`vname'" in ye _col(12) `df' /*
		 */ _col(20) %9.3f `devdeg' /*
		 */ _col(32) %9.3f `d' /*
		 */ _col(44) %6.3f `P' _col(55) "`pwrs'" _col(67) _cont
		if "`omit'"=="" & `P'>`alpha' {
/*
	`Dropping' RHS variable since 1 df test of beta=0 is non-sig.
*/
			local done 1
			local dev $S_E_d0
			local df 0
			local pwrs
			local n
		}
		if !`done' {
			local dev `devdeg'
		}
	}
	di `df'
	local rhs "`rhs' `n'"
	local i = `i'+1
}
di in gr _dup(67) "-"
global S_1 `rhs'
global S_2 `df'
global S_3 `dev'
global S_4 `pwrs'
global S_5 `n'
end

*! version 1.3.3 PR 02Dec96.
program define fracpoly
version 4.0
if "`1'" == "" | substr("`1'",1,1)=="," {
	if "$S_E_cmd"=="" 		{ error 301 }
	if "$S_E_fp"!="fracpoly"	{ error 301 }
	di in blue "->$S_E_cmd"
	$S_E_cmd
	_fraccmp `*'
	exit
}
loc cmd `1'
mac shift
if substr("`cmd'",1,3)=="reg" { local cmd regress }
if "`cmd'"!="regress" & "`cmd'"!="fit" & "`cmd'"!="qreg" & "`cmd'"!="rreg" /*
 */ & "`cmd'"!="poisson" & "`cmd'"!="glm" & "`cmd'"!="logit" & "`cmd'"!="clogit" /*
 */ & "`cmd'"!="logistic" & "`cmd'"!="cox" & "`cmd'"!="xtgee" {
	di in red "invalid or unrecognised command, `cmd'"
	exit 198
}
local glm="`cmd'"=="glm"
local qreg="`cmd'"=="qreg"
local xtgee="`cmd'"=="xtgee"
/*
    dist=0 (normal), 1 (binomial), 2 (poisson), 3 (cox), 4 (glm), 5 (xtgee).
*/
local dist=("`cmd'"=="logit"|"`cmd'"=="clogit"|"`cmd'"=="logistic") /*
*/ + 2*("`cmd'"=="poisson") + 3*("`cmd'"=="cox") + 4*`glm'+ 5*`xtgee'
local normal=("`cmd'"=="regress")|("`cmd'"=="fit")|("`cmd'"=="rreg")

loc star "`1' `2'"
mac shift 2
/*
	Look for fixpowers
*/
loc done 0
while !`done' {
	local s=index("`1'",",")
	local doshift=`s'==0	/* no comma found */
	if `doshift' {
		local pow `1'
	}
	else {
		local pow=substr("`1'",1,`s'-1)
	}
	cap confirm num `pow'
	if _rc==0 {	/* pow confirmed as a number */
		if "`fixpowe'"=="" { loc fixpowe `pow' }
		else loc fixpowe "`fixpowe' `pow'"
		local 1=substr("`1'",`s',.)
		if `doshift' { mac shift }
	}
	else loc done 1
}
local star "`star' `*'"

loc varlist "req ex min(2)"
loc if "opt"
loc in "opt"
loc weight "aweight fweight pweight iweight"
loc search="`fixpowe'"==""
if `search' {
	loc srchopt "ADDpowers(str) DEVthr(real 1e30) POwers(str) DEGree(int 2) LOg COMpare"
}
loc options "EXPx(str) ORIgin(str) ZERo noCONStant noSCAling `srchopt' *"
#delimit cr
loc small 1e-6
parse "`star'"
if `search' {
	if `degree'<1 | `degree'>9 {
		di in red "invalid degree()"
		exit 198
	}
	loc df=2*`degree'
	loc odddf=(2*int(`df'/2)<`df')
}
else loc df 1
loc lin="`fixpowe'"=="1"
if "`constan'"=="noconstant" {
	if "`cmd'"=="fit" | "`cmd'"=="cox" {
		di in red "noconstant invalid with `cmd'"
		exit 198
	}
	loc options "`options' nocons"
}
/*
	 Read powers to be searched into a variable, `p'
*/
if `search' {
	if "`powers'"=="" { loc powers "-2,-1,-.5,0,.5,1,2,3" }
	loc pwrs "`powers' `addpowe'"
	tempvar p
	qui gen `p'=.
	_fracpp `p' "`pwrs'" 1 1
	loc np $S_1
	if `np'<1 {
		di in red "no powers given"
		exit 2002
	}
	loc pwrs "$S_3"
}
gl S_E_cmd
parse "`varlist'", parse(" ")
loc y `1'
mac shift
loc rhs `1'
mac shift
loc base `*'
tempvar touse x lnx

quietly {
	mark `touse' `if' `in'
	markout `touse' `rhs' `y' `base'
	gen `x'=`rhs' if `touse'
	lab var `touse' "fracpoly sample"
/*
	Deal with weights.
*/
	loc mnlnwt 0 /* mean log normalized weights */
	if "`exp'" != "" { 
		tempvar s
		gen `s' `exp'
		replace `x'=. if `s'<=0|`s'==.
		replace `s'=. if `x'==.
		replace `touse'=0 if `x'==.
		if "`weight'"=="fweight" { 
			capture assert `s'==int(`s') if `touse'
			if _rc { noi error 401 } 
		}
		sum `s' 
		noi di "(sum of wgt is " _result(1)*_result(3) ")"
		replace `s'=log(`s'/_result(3))
		sum `s'
		loc mnlnwt=_result(3)
		drop `s'
		loc wgt "[`weight'`exp']"
	}
	gen `lnx'=.
	if `lin' { loc zero }	
	_fracxo `x' `lnx' `lin' "`expx'" "`origin'" "`zero'" "`scaling'"
	loc nobs $S_1
	if $S_3==0 {
		loc zeta $S_2
		loc shift 0
	}
	else loc shift $S_2
	loc kx "$S_4"
	loc scalfac $S_5
}
/*
	Determine residual and model df.
*/
qui reg `y' `base' `wgt' if `touse'
loc rdf=_result(5)+("`constan'"=="noconstant")
if `rdf'<(1+`df') { error 2001 }
/*
	Calc deviance=-2(log likelihood) for regression on covars only,
	allowing for possible weights.

	Note that for logit/clogit/logistic with nocons, must regress on zero,
	otherwise get r(102) error.
*/
if `dist'==1 & "`constan'"=="noconstant" {
	tempvar z0
	qui gen `z0'=0
}
qui `cmd' `y' `z0' `base' `wgt' if `touse', `options'
if `xtgee' & "`base'"=="" { global S_E_chi2 0 }
loc glmwarn 0
if `glm' {
 if abs($S_E_dc/$S_E_del-1)>`small' & abs($S_E_dd/$S_E_del-1)>`small' {
	loc scale $S_E_del
 }
 else loc glmwarn 1
}
_fracdv `normal' "`wgt'" `nobs' `mnlnwt' `dist' `glm' `xtgee' `qreg' "`scale'"
loc dev0 $S_1
if `normal' { loc rsd0=_result(9) }
if `search' {
	loc m `degree'
	if "`log'"!="" {
		di _n in gr "Model #      Deviance " _cont
		if `normal' { di in gr "  Res S.D. " _cont }
		loc i 1
		while `i'<=`m' {
			di in gr " Power " `i' _cont
			loc i=`i'+1
		}
		di _n
	}
/*
	Go!
*/
	loc i 1
	while `i'<=`np' {
		loc pi=`p'[`i']
		tempvar x`i'
		qui gen `x`i''=cond(`pi'==0, `lnx', /*
		 */ cond(`pi'==1, `x', cond(`x'==0,0,`x'^`pi')))
		loc i=`i'+1
	}
	loc i 0
	while `i'<`m' {
		loc j`i' 0
		loc i=`i'+1
		tempvar xp`i'
		qui gen `xp`i''=.
		loc j=`i'+`i'-1
		loc dev`j' `dev0' /* min deviance for df=j */
		loc j=`j'+1
		loc dev`j' `dev0'
	}
	loc j`m' 1
	loc i `m'
	loc kount 0  /* no. of models processed */
	loc kountd 0 /* no. of models with deviance<`devthr' */
	loc deg `j`m''
	loc done 0
	while !`done' {
		if `i'<`m' {
			loc ji `j`i''
			qui replace `xp`i''=`x`ji''
			loc l `i'
			while `l'<`m' {
				loc l=`l'+1
				loc j`l' `ji'
				loc l1=`l'-1
				qui replace `xp`l''=`xp`l1''*`lnx'
			}
			if "`log'"=="" { di "." _cont }
		}
		else qui replace `xp`m''=`x`j`m'''
/*
	Test for any power being 1 (i.e. odd-df model)
*/
		loc one 0
		loc l 0
		while `l'<`m' {
			loc l=`l'+1
			if `p'[`j`l'']==1 { loc one 1 }
		}
		if !`odddf' | (`deg'<`m') | `one' {
			loc dfi=2*`deg'-`one'
			loc i 1
			loc xlist
			while `i'<=`m' {
				if `j`i''>0 { loc xlist "`xlist' `xp`i''" }
				loc i=`i'+1
			}
			loc kount=`kount'+1
			qui `cmd' `y' `xlist' `base' `wgt' if `touse',  /*
			 */ `options'
			_fracdv `normal' "`wgt'" `nobs' `mnlnwt' `dist' /*
			 */ `glm' `xtgee' `qreg' "`scale'"
			loc dev $S_1
			if `normal' { loc rsd=_result(9) }
			if `dev'<`dev`dfi'' {
				loc dev`dfi' `dev'
				if `normal' { loc rsd`dfi' `rsd' }
				loc i 1
				loc pm`dfi'
				while `i'<=`m' {
					if `j`i''>0 {
						loc pm=`p'[`j`i'']
						loc pm`dfi' "`pm`dfi'' `pm'"
					}
					loc i=`i'+1
				}
			}
			if `dev'<`devthr' { loc kountd=`kountd'+1 }
			if "`log'"!="" {
				di in ye %5.0f `kount' /*
				 */ _col(10) %12.3f `dev' _cont
				if `normal' { di "   " %8.0g `rsd' _cont }
				loc i `m'
				while `i'>0 {
					if `j`i''==0 { di _skip(6) ". " _cont }
					else di in ye %8.1f `p'[`j`i''] _cont
					loc i=`i'-1
				}
				di
			}
		}
/*
	Increment the first possible index (of loop i) among indices of
	loops m, m-1, m-2, ..., 1
*/
		loc i `m'
/*
	Finish after all indices have achieved their upper limits (np).
*/
		while `j`i''==`np' {
			loc i=`i'-1
		}
		if `i'==0 { loc done 1 }
		else {
			if `j`i''==0 { loc deg = `m'-`i'+1 }
			loc j`i'=`j`i''+1
		}
	}
/*
	Update the results for even df to include odd df
*/
	loc i 2
	while `i'<=`df' {
		loc j=`i'-1
		if `dev`j''<`dev`i'' {
			loc dev`i' `dev`j''
			if `normal' { loc rsd`i' `rsd`j'' }
			loc pm`i' "`pm`j''"
		}
		loc i=`i'+2
	}
	if "`log'"=="" { di }
}
/*
	Store coefficients for linear fit on untransformed x
*/
qui `cmd' `y' `rhs' `base' `wgt' if `touse', `options'
cap local b0=_b[_cons]
if _rc { local b0 0 }
local b1=_b[`rhs']
/*
	Create FP transformation(s) of xvar for final model
*/
if `search' { loc pwrs `pm`df'' }
else loc pwrs `fixpowe'
if "`expx'"!=""   { local e "expx(`expx')" }
if "`origin'"!="" { local o "origin(`origin')" }
if !`lin' | (`lin' & "`e'`o'"!="") {
	fracgen `rhs' `pwrs', replace `zero' `scaling' `e' `o'
	loc xp $S_1
}
else loc xp `rhs'
if "`zeta'"!="" {
	version 2.1
	local zeta = `zeta'	/* round */
	version 4.0
	if `shift'!=0 {
		if `zeta'<0 {
			local z "-"
			local zeta = -`zeta'
		}
		else local z "+"
	}
	else {
		if `zeta'>0 { local z "-" }
		else {
			local zeta = -`zeta'
			local z "+"
		}
		di in bl "`rhs' -> (`rhs'`z'`zeta')/(max(`rhs')`z'`zeta')."
	}
}
/*
	Fit final model with permanent _sample filter
*/
cap drop _sample
rename `touse' _sample
`cmd' `y' `xp' `base' `wgt' if _sample, `options'
if !`search' {
/*
	Deviance for fixed-powers model.
	Note that `dev1' is stored in $S_E_dlin, deviance for a
	linear model, even when `fixpowe' is not 1; similarly 
	for `rsd1'.
*/
	_fracdv `normal' "`wgt'" `nobs' `mnlnwt' `dist' `glm' /*
	 */ `xtgee' `qreg' "`scale'"
	loc dev1 $S_1
	if `normal' { loc rsd1=_result(9) }
	loc pm1 `fixpowe'
}
di in gr "Deviance:" in ye %9.3f `dev`df'' in gr ". " _cont
if `search' {
	di in gr "Best powers of " in ye "`rhs'" in gr " among " /*
	 */ in ye `kount' in gr " models fit: " in ye "`pwrs'" _cont
	if `devthr'<1e30 { di _n in gr "Number of models with deviance < " /*
	 */ in ye `devthr' in gr ": " in ye `kountd' _cont }
	di in gr "." _cont
}
di
global S_E_d0 `dev0'
global S_E_s0 `rsd0'
global S_E_dlin `dev1'
if `normal' { global S_E_slin `rsd1' }
loc i 2
loc j 1
while `i'<=`df' {
	global S_E_d`j' `dev`i''
	global S_E_p`j' `pm`i''
	if `normal' { global S_E_s`j' `rsd`i'' }
	loc i=`i'+2
	loc j=`j'+1
}
global S_E_b0 `b0'
global S_E_b1 `b1'
global S_E_base `base'
global S_E_df `df'
global S_E_depv `y'
global S_E_dist `dist'
global S_E_wgt "`weight'"
global S_E_exp "`exp'"
global S_E_fprp no
global S_E_fp fracpoly
global S_E_glmw `glmwarn'
global S_E_nobs `nobs'
global S_E_nxp 0	/* redundant */
global S_E_pwrs `pwrs'
global S_E_rdf `rdf'
global S_E_rhs `rhs'
global S_E_srch `search'
global S_E_sfac `scalfac'
global S_E_xp `xp'
global S_E_xpx `kx'
global S_E_cmd `cmd'
_fraccmp, `compare'
end

program define _fracdv /* deviance calculation */
* 1=normal (1=yes), 2=weight type, 3=nobs, 4=mnlnwt, 5=dist, 6=glm (1=yes),
* 7=xtgee (1=yes), 8=qreg (1=yes), 9=scale (for glm).
* Returns deviance in $S_1.

local normal "`1'"
local e "`2'"
local nobs "`3'"
local mnlnwt "`4'"
local dist "`5'"
local glm "`6'"
local xtgee "`7'"
local qreg "`8'"
local scale "`9'"

if `normal' {
	if substr("`e'",1,2)=="iw" {
		global S_1=(`nobs'*(log(2*_pi)-`mnlnwt')+_result(4))
	}
	else {
		global S_1=`nobs'*(1+log(2*_pi*_result(4)/`nobs')-`mnlnwt')
	}
}
else {
	if `dist'==2 { 
		global S_1=-2*$S_E_ll
	}
	else if `glm' {
		global S_1=$S_3/`scale'
	}
	else if `xtgee' {
		global S_1=-$S_E_chi2
	}
	else if `qreg' {
		global S_1=`nobs'*log($S_E_msd/$S_E_rsd)
	}
	else global S_1=-2*_result(2)
}
end

program define _fracpv /* P-value calculation */
* 1=normal (1=yes), 2=weight type, 3=nobs, 4=d (deviance),
* 5=n1 (df for numerator), 6=n2 (df for denominator).
* Returns P-value in $S_1.

local normal "`1'"
local e "`2'"
local nobs "`3'"
local d "`4'"
local n1 "`5'"
local n2 "`6'"

if (`normal' & substr("`e'",1,2)=="iw") | !`normal' {
	global S_1=chiprob(`n1',`d')
}
else {
	global S_1=fprob(`n1',`n2',(exp(`d'/`nobs')-1)*`n2'/`n1')
}
end

program define _fraccmp /* report model comparisons */
local options "COMpare"
parse "`*'"
if "`compare'"=="" | $S_E_srch==0 { exit }
local normal=($S_E_dist==0)
local dash "     --"
local ddup=54+15*`normal'
di in gr _n "Fractional Polynomial Model Comparisons:"
di in gr _dup(`ddup') "-"
di in gr "$S_E_rhs" _col(18) "Deviance     " _cont
if `normal' { di in gr " Res. SD     " _cont }
di in gr " Gain      P      Powers" 
di in gr _dup(`ddup') "-"

di in gr "Not in model" in ye _col(13) %13.3f $S_E_d0 _cont
if `normal' { di in ye _skip(5) %8.0g $S_E_s0 _cont }
di in gr _skip(3) "`dash'`dash'"
local i 1
while `i'<=$S_E_df {
	local m=1+int((`i'-1)/2)
	if `i'==1 {
		di in gr "Linear" _col(13) _cont
		local dev $S_E_dlin
		local rsd $S_E_slin
		local d=$S_E_d0-`dev'
		local n1 1
		local n2=$S_E_rdf-1
		local pwr 1
	}
	else {
		di in gr "m = `m' (`i' df)" _col(13) _cont
		local dev ${S_E_d`m'}
		local rsd ${S_E_s`m'}
		local d=`devlast'-`dev'
		if `m'==1 { local n1 1 }
		else 	    local n1 2
		local n2=$S_E_rdf-`n1'
		local pwr ${S_E_p`m'}
	}
	_fracpv `normal' "$S_E_wgt" $S_E_nobs `d' `n1' `n2'
	local P $S_1
	if `i'==1 { global S_E_Plin `P' }
	else global S_E_P`m' `P'
	di in ye %13.3f `dev' _cont
	if `normal' {
		di in ye _skip(5) %8.0g `rsd' _cont
	}
	di in ye _skip(3) %7.3f $S_E_dlin-`dev' /*
	 */ %9.3f `P' _skip(2) "`pwr'"
	local devlast `dev'
	if `i'==1 { local i=`i'+1 }
	else 	    local i=`i'+2
}
di in gr _dup(`ddup') "-"
if $S_E_glmw { di in bl /*
 */ "[warning: GLM deviances are unscaled---ignore P-values]" }
end

program define _jprxrpa /* parses "clusters" */
local stuff "`1'"	/* string to be parsed */
local prefixs "`2'"	/* string of permitted one-character prefixes */
local items "`3'"	/* items in error message if #clusters > #prefixes */

local nprefix=length("`prefixs'")
local i 1
while `i'<=`nprefix' {
	local p`i'=substr("`prefixs'",`i',1)
	global S_p`i'
	global S_c`i'
	local i=`i'+1
}
parse "`stuff'", parse(",")
global S_1 0 /* # of comma-delimited clusters */
while "`1'"!="" {
	if "`1'"=="," { mac shift }
	global S_1 = $S_1+1
	local clust$S_1 "`1'"
	mac shift
}
if "`clust${S_1}'"=="" { global S_1=$S_1-1 }
if $S_1>`nprefix' {
	noi di in red "too many `items' specified"
	exit 198
}
/*
	Disentangle each prefix:string cluster
*/
local i 1
while `i'<=$S_1 {
	parse "`clust`i''", parse("=:")
	local prefix = lower(substr("`1'",1,1))
	local j 1
	local found 0
	while (`j'<=`nprefix') & !`found' {
		if "`prefix'"=="`p`j''" {
			global S_p`i' `prefix'
			global S_c`i' "`3'"
			local found 1
		}
		local j=`j'+1
	}
	if !`found' {
		noi di in red "invalid `clust`i''"
		exit 198
	}
	local i = `i'+1
}
end