*! version 1.0.0 18apr2012
program define xfrac_154, eclass
local vv : di "version " string(_caller()) ":"
version 11.0 /* 8 */
if replay() {
if "`e(cmd)'" == "" | "`e(fp_cmd2)'" != "fracpoly" {
error 301
}
syntax [, COMpare SEQuential *]
if `"`e(cmd2)'"' != "" {
di in blue `"->`e(cmd2)'"'
`e(cmd2)', `options'
}
else {
di in blue `"->`e(cmd)'"'
`e(cmd)', `options'
}
fraccomp, `compare' `sequential'
exit
}
gettoken cmd 0 : 0, parse(" ,")
xfrac_chk `cmd'
if `s(bad)' {
di as err "invalid or unrecognised command, `cmd'"
exit 198
}
/*
dist=0 (normal), 1 (binomial), 2 (poisson), 3 (cox), 4 (glm),
5 (xtgee), 6 (ereg/weibull), 7 (stcox/streg).
*/
local dist `s(dist)'
local glm `s(isglm)'
local qreg `s(isqreg)'
local xtgee `s(isxtgee)'
local normal `s(isnorm)'
if `dist' != 7 {
gettoken lhs 0 : 0, parse(" ,")
if `dist' == 8 {
gettoken lhs2 0 : 0, parse(" ,")
local lhs `lhs' `lhs2'
}
}
gettoken star 0 : 0, parse(" ,")
local star "`lhs' `star'"
/*
Look for fixpowers
*/
local done 0
gettoken nxt : 0, parse(" ,")
while !`done' {
local done 1
if "`nxt'"!="" & "`nxt'"!="," {
cap confirm num `nxt'
if _rc==0 {
local fixpowe `fixpowe' `nxt'
local done 0
gettoken nxt 0 : 0, parse(" ,")
gettoken nxt : 0, parse(" ,")
}
}
}
local 0 `"`star' `0'"'
local search = ("`fixpowe'"=="")
if `search' {
local srchopt ADDpowers(str) DEVthr(real 1e30) /*
*/ POwers(numlist) DEGree(int 2) LOg COMpare SEQuential
}
if `dist' == 8 local minv 3
else if `dist' != 7 local minv 2
else local minv 1
syntax varlist(min=`minv' fv) [if] [in] [aw fw pw iw] [, ///
restrict(string) DEAD(str) CATzero EXPx(str) ORIgin(str) ZERo ///
noCONstant noSCAling ADJust(str) CENTer(str) NAme(str) `srchopt' * ]
local small 1e-6
if ("`adjust'"=="") {
local adjust "`center'"
}
else if ("`center'"!="") {
di as err "may not specify both adjust() and center()"
exit 198
}
if "`constant'"=="" & "`cmd'"=="regress" & "`options'"!="" {
/* -regress- has syntax noConstant */
/* do not let it quietly passthru to -regress- */
/* in the options macro */
CheckForNoConstant, `options'
local constant `s(constant)'
local options `"`s(options)'"'
}
frac_cox "`dead'" `dist'
local cz="`catzero'"!=""
if `cz' local zero zero
if `search' {
if `degree'<1 | `degree'>9 {
di as err "invalid degree()"
exit 198
}
local df=2*`degree'
local odddf=(2*int(`df'/2)<`df')
}
else local df 1
local lin=("`fixpowe'"=="1") & ("`zero'"=="") & ("`catzero'"=="")
if "`constant'"=="noconstant" {
if "`cmd'"=="fit" | "`cmd'"=="cox" | "`cmd'"=="stcox" | /*
*/ "`cmd'"=="streg" | "`cmd'"=="clogit" | /*
*/ "`cmd'"=="logistic" {
di as err "noconstant invalid with `cmd'"
exit 198
}
local options "`options' nocons"
}
/*
Read powers to be searched into a variable, `p'
*/
if `search' {
if "`powers'"=="" local powers "-2,-1,-.5,0,.5,1,2,3"
local pwrs "`powers' `addpowers'"
frac_pq "`pwrs'" 1 1
local np `r(np)'
if `np'<1 {
di as err "no powers given"
exit 2002
}
local pwrs "`r(powers)'"
forvalues i=1/`np' {
local p`i' `r(p`i')'
}
}
ereturn clear /* do we want to do this ?? */
tokenize `varlist'
if `dist' == 8 { // intreg
local y1 `1'
local y2 `2'
local y `y1' `y2'
local rhs `3'
mac shift 3
}
else if `dist' != 7 {
local y `1'
local rhs `2'
mac shift 2
}
else { // dist = 7: stcox, streg, stpm
local y _t
local rhs `1'
mac shift
}
fvexpand `rhs'
if "`r(fvops)'" == "true" {
di as err "factor variables not allowed in this context"
exit 198
}
local base `*'
tempvar touse x lnx
quietly {
mark `touse' [`weight' `exp'] `if' `in'
if `dist' == 8 {
replace `touse' = 0 if `y1'>=. & `y2'>=.
markout `touse' `rhs' `base' `dead'
}
else markout `touse' `rhs' `y' `base' `dead'
lab var `touse' "fracpoly sample"
if "`dead'"!="" {
local options "`options' dead(`dead')"
}
if "`restrict'"!="" {
/*
-restrict()- becomes part of touse for estimation purposes;
only for fracgen is the original touse used.
*/
tempvar touse_user
gen byte `touse_user' = 1 `if' // ! PR bug fix
frac_restrict `touse' `restrict'
local restrict restrict(`restrict')
}
else local touse_user `touse'
/*
Deal with weights.
*/
frac_wgt `"`exp'"' `touse' `"`weight'"'
local mnlnwt = r(mnlnwt) /* mean log normalized weights */
local wgt `r(wgt)'
gen double `x'=`rhs' if `touse'
gen double `lnx'=.
frac_xo `x' `lnx' `lin' "`expx'" "`origin'" /*
*/ "`zero'" "`scaling'" `rhs' `touse'
local nobs = r(N)
if r(shifted)==0 {
local zeta `r(zeta)'
local shift 0
}
else local shift `r(zeta)'
local kx `"`r(expxest)'"'
local scalfac `r(scale)'
if `cz' {
tempvar catz
gen byte `catz'=`x'<=0
}
}
if `dist' != 7 local yvar `y'
/*
Determine residual and model df.
*/
qui regress `y' `base' `wgt' if `touse'
local rdf=e(df_r)+("`constant'"=="noconstant")
/*
Calc deviance=-2(log likelihood) for regression on covars only,
allowing for possible weights.
*/
if "`constant'"!="" & `dist'>0 {
if `dist' == 1 {
/* logit will allow zero vector to be dropped */
/* by _rmcoll and complete w/out error */
tempvar z0
qui gen `z0' = 0
}
else {
/* all others dev0 = 0 */
local scale = 1
local dev0 = .
}
}
if "`dev0'" == "" {
`vv' ///
qui `cmd' `yvar' `z0' `base' `wgt' if `touse', `options'
if `xtgee' & "`base'"=="" global S_E_chi2 0
if `glm' {
/* Note: with Stata 8 scale param is e(phi); */
/* was e(delta) in Stata 6 */
/* Also e(dispersp) has become e(dispers_p). */
local scale 1
if abs(e(dispers_p)/e(phi)-1)>`small' & /*
*/ abs(e(dispers)/e(phi)-1)>`small' {
local scale = e(phi)
}
}
frac_dv `normal' "`wgt'" `nobs' `mnlnwt' `dist' /*
*/ `glm' `xtgee' `qreg' "`scale'"
local dev0 = r(deviance)
if `normal' local rsd0= e(rmse)
}
if `search' {
local m `degree'
if "`log'"!="" {
di _n as txt "Model # Deviance " _cont
if `normal' di as txt " Res S.D. " _cont
forvalues i=1/`m' {
di as txt " Power " `i' _cont
}
di _n
}
/*
Go!
*/
forvalues i=1/`np' {
local pi `p`i''
tempvar x`i'
qui gen double `x`i''=cond(`pi'==0, `lnx', /*
*/ cond(`pi'==1, `x', cond(`x'==0,0,`x'^`pi')))
}
local i 0
while `i'<`m' {
local j`i' 0
local ++i
tempvar xp`i'
qui gen double `xp`i''=.
local j=`i'+`i'-1
local dev`j' `dev0' /* min deviance for df=j */
local ++j
local dev`j' `dev0'
}
local j`m' 1
local i `m'
local kount 0 /* no. of models processed */
local kountd 0 /* no. of models with deviance<`devthr' */
local deg `j`m''
local done 0
while !`done' {
if `i'<`m' {
local ji `j`i''
qui replace `xp`i''=`x`ji''
local l `i'
while `l'<`m' {
local l=`l'+1
local j`l' `ji'
local 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)
*/
local one 0
local l 0
while `l'<`m' {
local ++l
if "`p`j`l'''"=="1" local one 1
}
if !`odddf' | (`deg'<`m') | `one' {
local dfi=2*`deg'-`one'
local i 1
local xlist
while `i'<=`m' {
if `j`i''>0 {
local xlist "`xlist' `xp`i''"
}
local ++i
}
local ++kount
`vv' ///
cap `cmd' `yvar' `xlist' `catz' `base' `wgt' /*
*/ if `touse', `options'
if _rc == 0 {
frac_dv `normal' "`wgt'" `nobs' /*
*/ `mnlnwt' `dist' `glm' `xtgee' /*
*/ `qreg' "`scale'"
local dev = r(deviance)
if `normal' {
local rsd=e(rmse)
if `dfi'==1 local rsd1 `rsd'
}
}
else if _rc==1 {
error 1
}
else local dev = .
if `dev'<`dev`dfi'' {
local dev`dfi' `dev'
if `normal' local rsd`dfi' `rsd'
local i 1
local pm`dfi'
while `i'<=`m' {
if `j`i''>0 {
local pm `p`j`i'''
local pm`dfi' /*
*/" `pm`dfi'' `pm'"
}
local ++i
}
}
if `dev'<`devthr' local ++kountd
if "`log'"!="" {
di as res %5.0f `kount' /*
*/ _col(10) %12.3f `dev' _cont
if `normal' di " " %8.0g `rsd' _cont
local i `m'
while `i'>0 {
if `j`i''==0 di _skip(6) ". " _cont
else di as res %8.1f `p`j`i''' _cont
local --i
}
di
}
if `deg'==1 {
local fppow `p`j`m'''
local fpdev`j`m'' "`fppow' `dev'"
}
}
/*
Increment the first possible index (of loop i) among indices of
loops m, m-1, m-2, ..., 1
*/
local i `m'
/*
Finish after all indices have achieved their upper limits (np).
*/
while `j`i''==`np' {
local --i
}
if `i'==0 local done 1
else {
if `j`i''==0 local deg = `m'-`i'+1
local ++j`i'
}
}
/*
Update the results for even df to include odd df
*/
local i 2
while `i'<=`df' {
local j=`i'-1
if `dev`j''<`dev`i'' {
local dev`i' `dev`j''
if `normal' local rsd`i' `rsd`j''
local pm`i' "`pm`j''"
}
local i=`i'+2
}
if "`log'"=="" di
}
/*
Create FP transformation(s) of xvar for final model
*/
if `search' local pwrs `pm`df''
else local pwrs `fixpowe'
if "`expx'"!="" local e "expx(`expx')"
if "`origin'"!="" local o "origin(`origin')"
if "`adjust'"!="" local a "adjust(`adjust')"
if "`name'"!="" local n "name(`name')"
if !`lin' | (`lin' & "`e'`o'`a'"!="") {
fracgen `rhs' `pwrs' if `touse_user', `restrict' sayesamp replace /*
*/ `zero' `catzero' `scaling' `a' `e' `o' `n'
local xp `r(names)'
}
else local xp `rhs'
/*
Fit final model with permanent e(sample)=`touse' filter
*/
`vv' ///
`cmd' `yvar' `xp' `base' `wgt' if `touse', `options'
if !`search' {
/*
Deviance for fixed-powers model.
Note that `dev1' is stored in e(fp_dlin), deviance for a
linear model, even when `fixpowe' is not 1; similarly
for `rsd1'.
*/
frac_dv `normal' "`wgt'" `nobs' `mnlnwt' `dist' `glm' /*
*/ `xtgee' `qreg' "`scale'"
local dev1 = r(deviance)
if `normal' local rsd1=e(rmse)
local pm1 `fixpowe'
}
di as txt "Deviance:" as res %9.2f `dev`df'' as txt ". " _cont
local dev `dev`df''
if `search' {
di as txt "Best powers of " as res "`rhs'" as txt " among " /*
*/ as res `kount' as txt " models fit: " as res "`pwrs'" _cont
if `devthr'<1e30 {
di _n as txt "Number of models with deviance < " /*
*/ as res `devthr' as txt ": " as res `kountd' _cont
}
di as txt "." _cont
}
di
ereturn scalar fp_d0 = `dev0'
cap ereturn scalar fp_s0 = `rsd0'
ereturn scalar fp_dlin = `dev1'
global S_E_d0 `dev0' /* double save */
global S_E_s0 `rsd0' /* double save */
global S_E_dlin `dev1' /* double save */
if `normal' {
ereturn scalar fp_slin = `rsd1'
* global S_E_slin `rsd1' /* double save */
}
local i 2
local j 1
while `i'<=`df' {
ereturn scalar fp_d`j' = `dev`i''
ereturn local fp_p`j' `pm`i'' /* PR bug fix */
global S_E_d`j' `dev`i'' /* double save */
global S_E_p`j' `pm`i'' /* double save */
if `normal' {
ereturn scalar fp_s`j' = `rsd`i''
global S_E_s`j' `rsd`i'' /* double save */
}
local i=`i'+2
local j=`j'+1
}
/*
New code in v 1.4.7 for consistency with mfracpol
*/
ereturn local fp_x1 `rhs'
ereturn local fp_k1 `pwrs'
local nbase: word count `base'
local i 0
while `i'<`nbase' {
local i=`i'+1
local j=`i'+1
ereturn local fp_x`j' : word `i' of `base'
ereturn local fp_k`j' 1
}
/*
End of new code in v 1.4.7 for consistency with mfracpol
*/
if `search' {
local i `np'
while `i'>0 {
ereturn local fp_bt`i' `fpdev`i''
local i=`i'-1
}
}
ereturn scalar fp_dev = `dev'
ereturn scalar fp_catz = `cz'
ereturn scalar fp_nx = `nbase'+1
ereturn scalar fp_df = `df'
ereturn scalar fp_rdf = `rdf'
ereturn scalar fp_N = `nobs'
ereturn local fp_opts `options'
ereturn local fp_t1t "Fractional Polynomial"
ereturn local fp_pwrs `pwrs'
ereturn local fp_wgt "`weight'"
ereturn local fp_wexp "`exp'"
ereturn local fp_xp `xp'
ereturn local fp_base `base'
ereturn local fp_rhs `rhs'
ereturn local fp_depv `y'
ereturn local fp_fvl `xp' `base'
ereturn scalar fp_dist = `dist'
ereturn local fp_fprp "no"
ereturn scalar fp_srch = `search'
ereturn scalar fp_sfac = `scalfac'
ereturn scalar fp_shft = `shift'
ereturn local fp_xpx `kx'
ereturn local fp_cmd "fracpoly"
ereturn local fp_cmd2 "fracpoly"
global S_E_base `base' /* double save */
global S_E_df `df' /* double save */
global S_E_depv `y' /* double save */
global S_E_wgt "`weight'" /* double save */
global S_E_exp "`exp'" /* double save */
global S_E_fp fracpoly /* double save */
global S_E_fp2 fracpoly /* double save */
global S_E_nobs `nobs' /* double save */
global S_E_pwrs `pwrs' /* double save */
global S_E_rdf `rdf' /* double save */
global S_E_rhs `rhs' /* double save */
global S_E_xp `xp' /* double save */
global S_E_cmd `cmd' /* double save */
fraccomp, `compare' `sequential'
end
program define fraccomp /* report model comparisons */
syntax [, COMpare SEQuential ]
if "`compare'`sequential'"=="" | e(fp_srch)==0 exit
local normal=(e(fp_dist)==0)
local cz = e(fp_catz)
if `cz' local catz " (+)"
local dash " {hline 2}"
local ddup=63+15*`normal'
di as txt _n "Fractional polynomial model comparisons:"
di as txt "{hline `ddup'}"
di as txt abbrev("`e(fp_rhs)'",12) _col(18) "df Deviance " _cont
if `normal' di as txt " Res. SD " _cont
di as txt "Dev. dif. P (*) Powers"
di as txt "{hline `ddup'}"
local maxdf = e(fp_df)
local degree = 1+int((`maxdf'-1)/2)
// PR 22mar2009 -start-
frac_eqmodel k
if "`sequential'"=="" { // default is closed-test-style presentation
local dev0 = e(fp_d`degree')
local j 0
while `j'<=`maxdf' {
local m = 1 + int((`j' - 1) / 2) // degree of model being reported
if `j' == 0 local idf 0
else local idf = cond(`j' == 1, `k' * (1 + `cz'), (`k' + 1) * `m' + `k' * `cz')
if `j'<`maxdf' {
// denominator df for F-tests (`n1' is numerator)
local n2 = e(fp_rdf) - `idf'
if `j'==0 {
di as txt "Not in model" _col(18) as res %2.0f `idf' _col(22) _cont
local dev = e(fp_d0)
local rsd = e(fp_s0)
if `maxdf' == 1 local n1 = `k' * (1 + `cz')
else local n1 = (`k' + 1) * `degree' + `k' * `cz'
local pwr
}
else if `j'==1 {
di as txt "Linear`catz'" _col(18) as res %2.0f `idf' _col(22) _cont
local dev = e(fp_dlin)
local rsd = e(fp_slin)
local n1 = (`k' + 1) * `degree' - `k'
local pwr 1
}
else {
di as txt "m = `m'`catz'" _col(18) as res %2.0f `idf' _col(22) _cont
local dev = e(fp_d`m')
local rsd = e(fp_s`m')
local n1 = (`degree' - `m') * (`k' + 1)
local pwr `e(fp_p`m')'
}
local d = `dev'-`dev0'
frac_pv `normal' "`e(fp_wgt)'" `e(fp_N)' `d' `n1' `n2'
local P = r(P)
if `j'==1 global S_E_Plin `P'
else global S_E_P`m' `P'
di as res %13.3f `dev' _cont
if `normal' di as res _skip(5) %8.0g `rsd' _cont
di as res _skip(3) %7.3f `d' %9.3f `P' _skip(2) "`pwr'"
}
else {
di as txt "m = `degree'`catz'" _col(18) as res %2.0f `idf' /*
*/ _col(22) %13.3f `dev0' _cont
if `normal' di as res _skip(5) %8.0g e(fp_s`degree') _cont
di as txt _skip(1) "`dash'`dash'" as res _skip(2) "`e(fp_p`degree')'"
}
local j = `j' + cond(`j'<2, 1, 2)
}
}
else {
di as txt "Not in model" as res _col(18) " 0" /*
*/ _col(22) %13.3f e(fp_d0) _cont
if `normal' di as res _skip(5) %8.0g e(fp_s0) _cont
di as txt _skip(1) "`dash'`dash'"
local i 1
while `i'<=`maxdf' {
local m = 1 + int((`i' - 1) / 2) // degree of model being reported
local idf = cond(`i' == 1, `k' * (1 + `cz'), (`k' + 1) * `m' + `k' * `cz')
if `i'==1 {
di as txt "Linear`catz'" _col(18) as res /*
*/ %2.0f `idf' _col(22) _cont
local dev = e(fp_dlin)
local rsd = e(fp_slin)
local d = e(fp_d0)-`dev'
local n1 `idf'
local n2 = e(fp_rdf) - `idf'
local pwr 1
}
else {
di as txt "m = `m'`catz'" _col(18) as res /*
*/ %2.0f `idf' _col(22) _cont
local dev = e(fp_d`m')
local rsd = e(fp_s`m')
local d = `devlast'-`dev'
local n1 = cond(`m'==1, 1, `k' + 1)
local n2 = e(fp_rdf) - `idf'
local pwr `e(fp_p`m')'
}
// PR 22mar2009 -end-
frac_pv `normal' "`e(fp_wgt)'" `e(fp_N)' `d' `n1' `n2'
local P = r(P)
if `i'==1 global S_E_Plin `P'
else global S_E_P`m' `P'
di as res %13.3f `dev' _cont
if `normal' di as res _skip(5) %8.0g `rsd' _cont
di as res _skip(3) %7.3f `d' %9.3f `P' _skip(2) "`pwr'"
local devlast `dev'
local i = `i' + cond(`i'==1, 1, 2)
}
}
di as txt "{hline `ddup'}"
if "`sequential'"=="" ///
di as txt "(*) P-value from deviance difference comparing reported model with m = `degree' model"
else ///
di as txt "(*) P-value from dev. dif. comparing reported model with next lower df model"
if `cz' ///
di as txt "(+) indicates dummy variable for `e(fp_rhs)'==0 included in model"
end
program CheckForNoConstant, sclass
syntax [, noConstant * ]
sreturn local constant `constant'
sreturn local options `"`options'"'
end
exit
1.2.0/PR/JSP:- replace -all- option with -restrict()- option
1.1.13/PR:- mods in output for compare option for Stata 10.
1.1.12/PR:- translated to Stata 8
1.1.11/JML:- support for extended missing values
1.1.10/PR:- fixed bug in support for streg, stcox
1.1.9/RG:- added support for streg, stcox
1.1.8/WG:- %9.2g changed to %9.2f
1.1.4/WG:- minor title changes
1.1.3/PR:- saves adjustments stored by -fracgen-.
- change frac_pp to frac_pq and avoid storing powers in a variable.
1.1.2/PR:- kludge to avoid bug in -glm- with no covariate and noconstant.
1.1.1/PR:- bug fixed, not saving fp_p`j'. Flag fp_cmd2 added and checked.
1.1.0/WG:- translated to Stata 6
1.0.2/PR:- `all' option added to facilitate out-of-sample predictions.
Double precision for transformations of xvar.
1.0.1/PR:- added name() option to be passed to fracgen to enable more
careful name control by calling program (fracpoly or mfracpol).
Version 1.5.4 of old fracpoly (12Jul98).
Now used as slave routine for fracpoly.ado and mfracpol.ado.