*! version 1.0.2 PR 09ap2015 program define frac_154_wald, eclass version 9.0 /* Now uses mim rather than micombine, PR 09apr2015 Support for -svy- added 14sep2006. Modified for wald testing, PR 28mar2006. Version 1.5.4 of old fracpoly (12Jul98). Now used as slave routine for fracpoly.ado and mfracpol.ado. */ if replay() { if "`e(cmd)'" == "" | "`e(fp_cmd2)'" != "fracpoly" { error 301 } syntax [, COMpare *] if `"`e(cmd2)'"' != "" { di in blue `"->`e(cmd2)'"' `e(cmd2)', `options' } else { di in blue `"->`e(cmd)'"' `e(cmd)', `options' } fraccomp, `compare' exit } // Note that syntax is "mim ...". mim invokes mim, storebv: gettoken micmb 0 : 0 if "`micmb'"!="mim" { local cmd `micmb' local micmb } else { local micmb mim, storebv: gettoken cmd 0 : 0, parse(" ,") } 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), 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(" ,") } 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) POwers(numlist) DEGree(int 2) LOg COMpare } if `dist' != 7 { local minv 2 } else local minv 1 syntax varlist(min=`minv') [if] [in] [aw fw pw iw] [, /* */ ALL DEAD(str) CATzero EXPx(str) ORIgin(str) ZERo /* */ noCONStant noSCAling ADJust(str) NAme(str) `srchopt' /* */ SVYalone svy(string) * ] local small 1e-6 frac_cox "`dead'" `dist' local cz="`catzero'"!="" if `cz' local zero zero if `search' { if `degree'<1 | `degree'>9 { di in red "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" { di in red "noconstant invalid with `cmd'" exit 198 } local options "`options' nocons" } if "`svyalone'"!="" local svy svy: else if `"`svy'"'!="" local svy svy, `svy': /* 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' `addpowe'" frac_pq "`pwrs'" 1 1 local np `r(np)' if `np'<1 { di in red "no powers given" exit 2002 } local pwrs "`r(powers)'" local i 1 while `i'<=`np' { local p`i' `r(p`i')' local i=`i'+1 } } ereturn clear /* do we want to do this ?? */ tokenize `varlist' if `dist' != 7 { local y `1' local rhs `2' mac shift 2 } else { local y _t local rhs `1' mac shift } local base `*' tempvar touse x lnx quietly { mark `touse' [`weight' `exp'] `if' `in' markout `touse' `rhs' `y' `base' `dead' lab var `touse' "fracpoly sample" if "`dead'"!="" { local options "`options' dead(`dead')" } /* 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 } } /* Store coefficients for linear fit on untransformed x */ if `dist' != 7 { local yvar `y' } qui `svy'`cmd' `yvar' `rhs' `base' `wgt' if `touse', `options' cap local b0=_b[_cons] if _rc local b0 0 cap local b1=_b[`rhs'] if _rc { di in red "could not fit linear model for `rhs'" exit 2001 } /* Determine residual and model df. */ qui reg `y' `base' `wgt' if `touse' local rdf=e(df_r)+("`constant'"=="noconstant") if `search' { local m `degree' if "`log'"!="" { di _n in gr "Model # Wald stat" _cont if `normal' di in gr " Res S.D. " _cont local i 1 while `i'<=`m' { di in gr " Power " `i' _cont local i=`i'+1 } di _n } /* Go! */ local i 1 while `i'<=`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=`i'+1 } local i 0 while `i'<`m' { local j`i' 0 local i=`i'+1 tempvar xp`i' qui gen double `xp`i''=. local j=`i'+`i'-1 local wald`j' 0 /* min Wald stat for df=j */ local j=`j'+1 local wald`j' 0 } local j`m' 1 local i `m' local kount 0 /* no. of models processed */ 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=`l'+1 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=`i'+1 } local kount=`kount'+1 qui `micmb' `svy'`cmd' `yvar' `xlist' `catz' `base' `wgt' /* */ if `touse', `options' capture frac_wald `xlist' if c(rc) { // can happen if a variable is omitted from `xlist' due to collinearity if `normal' local rsd . local wald . } else { local wald = r(wald) if `normal' local rsd=e(rmse) } if (`wald' < .) & (`wald' > `wald`dfi'') { local wald`dfi' `wald' 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=`i'+1 } } if "`log'"!="" { di in ye %5.0f `kount' /* */ _col(10) %12.3f `wald' _cont if `normal' { di " " %8.0g `rsd' _cont } local i `m' while `i'>0 { if `j`i''==0 { di _skip(6) ". " _cont } else /* */ di in ye %8.1f `p`j`i''' _cont local i=`i'-1 } di } if `deg'==1 { local fppow `p`j`m''' local fpwald`j`m'' "`fppow' `wald'" } } /* 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=`i'-1 } if `i'==0 local done 1 else { if `j`i''==0 local deg = `m'-`i'+1 local j`i'=`j`i''+1 } } /* Update the results for even df to include odd df */ local i 2 while `i'<=`df' { local j=`i'-1 if `wald`j''>`wald`i'' { local wald`i' `wald`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', `all' sayesamp replace /* */ `zero' `catzero' `scaling' `a' `e' `o' `n' local xp `r(names)' if "`adjust'"!="" { local j 1 local Np: word count `pwrs' while `j'<=`Np' { local a`j'=r(adj`j') local j=`j'+1 } } } else local xp `rhs' /* Fit final model with permanent e(sample)=`touse' filter */ `micmb' `svy'`cmd' `yvar' `xp' `base' `wgt' if `touse', `options' if !`search' { /* Fixed-powers model. Note that `wald1' is stored in e(fp_dlin), deviance for a linear model, even when `fixpowe' is not 1; similarly for `rsd1'. */ frac_wald `xp' local wald1 = r(wald) if `normal' local rsd1=e(rmse) local pm1 `fixpowe' } di in gr "Wald stat:" in ye %9.2f `wald`df'' in gr ". " _cont local wald `wald`df'' 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 di in gr "." _cont } di ereturn scalar fp_d0 = 0 cap ereturn scalar fp_s0 = `rsd0' ereturn scalar fp_dlin = `wald1' if `normal' { ereturn scalar fp_slin = `rsd1' } local i 2 local j 1 while `i'<=`df' { ereturn scalar fp_d`j' = `wald`i'' ereturn local fp_p`j' `pm`i'' /* PR bug fix */ if `normal' { ereturn scalar fp_s`j' = `rsd`i'' } 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 } } if "`adjust'"!="" { local j 1 while `j'<=`Np' { ereturn scalar fp_a`j'=`a`j'' local j=`j'+1 } } ereturn scalar fp_wald = `wald' 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 scalar f_b0 = `b0' ereturn scalar f_b1 = `b1' 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" fraccomp, `compare' end program define fraccomp /* report model comparisons */ syntax [, COMpare] if "`compare'"=="" | e(fp_srch)==0 exit local normal=(e(fp_dist)==0) local cz = e(fp_catz) if `cz' local catz " + 0" local dash " {hline 2}" local ddup=63+15*`normal' di in gr _n "Fractional polynomial model comparisons:" di in smcl in gr "{hline `ddup'}" di in gr abbrev("`e(fp_rhs)'",12) _col(18) "df Wald stat " _cont if `normal' di in gr " Res. SD " _cont di in gr " Gain P(term) Powers" di in smcl in gr "{hline `ddup'}" di in gr "Not in model" in ye _col(18) " 0" /* */ _col(22) %13.3f e(fp_d0) _cont if `normal' di in ye _skip(5) %8.0g e(fp_s0) _cont di in smcl in gr _skip(3) "`dash'`dash'" local i 1 while `i'<=e(fp_df) { local m=1+int((`i'-1)/2) local idf=`i'+`cz' if `i'==1 { di in gr "Linear`catz'" _col(18) in ye /* */ %2.0f `idf' _col(22) _cont local wald = e(fp_dlin) local rsd = e(fp_slin) local d=-(e(fp_d0)-`wald') local n1=1+`cz' local n2=e(fp_rdf)-1-`cz' local pwr 1 } else { di in gr "m = `m'`catz'" _col(18) in ye /* */ %2.0f `idf' _col(22) _cont local wald = e(fp_d`m') local rsd = e(fp_s`m') local d=-(`waldlast'-`wald') local n1 = cond(`m'==1, 1, 2) local n2=e(fp_rdf)-2*`m'-`cz' *local pwr ${S_E_p`m'} local pwr `e(fp_p`m')' } 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 in ye %13.3f `wald' _cont if `normal' { di in ye _skip(5) %8.0g `rsd' _cont } di in ye _skip(3) %7.3f -(e(fp_dlin)-`wald') /* */ %9.3f `P' _skip(2) "`pwr'" local waldlast `wald' local i = `i' + cond(`i'==1,1,2) } di in smcl in gr "{hline `ddup'}" if `cz' { di in bl /* */ "[Note:`catz' indicates dummy variable for `e(fp_rhs)'=0 included in model]" } end