frac_154.ado

是一个经济学管理应用软件 很难找的 但是经济学学生又必须用到

*! version 1.1.11 PR  17sep2004
/*
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).
*/
program define frac_154, eclass
	version 6, missing
/*
	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
	}

	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).
	*/
	loc dist `s(dist)'
	loc glm `s(isglm)'
	loc qreg `s(isqreg)'
	loc xtgee `s(isxtgee)'
	loc normal `s(isnorm)'

	if `dist' != 7 {
		gettoken lhs 0 : 0, parse(" ,")
	}
	gettoken star 0 : 0, parse(" ,")
	local star "`lhs' `star'"
	/*
		Look for fixpowers
	*/
	loc 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'"'
	loc search = ("`fixpowe'"=="")
	if `search' {
		local srchopt ADDpowers(str) DEVthr(real 1e30) /*
			*/ 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' * ]
	loc small 1e-6

	frac_cox "`dead'" `dist'
	loc cz="`catzero'"!=""
	if `cz' { loc zero zero }
	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") & ("zero'"=="") & ("catzero'"=="")

	if "`constan'"=="noconstant" {
		if "`cmd'"=="fit" | "`cmd'"=="cox" | "`cmd'"=="stcox" | /*
			*/ "`cmd'"=="streg" {
			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'"
		frac_pq "`pwrs'" 1 1
		loc np `r(np)'
		if `np'<1 {
			di in red "no powers given"
			exit 2002
		}
		loc pwrs "`r(powers)'"
		local i 1
		while `i'<=`np' {
			local p`i' `r(p`i')'
			local i=`i'+1
		}
	}
	est 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
	}
	loc 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'"'
		loc mnlnwt = r(mnlnwt) /* mean log normalized weights */
		loc wgt `r(wgt)'
		gen double `x'=`rhs' if `touse'
		gen double `lnx'=.
		frac_xo `x' `lnx' `lin' "`expx'" "`origin'" /*
			*/ "`zero'" "`scaling'" `rhs' `touse'
		loc nobs = r(N)
		if r(shifted)==0 {
			loc zeta `r(zeta)'
			loc shift 0
		}
		else 	loc shift `r(zeta)'
		loc kx `"`r(expxest)'"'
		loc 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 `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'
	loc 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
		qui gen `z0'=0
	}
	qui `cmd' `yvar' `z0' `base' `wgt' if `touse', `options'
	if `xtgee' & "`base'"=="" { global S_E_chi2 0 }
	loc glmwarn 0
	if `glm' {
 		loc scale 1
 		if abs(e(dispersp)/e(delta)-1)>`small' & /*
		*/ abs(e(dispers)/e(delta)-1)>`small' {
			loc scale = e(delta)
 		}
 		else 	loc glmwarn 1
	}
	frac_dv `normal' "`wgt'" `nobs' `mnlnwt' `dist' /*
		*/ `glm' `xtgee' `qreg' "`scale'"
	loc dev0 = r(deviance)
	if `normal' { loc rsd0= e(rmse) }
	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 double `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 double `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' `yvar' `xlist' `catz' `base' `wgt' /*
					*/ if `touse', `options'
				frac_dv `normal' "`wgt'" `nobs' `mnlnwt' /*
					*/ `dist' `glm' `xtgee' `qreg' "`scale'"
				loc dev = r(deviance)
				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''' _c
						loc i=`i'-1
					}
					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
	*/
			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 }
	}
	/*
		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 "`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'
		loc 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 loc xp `rhs'
	/*
		Fit final model with permanent e(sample)=`touse' filter
	*/
	`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'"
		loc dev1 = r(deviance)
		if `normal' { loc rsd1=_result(9) }
		loc pm1 `fixpowe'
	}
	di in gr "Deviance:" in ye %9.2f `dev`df'' in gr ". " _cont
	loc dev `dev`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
		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
	est scalar fp_d0 = `dev0'
	cap est scalar fp_s0 = `rsd0'
	est 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' {
		est scalar fp_slin = `rsd1'
		global S_E_slin `rsd1' 			/* double save */
	}
	loc i 2
	loc j 1
	while `i'<=`df' {
		est scalar fp_d`j' = `dev`i''
		est 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' {
			est scalar fp_s`j' = `rsd`i''
			global S_E_s`j' `rsd`i''	/* double save */
		}
		loc i=`i'+2
		loc j=`j'+1
	}
	/*
		New code in v 1.4.7 for consistency with mfracpol
	*/
	est local fp_x1 `rhs'
	est local fp_k1 `pwrs'
	local nbase: word count `base'
	local i 0
	while `i'<`nbase' {
		local i=`i'+1
		local j=`i'+1
		est local fp_x`j' : word `i' of `base'
		est 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 {
			est local fp_bt`i' `fpdev`i''
			local i=`i'-1
		}
	}
	if "`adjust'"!="" {
		local j 1
		while `j'<=`Np' {
			est scalar fp_a`j'=`a`j''
			local j=`j'+1
		}
	}
	est scalar fp_dev = `dev'
	est scalar fp_catz = `cz'
	est scalar fp_nx = `nbase'+1
	est scalar fp_df = `df'
	est scalar fp_rdf = `rdf'
	est scalar fp_N = `nobs'
	est local  fp_opts `options'
	est local  fp_t1t "Fractional Polynomial"
	est local  fp_pwrs `pwrs'
	est local  fp_wgt "`weight'"
	est local  fp_wexp "`exp'"
	est local  fp_xp `xp'
	est local  fp_base `base'
	est local  fp_rhs `rhs'
	est local  fp_depv `y'
	est local  fp_fvl `xp' `base'
	est scalar fp_dist = `dist'

	est scalar f_b0 = `b0'
	est scalar f_b1 = `b1'

	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 */

	est local  fp_fprp "no"
	est scalar fp_glmw = `glmwarn'
	est scalar fp_srch = `search'
	est scalar fp_sfac = `scalfac'
	est scalar fp_shft = `shift'
	est local  fp_xpx `kx'

	global S_E_cmd `cmd'			/* double save */
	est local  fp_cmd "fracpoly"
	est 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       Deviance     " _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 dev = e(fp_dlin)
			local rsd = e(fp_slin)
			local d=e(fp_d0)-`dev'
			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 dev = e(fp_d`m')
			local rsd = e(fp_s`m')
			local d=`devlast'-`dev'
			local n1 = cond(`m'==1, 1, 2)
			local n2=e(fp_rdf)-2*`m'-`cz'
			local pwr ${S_E_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 `dev' _cont
		if `normal' {
			di in ye _skip(5) %8.0g `rsd' _cont
		}
		di in ye _skip(3) %7.3f e(fp_dlin)-`dev' /*
	 	*/ %9.3f `P' _skip(2) "`pwr'"
		local devlast `dev'
		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]"
	}
	if e(fp_glmw) { 
		di in bl /*
 	*/ "[warning: GLM deviances are unscaled---ignore P-values]" }
end