glm_6.ado

来自「是一个经济学管理应用软件 很难找的 但是经济学学生又必须用到」· ADO 代码 · 共 831 行 · 第 1/2 页

*! version 4.1.9  20dec2004
program define glm_6, eclass byable(recall)
	version 6.0, missing
	if replay() {
		if _by() { error 190 }
		Playback `0'
		exit
	}

	syntax varlist [if] [in] [aw fw iw] [, LEvel(cilevel) EForm /*
		*/ LTol(real -1) ITerate(int 0) /*
		*/ Family(string) Link(string) DISP(real 1) Scale(string) /*
		*/ FRVars Offset(varname numeric) LNOffset(varname numeric) /*
		*/ noCONStant INIt(varname numeric) noLOg ]


	marksample touse 		/* must be done early */

	if `"`log'"'!="" { local iflog `"*"' }
	else	local iflog "noi"

	if `"`constan'"'!="" { local cons "nocons" }
	local small 1e-6
	if `ltol'<0 { local ltol 1e-6 } /* was 1e-8 in glmr (sg22)-too small */
	if `iterate'<=0 { local iterate 50 }
/*
	Validate family and set up link indicator
	(MapFL is why touse must be resolved so soon)
*/
	MapFL `"`family'"' `"`link'"' 
	local family 	`"`r(family)'"'
	local link 	`"`r(link)'"'
	local pow 	`r(power)'
	local scale1 	`r(scale)' /* 1 for fams with default scale param 1 */
	local m 	`r(m)'
	local mfixed 	`r(mfixed)'
	local k 	`r(k)'

	local bernoul = `mfixed' & `"`m'"'=="1" /* Bernoulli dist */
	local reg = `"`family'"'=="gau" & `"`pow'"'=="1"
	if `reg' { local iterate 1 }
/*
	Values of `scale' used:
	-------------------------------------------
	"`scale'"    known scale      unknown scale
	-------------------------------------------
	"" (default)	"" => 1		"" => dispc
	dev		-1 => dispd	-1 => dispd
	x2		 0 => dispc	 0 => dispc
	#		 #		 #
	-------------------------------------------
*/
	if `"`scale'"'!="" {
		if `"`scale'"'=="x2" { local scale 0 }
		else if `"`scale'"'=="dev" { local scale -1 }
		else {
			capture confirm number `scale'
			if _rc {
				di in red "invalid scale()"
				exit 198
			}
		}
	}
	local efopt `"`eform'"'
	if `"`eform'"'!="" {
		MapEform `family' `"`link'"'
		local eform `"`r(eform)'"'
	}
 	tempvar mu eta W z V dev wt
	quietly {
		tokenize `varlist'
		local y `"`1'"'
		mac shift
/*
	Deal with missing values of binomial denominator variable (m).
*/
		if !`mfixed' { local mmiss `"`m'"' }
/*
	Offset / log offset
*/
		if `"`lnoffse'"'!="" {
			if `"`offset'"'!="" { error 198 }
			tempvar offvar
			local offstr `"ln(`lnoffse')"'
			gen double `offvar' = `offstr'
		}
		else if `"`offset'"'!="" {
			local offvar `"`offset'"'
			local offstr `"`offset'"'
		}
		if `"`offstr'"'!="" {
			local offopt `", offset `offstr'"'
		}

		markout `touse' `offvar' `mmiss'

		sum `y' if `touse'	/* purposely not weighted */
		if r(N)==0 { noisily error 2000 } 
		if r(N)==1 { noisily error 2001 }
		if r(min)==r(max) & `"`offstr'"'=="" {
			noi di _n in bl "outcome does not vary"
			exit 2000
		}

		local nobs = r(N) /* nobs reset below when fweights */
		if `"`weight'"'!="" { 
			gen double `wt' `exp' if `touse'
			if `"`weight'"'=="aweight" {
				summ `wt'
				replace `wt' = `wt'/r(mean)
			}
			else if `"`weight'"'=="fweight" {
				summ `touse' [fw=`wt'] if `touse'
				local nobs = round(r(N),1)
			}
		}
		else gen byte `wt' = `touse' if `touse'

/* check family(binomial) for sensibility of dependent variable */

	if `"`family'"'=="bin" {
		capture assert `y'>=0 if `touse'
		if _rc { 
			di in red `"dependent variable `y' has negative values"'
			exit 499
		}
		if `mfixed'==0 { 
			capture assert `m'>0 if `touse'		/* sic, > */
			if _rc { 
				di in red `"`m' has nonpositive values"'
				exit 499
			}
			capture assert `m'==int(`m') if `touse'
			if _rc {
				di in red `"`m' has noninteger values"'
				exit 499
			}
		}
		capture assert `y'<=`m' if `touse'
		if _rc {
			di in red `"`y' > `m' in some cases"'
			exit 499
		}
	}
/*
	Initialization: tolerance, iterations, mean (mu), working vectors
*/
		if `"`init'"'!="" {
			gen double `mu' = `init' if `touse'
		}
		else {
			sum `y' [aw=`wt'] if `touse'
			gen double `mu' = (`y'+r(mean))/(`m'+1)
		}
		gen double `W' = . in 1
		gen double `z' = . in 1
		gen double `V' = . in 1
		gen double `dev' = . in 1
/*
	Initialise linear predictor
*/
		if `"`link'"'!="pow" {
			if `"`family'"'=="bin" {
				if `"`init'"'=="" {
					replace `mu' = `m'*(`y'+.5)/(`m'+1) /*
					 */ if `touse'
				}
				if `"`link'"'=="l" {
					gen double `eta' = ln(`mu'/(`m'-`mu'))
				}
				else if `"`link'"'=="p" {
					gen double `eta' = invnorm(`mu'/`m')
				}
				else if `"`link'"'=="c" {
					gen double `eta' = ln(-ln(1-`mu'/`m'))
				}
				if `"`link'"'=="opo" {
					gen double `eta' = /*
					 */  cond(abs(`pow')<`small', /*
					 */ ln(`mu'/(`m'-`mu')), /*
					 */ ((`mu'/(`m'-`mu'))^`pow'-1)/`pow')
				}
			}
			if `"`family'"'=="gau" {
				gen double `eta' = `mu'
			}
			else if `"`family'"'=="poi" {
				gen double `eta' = ln(`mu')
			}
			else if `"`family'"'=="gam" {
				gen double `eta' = 1/`mu'
			}
			else if `"`family'"'=="ivg" {
				gen double `eta' = 1/`mu'^2
			}
			else if `"`family'"'=="nb" {
				gen double `eta' = -ln(1 + 1/(`k'*`mu'))
			}
		}
		else {
			if abs(`pow')<`small' {
				gen double `eta' =  ln(`mu')
			}
			else gen double `eta' = `mu'^`pow'
		}
