mfp.ado

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

*! version 1.1.10  05nov2004
program define mfp, eclass sortpreserve
version 8
if "`1'" == "" | "`1'"=="," {
	if "`e(fp_cmd2)'"!="mfp" {
		error 301
	}
	FracRep "fractional polynomial" "  df  " "Powers"
	exit
}
gettoken cmd 0 : 0
frac_chk `cmd' 
if `s(bad)' {
	di in red "invalid or unrecognized command, `cmd'"
	exit 198
}
/*
	dist=0 (normal), 1 (binomial), 2 (poisson), 3 (cox), 4 (glm),
	5 (xtgee), 6(ereg/weibull), 7(streg,stcox).
*/
local dist `s(dist)'
local glm `s(isglm)'
local qreg `s(isqreg)'
local xtgee `s(isxtgee)'
local normal `s(isnorm)'

global MFpdist `dist'

* Disentangle
GetVL `0'
local 0 `s(nought)'

syntax [if] [in] [aw fw pw iw] , /* 
 */ [ ADDpowers(str) ADJust(str) AIC ALpha(str) ALL CATzero(str) DF(str) /*
 */ DFDefault(int 4) CYCles(int 5) DEAD(str) FIXpowers(str) FP01 noSCAling /*
 */ noCONStant POwers(str) SELect(str) SEQuential noTESTLinear XOrder(str) /*
 */ XPowers(str) ZERo(str) * ]
