*! version 1.4.0 PR 06Mar1999.    STB-49 sg112, STB-50 sg112.1
program define boxtid, eclass
version 6
if replay() {
	if "`e(cmd)'"=="" | "`e(fp_cmd2)'"!="boxtid" {
		error 301
	}
	di in blue "->`e(cmd)'"
	_rslt
	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] [, /* 
 */ ADJust(str) ALL noCONStant DEAD(varname) DF(str) DFDefault(int 2) /*
 */ INit(str) ITer(int 100) LTOLerance(real .001) EXPon(varlist) /*
 */ POWers(numlist) TRace ZERo(varlist) * ]
frac_cox "`dead'" `dist'
if "`trace'"!="" { local trace noi }
if "`constan'"=="noconstant" {
	if "`cmd'"=="fit" | "`cmd'"=="cox" {
		di in red "noconstant invalid with `cmd'"
		exit 198
	}
	if "`adjust'"=="" {
		local adjust no
		di in bl "[Note: adjust(no) assumed with `constan']"
	}
	else {
		di in bl "[Note: adjust(no) advised with `constan']"
	}
	local options "`options' nocons"
}
if "`powers'"=="" { local powers 2 -1 -0.5 0 0.5 1 2 3 }
tempname small
scalar `small'=1e-6
gettoken y xvars : varlist
tempvar touse
quietly {
	marksample touse
	markout `touse' `y' `xvars' `dead'
/*
	Deal with weights.
*/
	frac_wgt `"`exp'"' `touse' `"`weight'"'
	local mnlnwt = r(mnlnwt) /* mean log normalized weights */
	local wgt `r(wgt)'
	if "`dead'"!="" { local dead "dead(`dead')" }

	local nxvars: word count `xvars'
/*
	Adjustment
*/
	frac_adj "`adjust'" "`xvars'" `touse'
	local i 1
	while `i'<=`nxvars' {
		if "`r(adj`i')'"!="" {
			local Adj`i' adjust(`r(adj`i')')
		}
		local uniq`i'=r(uniq`i')
		local i=`i'+1
	}
/*
	Degrees of freedom
*/
	if "`df'"!="" {
		frac_dis "`df'" df 1 . "`xvars'"
		local i 1
		while `i'<=`nxvars' {
			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.

	Also, initialisations for vars in non-linear part of model.
*/
	local rhs
	local covars
	local nx 0
	local nbase 0
	local i 1
	tokenize `xvars'
	while "``i''"!="" {
		frac_mun ``i'' purge
		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' }
		}
		if `df`i''==1 {
			local nbase=`nbase'+1
			local covars `covars' ``i''
		}
		else {
			local nx=`nx'+1
			local rhs `rhs' ``i''
			local v`nx' ``i''
			local deg`nx'=int(.01+.5*`df`i'')
			local xnz`nx' `touse'
			local ifxnz`nx' "if `touse'"
			local lnx`nx' 1
		}
		local i=`i'+1
	}
/*
	Assign adjustment macros adj for covars (df=1) and rhs (df>1).
	covars follow rhs vars in adj`i' macros.
*/
	local ix 0
	local ic `nx'
	local i 1
	while `i'<=`nxvars' {
		if `df`i''==1 {
			local ic=`ic'+1
			local adj`ic' `Adj`i''
		}
		else {
			local ix=`ix'+1
			local adj`ix' `Adj`i''
		}
		local Adj`i'
		local i=`i'+1
	}