/*
	Local scoring algorithm: IRLS
*/
		tempname newdev oldev Wscale
		scalar `oldev' =  0
		scalar `newdev' =  1  
		local i 1
		`iflog' di
		while abs(`newdev'-`oldev')>`ltol' & `i'<=`iterate' {
			scalar `oldev' = `newdev'
/*
	Get 1/(d(eta)/d(mu)) and iterative weights, W.
	(Up to a constant, these are the same for canonical link functions.)
	Also calculate the adjusted independent variable, z.
*/
			CalcV `touse' `y' `family' `"`link'"' `pow' `m' /*
			 */ `disp' `k' `eta' `mu' `W' `z' `V' `"`offvar'"'
/*
	Weighted regression
		`Wscale' added by wms 10 Apr 1995
*/
                        summarize `W' [aw=`wt']
                        scalar `Wscale' = r(mean)
                        regress `z' `*' [iw=`W'*`wt'/`Wscale'], mse1 `cons'
/*
	Get linear predictor
*/
			cap drop `eta'
			_predict double `eta' if `touse' 
			if `"`offstr'"'!="" { replace `eta' = `eta'+`offvar' }
/*
	Get inverse link
*/
			est local family `"`family'"'
			est local link `"`link'"'
			global S_E_fam `"`family'"'   /* double save */
			global S_E_link `"`link'"'
			_crcglil `eta' `mu' `pow' `m' `k' `bernoul'
			count if `mu'>=. & `eta' <  . & `touse'
			if r(N) != 0 {
				noi di in red "estimates diverging"
				exit 430
			}
/*
	Get weighted squared deviance contributions (residuals)
*/
			_crcgldv `y' `family' `bernoul' `m' `k' `mu' `wt' `dev'
/*
	Get deviance (note: analytic weights are already normalized).
*/
			replace `z' = sum(`dev')
			scalar `newdev' = `z'[_N]/`disp'
			`iflog' di in gr `"Iteration `i' : deviance = "' /*
			 */ in ye %9.4f `newdev'
			local i = `i'+1
		}
		local conrc = cond(abs(`newdev'-`oldev')>`ltol' & !`reg',430,0)
		if `conrc' { 
			noisily di in ye "(convergence not achieved)"
		}
/*
	Get weighted squared contributions (residuals) to Pearson X2
*/
		replace `z' = sum(`wt'*(`y'-`mu')^2/`V')
		local chisq = `z'[_N]/`disp'
		local df = `nobs'-e(df_m)-(`"`cons'"'=="")
		local dispc = `chisq'/`df'
		local dispd = `newdev'/`df'
	}