frac_cox "`dead'" `dist'
/*
	Process options
*/

local regopt `options' `constant'
if "`dead'"!="" {
	local regopt "`regopt' dead(`dead')"
}
if "`powers'"==""  {
	local powers "-2,-1,-.5,0,.5,1,2,3"
}
if "`addpowers'"!="" {
	local fpopt "`fpopt' addp(`addpowers')"
}
if "`fixpowers'"!="" {
	local fpopt "`fpopt' fixp(`fixpowers')"
}
if "`aic'"!="" { /* aic selection for vars and functions */
	if "`alpha'`select'"!="" {
		noi di as err "alpha() and select() invalid with aic"
		exit 198
	}
	local alpha -1
	local select -1
}/*
	Check for missing values in lhs, rhs and model vars.
*/
tempvar touse
quietly {
	marksample touse
	markout `touse' $MFP_cur `dead'
	frac_wgt "`exp'" `touse' "`weight'"
	local wgt `r(wgt)'				/* [`weight'`exp'] */
	count if `touse'
	local nobs = r(N)
}
/*
	Detect collinearity among covariates, and fail if found.
*/
local ncur: word count $MFP_cur
_rmcoll $MFP_cur `if' `in' [`weight' `exp'], `constant'
local ncur2: word count `r(varlist)'
if `ncur2'<`ncur' {
	local ncoll=`ncur'-`ncur2'
	if `ncoll'>1 {
		local s ies
	}
	else local s y
	di in red `ncoll' " collinearit`s' detected among covariates"
	exit 198
}
/*
	Rearrange order of variables in varlist
*/
if "`xorder'"=="" {
	local xorder "+"
}
/*
	Apply fracord to data in first imputation, using micombine to get param estimates
*/
FracOrd `wgt' if `touse', order(`xorder') `regopt' cmd(`cmd')
local nx $MFP_n	/* number of clusters, <= number of predictors */
local lhs $MFP_dv
local i 1
while `i'<=`nx' {
/*
	Store original order and reverse order
	of each RHS variable/variable set
*/
	local r`i' `s(ant`i')'
	local i = `i'+1
}
/*
	Initialisation.
*/
local i 1
while `i'<=`nx' {
	local x ${MFP_`i'}
	local nx`i': word count `x'
	local alp`i' .05	/* default FP selection level */
	local h`i' `x'		/* names of H(xvars) 	*/
	local n`i' `x'		/* names of xvars 	*/
	local po`i' 1 		/* to be final power	*/
	local sel`i' 1		/* default var selection level */
/*
	Remove old I* variables
*/
	if `nx`i''==1 {
		frac_mun `n`i'' purge
	}
	local i=`i'+1
}
/*
	Adjustment
*/
FracAdj "`adjust'" `touse'
local i 1
while `i'<=`nx' {
	if "`r(adj`i')'"!="" {
		local adj`i' adjust(`r(adj`i')')
	}
	local uniq`i'=r(uniq`i')
	local i=`i'+1
}
/*
	Set up degrees of freedom for each variable
*/
if "`df'"!="" {
	FracDis "`df'" df 1 .
	local i 1
	while `i'<=`nx' {
		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.
*/
local i 1
while `i'<=`nx' {
	if `nx`i''>1 {
		* over-ride all suggestions that df>1 for grouped vars
		local df`i' 1
	}
	else {
		if "`df`i''"=="" {
			if `uniq`i''<=3 {
				local df`i' 1
			}
			else if `uniq`i''<=5 {
				local df`i'=min(2,`dfdefault')
			}
			else local df`i' `dfdefault'
		}
	}
	local i=`i'+1
}
/*
	Set up FP selection level (alpha) for each variable
*/
if "`alpha'"!="" {
	FracDis "`alpha'" alpha -1 1
	local i 1
	while `i'<=`nx' {
		if "${S_`i'}"!="" {
			local alp`i' ${S_`i'}
			if `alp`i''<0 {
				local alp`i' -1  /* AIC */
			}
		}
		local i=`i'+1
	}
}
/*
	Set up selection level for each variable
*/
if "`select'"!="" {
	FracDis "`select'" select -1 1
	local i 1
	while `i'<=`nx' {
		if "${S_`i'}"!="" {
			local sel`i' ${S_`i'}
			if `sel`i''<0 {
				local sel`i' -1  /* AIC */
			}
		}
		local i=`i'+1
	}
}
/*
	Rationalise select() and alpha() in cases of aic
*/
local i 1
while `i'<=`nx' {
	if `sel`i''==-1 & `alp`i''!=1 {
		local alp`i' -1
	}
	if `alp`i''==-1 & `sel`i''!=1 {
		local sel`i' -1
	}
	local i=`i'+1
}
/*
	Individual FP powers for variables.
*/
if "`xpowers'"!="" {
	FracDis "`xpowers'" xpowers
	local i 1
	while `i'<=`nx' {
		if "${S_`i'}"!="" & `nx`i''==1 {
			local xpow`i' ${S_`i'}
		}
		local i=`i'+1
	}
}
/*
	Vars with zero option
*/
if "`zero'"!="" {
	tokenize `zero'
	while "`1'"!="" {
		FracIn `1'
		local i `s(k)'
		if `nx`i''==1 {
			local zero`i' "zero"
		}
		mac shift
	}
}
/*
	Vars with catzero option
*/
if "`catzero'"!="" {
	tokenize `catzero'
	while "`1'"!="" {
		FracIn `1'
		local i `s(k)'
		if `nx`i''==1 {
			local catz`i' "catzero"
			local zero`i' "zero"		/* catzero implies zero */
		}
		mac shift
	}
}
/*
	Reserve names for H(predictors) by creating a dummy variable
	for each predictor which potentially needs transformation.
*/
local i 1
while `i'<=`nx' {
	if `df`i''>1 | "`zero`i''`catz`i''"!="" {
		frac_mun `n`i''
		local stub`i' `s(name)'
		qui gen byte `stub`i''_1=.
	}
	local i=`i'+1
}
/*
	Build FP model.
	`r*' macros present predictors according to FracOrd ordering,
	e.g. i=1, r`i'=3 means most sig predictor is third in user's xvarlist.
*/
local it 0
local initial 1
local stable 0 /* convergence flag */
while !`stable' & `it'<=`cycles' {
	local it = `it'+1
	local pwrs
	local rhs1
	local stable 1 /* later changed to 0 if any power or status changes */
	local lastch 0 /* becomes index of last var which changed status */
	local i 1
	while `i'<=`nx' {
		local r `r`i''
		local ni `n`r''
		local dfi df(`df`r'')
/*
	Build up RHS2 from the i+1th var to the end 
*/
		local rhs2
		local j `i'
		while `j'<`nx' {
			local j = `j'+1
			local rhs2 `rhs2' `h`r`j'''
		}
		if `initial' {
			if "`rhs2'"!="" {
				local fixed "base(`rhs2')"
			}
			else local fixed
			qui FracSel `cmd' `lhs' `ni' `wgt' if `touse', /*
			 */ df(1) `fixed' select(1) `regopt'
			local dev=r(dev)
			di in gr _n /*
			 */ "Deviance for model with all terms " /*
			 */ "untransformed = " in ye %9.3f r(dev) in gr ", " /*
			 */ in ye `nobs' in gr " observations"
		}
		if "`rhs1'`rhs2'"!="" {
			local fixed "base(`rhs1' `rhs2')"
		}
		else local fixed