/*
	Vars with expon option
*/
	if "`expon'"!="" {
		tokenize `expon'
		while "`1'"!="" {
			frac_in `1' "`rhs'"
			local lnx`s(k)' 0
			mac shift
		}
	}
/*
	Vars with zero option
*/
	if "`zero'"!="" {
		tokenize `zero'
		while "`1'"!="" {
			frac_in `1' "`rhs'"
			local j `s(k)'
			local zero`j' "zero"
			tempvar xnz`j'
			gen byte `xnz`j''=(`v`j''>0 & `touse'==1)
			local ifxnz`j' "if `xnz`j''"
			mac shift
		}
	}
/*
	Put initial values into init`v'.
	(Note: correct for expx(sd) transformation in estimate of power)
*/
	if "`init'"!="" {
		frac_ext "`rhs'" "`init'" init
		local v 0
		while `v'<`nx' {
			local v=`v'+1
			local init`v' ${S_`v'}
		}
	}
	local v 0
	while `v'<`nx' {
		local v=`v'+1
		local adj
		local l 0	
		while `l'<`nx' {
			local l=`l'+1
			if `l'!=`v' { local adj `adj' `v`l'' }
		}
		`trace' init `cmd' `y' `v`v'' `lnx`v'' "`adj'" /*
		 */ "`covars'" "`wgt'" `touse' `deg`v'' "`dead'" /*
		 */ "`zero`v''" "`powers'" "`options'" "`init`v''"
		local pwrs `r(pwrs)'
		local f=r(f)
		initial `deg`v'' `f' "`pwrs'"
		local i 1
		while `i'<=`deg`v'' {
			local init`v'`i' ${S_`i'}
			local i=`i'+1
		}
	}
	count if `touse'
	local nobs=r(N)
/*
	Determine residual and model df.
*/
	regress `y' `covars' `wgt' if `touse'
	local rdf=e(df_r)+("`constan'"=="noconstant")
/*
	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 (`glm' | `dist'==1) & "`constan'"=="noconstant" {
		tempvar z0
		gen `z0'=0
	}
	`cmd' `y' `z0' `covars' `wgt' if `touse', `dead' `options'
	if `xtgee' & "`covars'"=="" { global S_E_chi2 0 }
	local glmwarn 0
	if `glm' {
	 local scale 1
	 if abs($S_E_dc/$S_E_del-1)>`small' & abs($S_E_dd/$S_E_del-1)>`small' {
		local scale $S_E_del
	 }
	 else local glmwarn 1
	}
	frac_dv `normal' "`wgt'" `nobs' `mnlnwt' `dist' /*
	 */ `glm' `xtgee' `qreg' "`scale'"
	local dev0=r(deviance)
/*
	Initial fit for each predictor
*/
	local wvars
	local v 0
	while `v'<`nx' {
		local v=`v'+1
		tempvar x`v'
		tempname range`v' tmp
		if `lnx`v'' {
                        fracgen `v`v'' 0 if `touse', name(`tmp') `zero`v'' /*
			 */ noscaling
                        rename `r(names)' `x`v''
		}
		else gen `x`v''=`v`v'' if `touse'