/*
	Store results
*/

	DescFL `"`family'"' `"`link'"' `pow' `bernoul' `k' `m'
	local o `"`r(dist)', `r(link)'`offopt'"'

	if `"`scale'"'=="" {	/* default */
		local scale `scale1'
		if `scale1' { local delta 1 }
		else local delta `dispc'
	}
	else if `scale'==0 {	/* Pearson X2 scaling */
		local delta `dispc'
		if `scale1' {
			local cd "square root of Pearson X2-based dispersion"
		}
	}
	else if `scale'==-1 {	/* deviance scaling */
		local delta `dispd'
		local cd "square root of deviance-based dispersion"
	}
	else {			/* user's scale parameter */
		local delta `scale'
		if !`scale1' | (`scale1' & `scale'!=1) {
			local cd `"dispersion equal to square root of `delta'"'
		}
	}
	if !`scale1' | (`scale1' & `scale'!=1) {
		if `scale1' { local dof 100000 }
		else local dof `df'
		if `delta'>=. { local zapse "yes" }
                else scalar `Wscale' = `Wscale'/`delta' /* scale to `delta' */
	} 

	tempname b V
	if `reg' {
		local dofopt `"dof(`df')"'
	}
	mat `b' = get(_b)
	mat `V' = get(VCE)
        scalar `Wscale' = 1/`Wscale'
        mat `V' = `Wscale'*`V' /* get rid of `Wscale' scaling */
	if `"`zapse'"'=="yes" { 
		local i 1
		while `i'<=rowsof(`V') {
			mat `V'[`i',`i'] = 0 
			local i=`i'+1
		}
	}
	tempvar mysamp
	qui gen byte `mysamp' = `touse'
	est post `b' `V', depname(`y') obs(`nobs') `dofopt' esample(`mysamp')

/* results can be saved now */
/* double save in S_E_ and e() */
	if `"`cd'"'!="" {
		est local msg_1 `"(Standard errors scaled using `cd')"'
		global S_E_msg1 `"(Standard errors scaled using `cd')"'
	}
	if `reg' {
		est local msg_2 "(Model is ordinary regression, use regress instead)"
		global S_E_msg2 "(Model is ordinary regression, use regress instead)"
	}
	if `disp'!=1 {
		est local msg_3 `"(Quasi-likelihood model with dispersion `disp')"'
		global S_E_msg3 `"(Quasi-likelihood model with dispersion `disp')"'
	}
	if `"`frvars'"'!="" {
		cap drop _mu
		cap drop _eta
		cap drop _dres
		gen _dres = sign(`y'-`mu')*sqrt(`wt'*`dev'/`delta') if `touse'
		rename `mu' _mu
		rename `eta' _eta
		lab var _mu `"E(`y')"'
		lab var _eta "Linear predictor"
		lab var _dres "Scaled deviance-residuals"
	}

/*
	Save results -- new  gould
*/

/* double save in S_E_ and e() (actually triple save in S_#) */
	est scalar df_pear = `df'
	est scalar N = `nobs'
	est scalar chi2 = `chisq'
	est scalar deviance = `newdev'
	est scalar dispersp = `dispc'
	est scalar dispers = `dispd'

	global S_E_rdf `df'
	global S_E_nobs `nobs'
	global S_E_chi2 `chisq'
	global S_E_dev = `newdev'
	global S_E_dc `dispc'
	global S_E_dd `dispd'

	global S_1 `nobs'
	global S_2 `df'
	global S_3 = `newdev'
	global S_4 `chisq'

/*
	Save more results (for glmrpred) -- new Royston 4/94.
*/

/* more double saves in S_E_ and e() */
	est local wtype `"`weight'"'
	est local wexp `"`exp'"'
	est local family `"`family'"'
	est local link `"`link'"'
	est local title_fl `"`o'"'
	est scalar bernoul = `bernoul'
	est local m `"`m'"'
	est local depvar `"`y'"'
	est scalar power = `pow'
	est local offset `"`offstr'"'	        /* offset/lnoffset */
	est scalar lnoffset = `"`lnoffse'"'!="" /* 1 if lnoffset, 0 otherwise */
	est local eform `"`eform'"'
	est scalar k = `k'  
	est scalar disp = `disp'
	est scalar delta = `delta'
	est scalar rc = `conrc'

	global S_E_vl `"`varlist'"'
	global S_E_if `"`if'"'
	global S_E_in `"`in'"'
	global S_E_wgt `"`weight'"'
	global S_E_exp `"`exp'"'
	global S_E_fam  `"`family'"'
	global S_E_link `"`link'"'
	global S_E_flo  `"`o'"'
	global S_E_ber  `bernoul'
	global S_E_m    `"`m'"'
	global S_E_depv `"`y'"'
	global S_E_pow  `"`pow'"'
	global S_E_off  `"`offstr'"'		/* offset/lnoffset */
	global S_E_lno  = `"`lnoffse'"'!=""	/* 1 if lnoffset, 0 otherwise */
	global S_E_ef   `"`eform'"'
	global S_E_k    `"`k'"'