/*
	Vars with df(1) are straight-line
*/
		local pvalopt "alpha(`alp`r'') select(`sel`r'')"
		if `i'==1 {
			di
		}
		if `df`r''==1 {
			local pw
			local fpo
		}
		else {
			if "`xpow`r''"!="" {
				local pw "powers(`xpow`r'')"
			}
			else local pw "powers(`powers')"
			local fpo `fpopt'
		}
		if `df`r''==1 & `sel`r''==1 {	/* var is included anyway */
			local rhs1 `rhs1' `h`r''
			di in gr "[`ni' included with 1 df in model]" _n
		}
		else {
			if "`stub`r''"!="" {
				local n name(`stub`r'')
			}
			else local n
			FracSel `cmd' `lhs' `ni' `wgt' if `touse', `dfi' /*
			 */ `pw' `zero`r'' `catz`r'' `fixed' `h' /*
			 */ `fpo' `regopt' `pvalopt' `n' `sequential' `testlinear'
			local h`r' `r(rhs)'
			local dev=r(dev)
			local p `r(pwrs)'
			if "`p'"!="`po`r''" {
				if `nx'>1 {
					local stable 0
				}
				local po`r' "`p'"
				local lastch `i'
			}
			if "`pwrs'"=="" {
				local pwrs "`p'"
			}
			else local pwrs "`pwrs',`p'"
			if "`h`r''"!="" {
				local rhs1 `rhs1' `h`r''
			}
		}
		if `initial' {
			local h "nohead"
			local initial 0
		}
		local i = `i'+1
	}
	if `lastch'==1 {
		local stable 1
	} /* 1 change only, at i=1 */
	if !`stable' {
		di in smcl in gr "{hline 70}" _n "End of Cycle " in ye `it' in gr /*
		 */ ": deviance =   " in ye %9.3f `dev' _n in gr "{hline 70}"
	}
}
if `nx'>1 {
	local s
	if `it'!=1 {
		local s "s"
	}
	if !`stable' {
		di _n in gr "No convergence" _cont
	}
	else di _n in gr /*
	 */ "Fractional polynomial fitting algorithm converged" _cont
	di in gr " after " in ye `it' in gr " cycle`s'."
}
if `stable' {
	di _n in gr "Transformations of covariates:" _n
}
/*
	Remove variables left behind by frac_154
*/
local i 1
while `i'<=`nx' {
	if "`stub`i''"!="" {
		cap drop `stub`i''*
	}
	local i=`i'+1
}
/*
	Store results
*/
local finalvl	/* predictors in final model */
local i 1
while `i'<=`nx' {
	local p `po`i''
	local x `n`i''
	local z `zero`i''
	local c `catz`i''
	local a `adj`i''
/*
	Create FP vars as necessary, with new unique names
*/
	if trim("`p'")!="." {
		local nxi: word count `x'
		if (trim("`p'")=="1" & "`z'`c'`a'"=="") | (`nxi'>1) {
			local namex `x'
		}
		else {
			frac_mun `x'
			local vn `s(name)'
			fracgen `x' `p' if e(sample), `all' `scaling' /*
			 */ `z' `c' `a' name(`vn')
			local namex `r(names)'
		}
		local finalvl `finalvl' `namex'
	}
/*
	Save stuff for FracRep to display
*/
	local npows: word count `p'
	local fd`i'=2*`npows'-(`npows'==1 & "`p'"=="1")+("`c'"!="")
	local id`i'=`df`i''+("`c'"!="")
	local i=`i'+1
}
/*
	Estimate final (conditionally linear) model.
*/
quietly `cmd' `lhs' `finalvl' `wgt' if e(sample), `regopt'
/*
	Store results for each predictor
	(requires expansion if grouped variables present)
*/
global S_1 `finalvl'
global S_2 `dev'
local i 1
local nx2 0	/* number of predictors after expansion of groups (if any) */
while `i'<=`nx' {
	ereturn scalar Fp_fd`i'=`fd`i''	/* final degrees of freedom */
	ereturn scalar Fp_id`i'=`id`i''	/* initial degrees of freedom */
	ereturn scalar Fp_al`i'=`alp`i''	/* FP selection level */
	ereturn scalar Fp_se`i'=`sel`i''	/* var selection levl */
	ereturn local Fp_k`i' `po`i'' 	/* "powers" for ith predictor */

	tokenize `n`i''
	while "`1'"!="" {
		local nx2=`nx2'+1
		ereturn local fp_x`nx2' `1'		/* name of ith predictor in user order */
		ereturn local fp_k`nx2' `po`i'' 	/* "powers" for ith predictor */
		if "`catz`i''"!="" {
			ereturn local fp_c`nx2' 1
		}
		mac shift
	}
	local i=`i'+1