/*
	Standardize useable values of x`v' to [0,1].
*/
		sum `x`v'' `ifxnz`v''
		scalar `range`v''=r(max)-r(min)
		replace `x`v''=(`x`v''-r(mean))/`range`v'' `ifxnz`v''
		local i 1
		while `i'<=`deg`v'' {
			tempname g`v'`i'
			scalar `g`v'`i''=`init`v'`i''*`range`v''
			tempvar w`v'`i' z`v'`i'
			gen double `w`v'`i''= /*
			 */ cond(`xnz`v'',exp(`g`v'`i''*`x`v''),0) if `touse'
			sum `w`v'`i''
			if r(min)==r(max) | r(N)<`nobs' {
	noi di in red "boxtid fit failed, derived variable has zero variance"
				exit 498
			}
			gen double `z`v'`i''=.
			local wvars `wvars' `w`v'`i''
			local i=`i'+1
		}
	}
/*
	Initial fit of whole model
*/
	`cmd' `y' `wvars' `covars' `wgt', `dead' `options'
	frac_dv `normal' "`wgt'" `nobs' `mnlnwt' `dist' /*
	 */ `glm' `xtgee' `qreg' "`scale'"
	local devnow=r(deviance)
	noi di in gr _n "Iteration 0:  Deviance = " in ye %9.0g `devnow'
	local v 0
	while `v'<`nx' {
		local v=`v'+1
		local i 1
		while `i'<=`deg`v'' {
			tempname b`v'`i'
			scalar `b`v'`i''=_b[`w`v'`i'']
			local i=`i'+1
		}
	}
/*
	Iterate
*/
	tempname gold ratio bz

	local devinit `devnow'
	local devprev `devnow'
	local devcon `ltolera' /* deviance convergence criterion */
	local converg 1
	local totalit 0
	local done 0
	local j 1
	while `j'<=`iter' & !`done' {
	    local v 0
	    while `v'<`nx' & !`done' {
		    local v=`v'+1
/*
	First loop: update g2 conditional on g1
	(g1 uses w1 & z1, g2 uses w2 & z2)
*/
		    local i 1
		    while `i'<=`deg`v'' {
			local devpr2 `devnow'
			replace `z`v'`i''=`w`v'`i''*`x`v''
			`cmd' `y' `wvars' `z`v'`i'' `covars' `wgt', /*
			 */ `dead' `options'
			scalar `gold'=`g`v'`i''
			cap scalar `bz'=_b[`z`v'`i'']
			if _rc { scalar `bz'=`b`v'`i'' }
			scalar `ratio'=`bz'/`b`v'`i''
			local dun 0
			local s 1
			while !`dun' & `s'>-2 {
				if `s'<0 { noi di in gr "(step sign changed)" }
				local gf 1
				while !`dun' & `gf'<=1e6 {
					scalar `g`v'`i''=`gold'+`s'*`ratio'/`gf'
					replace `w`v'`i''=cond(`xnz`v'', /*
					 */ exp(`g`v'`i''*`x`v''),0) if `touse'
					count if `w`v'`i''!=.
					if _result(1)<`nobs' {
						local gf=`gf'*10
						noi di in gr "(invalid transformation attempt, step length divided by 10)"
					}
					else {
						`cmd' `y' `wvars' `covars' `wgt', /*
						 */ `dead' `options'
						frac_dv `normal' "`wgt'" `nobs' `mnlnwt' /*
						 */ `dist' `glm' `xtgee' `qreg' "`scale'"
						local devdiff=(r(deviance)-`devpr2')
						if `devdiff'<0 { /* accept new parameters */
							local devnow=r(deviance)
							local i2 1
							while `i2'<=`deg`v'' {
								scalar `b`v'`i2''=_b[`w`v'`i2'']
								local i2=`i2'+1
							}
							local dun 1
						}
						else {
							local gf=`gf'*10
							noi di in gr "(unprofitable step attempted, step length divided by 10)"
						}
					}
				}
				local s=`s'-2
			}
			if "`trace'"!="" {
			    noi di in gr "`v`v''(`i')" %4.0f `j2' /*
			     */ _col(16) in ye %9.0g /*
			     */ `g`v'`i''/`range`v'' /*
			     */ _col(28) %12.6f r(deviance) /*
			     */ _col(44) %9.6f `devdiff'
			}
			local devpr2 `devnow'
			local totalit=`totalit'+1
			local i=`i'+1
		    }	/* i loop (term within variable) */
	    }
	    if `nx'>1 { local devdiff=(`devnow'-`devprev') }
	    if "`devdiff'"=="" { local devdiff 0 }
	    if `devdiff'>10 {
		local done 1
		local converg 0
	    }
	    else if `devdiff'>-`devcon' {
		local done 1
	    }
	    if "`trace'"!="" { noi di }
	    noi di in gr "Iteration `j':  Deviance = " in ye %9.0g `devnow' /*
	     */ in gr " (change = " in ye %9.0g `devdiff' in gr ")"
	    local devprev `devnow'
	    local j=`j'+1
	}
	if !`done' | (`devnow'>`devinit') { local converg 0 }
/*
	Tests for linearity for each nonlinear term
	(could make this an option)
*/
	local terms: word count `wvars'
	local mdf=2*`terms'	/* model df for power terms + Taylor terms */
	local v 0
	while `v'<`nx' {
		local v=`v'+1
/*
	Replace wv with power terms for all but v'th predictor.
	Delete terms recursively from wvars.
*/
		local wv `wvars'
		local i 1
		while `i'<=`deg`v'' {
			frac_str `w`v'`i'' `wv'
			local wv `s(new)'
			local i=`i'+1
		}
		`cmd' `y' `v`v'' `wv' `covars' `wgt' if `touse', /*
		 */ `dead' `options'
		frac_dv `normal' "`wgt'" `nobs' `mnlnwt' `dist' /*
		 */ `glm' `xtgee' `qreg' "`scale'"
		local lin`v'=_b[`v`v'']
		local se`v'=_se[`v`v'']
		local gain`v'=r(deviance)-`devnow'
		if `gain`v''>0 {
			local d `gain`v''
			local n1=2*`deg`v''-1
			local n2=`rdf'-`mdf'
			frac_pv `normal' "`wgt'" `nobs' `d' `n1' `n2'
			local P`v'=r(P)
		}
		else {
			local gain`v' 0
			local P`v' 1
		}
	}
/*
	Final model.
	Approximation to SE(k): see Snedecor and Cochran 1967 p 470.
*/
	local zvars
	local v 0
	while `v'<`nx' {
		local v=`v'+1
		local i 1
		while `i'<=`deg`v'' {
			tempname sek`v'`i'
			scalar `sek`v'`i''=.
			local zvars `zvars' `z`v'`i''
			local i=`i'+1
		}
	}
	`cmd' `y' `wvars' `zvars' `covars' `wgt', `dead' `options'
	local xp
	local pwrs
	local rc 0
	local zvars
	local fvl
	local v 0
	while `v'<`nx' {
		local v=`v'+1
		local xp`v'
		local pwrs`v'
		local i 1
		frac_mun `v`v''
		local vn `s(name)'
		gen byte `vn'_1=.   /* reserves name `vn'_1 */
		while `i'<=`deg`v'' {
			cap scalar `sek`v'`i''=(_se[`z`v'`i'']/ /*
			 */ abs(`b`v'`i''))/`range`v''
			if _rc { scalar `sek`v'`i''=. }
/*
	Rescale auxiliary variables to give their se=se(power param)
*/
*			replace `z`v'`i''=`z`v'`i''*_se[`z`v'`i'']/`sek`v'`i''
			replace `z`v'`i''=`z`v'`i''*abs(`b`v'`i'')*`range`v''
			local k=`g`v'`i''/`range`v''
			local pwrs`v' `pwrs`v'' `k'
			if `lnx`v'' {
				frac_ddp `k' 6
				local K `r(ddp)'
				local e
			}
			else {
				local K 1
				local e "expx(`k')"
			}
			noi fracgen `v`v'' `K' `K' if e(sample),`all' `e' /*
			 */ `adj`v'' replace index(`i') `zero`v'' name(`vn')
			local xx1: word 1 of `r(names)'
			local xx2: word 2 of `r(names)'
			count if `xx1'!=.
			if _result(1)<`nobs' {
				local rc 2002
			}
			local xp`v' `xp`v'' `xx1'
/*
	Rename auxiliary ("zero") variable.
*/
			rename `xx2' `vn'p`i'
			local zp`v' `zp`v'' `vn'p`i'
			lab var `vn'p`i' "boxtid auxil p`i'(`v`v'')"
			local zvars `zvars' `vn'p`i'
			local i=`i'+1
		}
		local xp `xp' `xp`v''
		local fvl `fvl' `xp`v'' `zp`v''
		local pwrs `pwrs' `pwrs`v''
	}
}
/*
	Adjust covariates if necessary, and create final var list (`fvl')
*/
local i 0
while `i'<`nbase' {
	local i=`i'+1
	local j=`i'+`nx'
	local v: word `i' of `covars'
	if "`adj`j''"!="" {
		frac_mun `v'
		local n `s(name)'
		fracgen `v' 1 if e(sample),`all' replace noscaling `adj`j'' /*
		 */ name(`n')
		local v `r(names)'
	}
	local fvl `fvl' `v'
}
if `rc'==0 { qui `cmd' `y' `fvl' `wgt', `dead' `options' }
di in gr _n "[Total iterations: `totalit']"
/*
	New code in v 1.1.6 for consistency with mfracpol
	(Differs from code in fracpoly because of multiple
	transformed predictors--offset to nbase is `nx' not 1)
*/
local i 0
while `i'<`nbase' {
	local i=`i'+1
	local j=`i'+`nx'
	global S_E_x`j': word `i' of `covars'
	global S_E_k`j' 1
	global S_E_a`j'
	est local fp_x`j': word `i' of `covars'
	est local fp_k`j' 1
	est local fp_a`j'
}
/*
	End of new code in v 1.1.6 for consistency with mfracpol
*/

global S_1 `nobs'
global S_2 `devnow'

global S_E_xp `xp'
global S_E_pwrs `pwrs'
global S_E_x `rhs'
global S_E_dev `devnow'
global S_E_cnvg `converg'
est local fp_xp `xp'
est local fp_pwrs `pwrs'
est local fp_x `rhs'
est scalar fp_dev=`devnow'
est scalar fp_cnvg=`converg'

local v 0
while `v'<`nx' {
	local v=`v'+1
	global S_E_x`v' `v`v''		/* Double save */
	global S_E_k`v' `pwrs`v''	/* Double save */
	global S_E_g`v' `gain`v''	/* Double save */
	global S_E_P`v' `P`v''		/* Double save */
	global S_E_l`v' `lin`v''	/* Double save */
	global S_E_sl`v' `se`v''	/* Double save */
	global S_E_a`v' `deg`v''	/* Double save */

	est local fp_x`v' `v`v''	/* name of v`th predictor */
	est local fp_k`v' `pwrs`v''	/* powers for v'th predictor */
	est scalar fp_g`v'=`gain`v''	/* gain for v'th predictor */
	est scalar fp_P`v'=`P`v''	/* P-value for nonlinearity */
	est scalar fp_l`v'=`lin`v''	/* linear coefficient */
	est scalar fp_sl`v'=`se`v''	/* SE of linear coefficient */
	est scalar fp_a`v'=`deg`v''	/* #auxiliaries for this var */
	local i 1
	local S
	while `i'<=`deg`v'' {
		local s=`sek`v'`i''
		local S `S' `s'
		local i=`i'+1
	}
	global S_E_s`v' `S'		/* Double save */
	est local fp_s`v' `S'		/* se(powers) for v'th predictor */
}
global S_E_base `covars'		/* Double save */
global S_E_depv `y'			/* Double save */
global S_E_dist `dist'			/* Double save */
global S_E_fvl `fvl'			/* Double save */
global S_E_wgt "`weight'"		/* Double save */
global S_E_exp "`exp'"			/* Double save */
global S_E_glmw `glmwarn'		/* Double save */
global S_E_nobs `nobs'			/* Double save */
global S_E_sfac `scalfac'		/* Double save */
global S_E_nx=`nx'+`nbase'		/* Double save */
global S_E_opts `dead' `options'	/* Double save */
global S_E_t1t "Box-Tidwell model"	/* Double save */
global S_E_cmd `cmd'			/* Double save */
global S_E_fp fracpoly			/* Double save */
global S_E_fp2 boxtid			/* Double save */

est local fp_base `covars'
est local fp_depv `y'
est scalar fp_dist=`dist'
est local fp_fvl `fvl'
est local fp_wgt "`weight'"
est local fp_exp "`exp'"
est local fp_glmw `glmwarn'
est scalar fp_nobs=`nobs'
est local fp_sfac `scalfac'
est scalar fp_nx=`nx'+`nbase'
est local fp_opts `dead' `options'
est local fp_t1t "Box-Tidwell model"
est local fp_cmd fracpoly
est local fp_cmd2 boxtid

if `rc'==0 { _rslt }
else error `rc'
end

program define _rslt
di in gr _n "Box-Tidwell regression model"
`e(cmd)'

local nx: word count `e(fp_x)'
local v 0
while `v'<`nx' {
	local v=`v'+1
	local var: word `v' of `e(fp_x)'
	local skip=9-length("`var'")

	cap local b1=e(fp_l`v')
	if _rc { local b1 . }
	cap local seb1=e(fp_sl`v')
	if _rc { local seb1 . }
	cap local pnlin=e(fp_P`v')
	if _rc { local pnlin . }

	frac_ddp e(fp_g`v') 3
	local nlindev `r(ddp)'

	di in gr "`var'" _skip(`skip') "|" /*
	 */ in ye _col(13) %9.0g `b1' "  " %9.0g `seb1' /*
	 */ " " %10.3f `b1'/`seb1' /*
	 */ in gr _col(47) "Nonlin. dev. " in ye "`nlindev'" /*
	 */ in gr _col(68) "(P = " in ye %5.3f `pnlin' in gr ")"

	local deg: word count `e(fp_k`v')'
	local i 1
	while `i'<=`deg' {
		local k: word `i' of `e(fp_k`v')'
		local s: word `i' of `e(fp_s`v')'
		if "`k'"=="" { local k . }
		if "`s'"=="" { local s . }
		di in gr "      p`i' |  " in ye %9.0g `k' "  " %9.0g /*
		 */ `s' " " %10.3f `k'/`s'
		local i=`i'+1
	}
	di in gr _dup(78) "-"
}
di in gr "Deviance:" in ye %9.3f e(fp_dev) in gr "."
if e(fp_cnvg)==0 {
	di in red "divergence, or only partial convergence achieved"
	exit 430
}
end

program define init, rclass
	args cmd y vv lnx adj covars wgt touse deg dead zero powers options init
	local wi: word count `init'	/* initial powers can be null */
	if `wi'>0 {
		if `wi'!=`deg' {
			noi di in red /*
			 */ "incorrect number of initial values for `vv'"
			exit 198
		}
	}
	if !`lnx' 		{ local isexpx expx(sd) }
	if "`covars'`adj'"!="" 	{ local cov covars(`covars' `adj') }
	if "`init'"=="" {
		tempname vn /* creates temporary name */
		frac_154 `cmd' `y' `vv' `covars' `adj' `wgt' if `touse', /*
		 */ degree(`deg') `dead' `zero' `isexpx' powers(`powers') `options' /*
		 */ name(`vn')
		cap drop `e(fp_xp)' /* remove var(s) with funny names */
		if !`lnx' { local xpx=`e(fp_xpx)' } /* factor if expx used */
		else local xpx 1
		if `deg'==1 {
			local dev=`e(fp_d1)'
			local pwr `e(fp_k1)'
* Beginning of code taken from previous subroutine "sorted".
			local mx 1
			local mp1 `e(fp_p1)'
			local i 1
			while "`e(fp_bt`i')'"!="" {
				local p`i': word 1 of `e(fp_bt`i')'
				local d`i': word 2 of `e(fp_bt`i')'
				local s`i'=(`p`i'')^2
				if `p`i''==`mp1' { local m1 `i' }
				if `d`i''>`d`mx'' { local mx `i' }
				local i=`i'+1
			}
			local np=`i'-1
/*
	find 2nd smallest deviance
*/
			local m2 `mx'
			local i 0
			while `i'<`np'{
				local i=`i'+1
				if `d`i''<`d`m2'' & `i'!=`m1' {
					local m2 `i'
				}
			}
/*
	find 3rd smallest deviance
*/
			local m3 `mx'
			local i 0
			while `i'<`np'{
				local i=`i'+1
				if `d`i''<`d`m3''&`i'!=`m1'&`i'!=`m2' {
					local m3 `i'
				}
			}
/*
	form matrices and do quadratic regression
*/
			tempname X Y XTX XTXi XTY B
			matrix `X' = (1,`p`m1'',`s`m1''\1,`p`m2'',`s`m2''\1,`p`m3'',`s`m3'')
			matrix `Y' = (`d`m1'' \ `d`m2'' \ `d`m3'')
			matrix `XTX' = `X''*`X'
			matrix `XTXi' = syminv(`XTX')
			matrix `XTY' = `X''*`Y'
			matrix `B' = `XTXi'*`XTY'
/*
	find initial value
*/
			local pwrs = -0.5*`B'[2,1]/`B'[3,1]
* End of code taken from previous subroutine "sorted".
			if abs(`pwrs')<.05 {
				local pwrs=sign(`pwrs')*.05
			}
			frac_154 `cmd' `y' `vv' `pwrs' `covars' `adj' `wgt' /*
			 */ if `touse', /*
			 */ `dead' `zero' `isexpx' `options' name(`vn')
			cap drop `e(fp_xp)'
			local devsrt=e(fp_dlin)
			local pwrsrt `e(fp_k1)'
			if `devsrt'>`dev' {
				local pwrs `pwr'
			}
			else local pwrs `pwrsrt'
		}
		else local pwrs `e(fp_pwrs)'
		tokenize `pwrs'
		local pwrs
		while "`1'"!="" {
			if abs(`1')<.05 { local 1 .05 }
			local pwrs `pwrs' `1'
			mac shift
		}
	}
	else {
/*
	Find initial power values in string `init'
*/
		local i 0
		local pwrs
		while `i'<`deg' {
			local i=`i'+1
			local p: word `i' of `init'
			cap confirm num `p'
			if _rc {
				noi di in red "invalid init() for var `vv'"
				exit 198
			}
			local pwrs `pwrs' `p'
		}
		local xpx 1
	}
	if `lnx' { local xpx 1 }
	return local pwrs `pwrs'	/* initialised powers */
	return scalar f = `xpx'		/* factor from expx transf */
end

program define initial
* initialize powers
	args deg f pwrs
	local i 1
	while `i'<=`deg' {
		local p: word `i' of `pwrs'
		if abs(`p')<.001 { local p=sign(`p')*.001 }
		if `i'==1 {
			local p1 `p'
			local padd 0
		}
		else {
			if abs(`p'-`p1')<.001 { /* tiebreak */
				local padd=`padd'+0.22
				local p=`p'+`padd'
			}
			else {
				local p1 `p'
				local padd 0
			}
		}
		if `f'!=1 { global S_`i'=`p'*`f' }
		else global S_`i' `p'
		local i=`i'+1
	}
end
*! version 1.0.0 PR 26Feb1999.
program define frac_str, sclass
version 6
gettoken target 0 : 0
if "`target'"=="" { error 198 }
tokenize `0'
local found 0
while "`1'"!="" {
	if "`1'"=="`target'" { local found = `found'+1 }
	else local new "`new' `1'"
	mac shift
}
sret local new `new'
sret local found `found'
end
*! version 1.0.1 PR 01Mar1999.
program define frac_adj, rclass
	version 6
	* Inputs: 1=macro `adjust', 2=varlist, 3=case filter.
	* Returns adjustment values in r(adj1),...
	* Returns number of unique values in r(uniq1),...

	args adjust rhs touse
	if "`adjust'"=="" { frac_dis mean adjust "`rhs'" }
	else		    frac_dis "`adjust'" adjust "`rhs'"
	tokenize `rhs'
	tempname u
	local i 1
	while "``i''"!="" {
		quietly inspect ``i'' if `touse'
		scalar `u'=r(N_unique)
		local a ${S_`i'}
		if "`a'"=="" | "`adjust'"=="" {	/* identifies default cases */
			if `u'==1 {		/* no adjustment */
				local a
			}
			else if `u'==2 {	/* adjust to min value */
				quietly summarize ``i'' if `touse', meanonly
				if r(min)==0 { local a }
				else local a=r(min)
			}
			else local a mean
		}
		else if "`a'"=="no" {
			local a
		}
		else if "`a'"!="mean" {
			confirm num `a'
		}
		return local adj`i' `a'
		return scalar uniq`i'=`u'
		local i=`i'+1
	}
end