glm.ado

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

*! version 5.6.13  01may2005
program glm, eclass byable(onecall) prop(ml_score swml)
	local version : di "version " string(_caller()) ", missing:"
	if _by() {
		local BY `"by `_byvars'`_byrc0':"'
	}
	local vceopts bootopts(noheader notable) jkopts(noheader notable) 
	`BY' _vce_parserun glm, `vceopts' robustok multivce : `0'
	if "`s(exit)'" != "" {
		// -bootstrap- and -jackknife- are not aware of -vfactor()-
		if e(vf) != 1 {
			tempname V
			matrix `V' = e(vf)*e(V)
			version 9: ereturn repost V = `V'
		}
		version 9: ereturn local clustvar `"`e(cluster)'"'
		version 9: ereturn local cluster
		quietly syntax [anything] [fw aw iw pw] [if] [in] [,	///
			noHEADer noTABLE Level(passthru) eform * ]
		version 7: Display, `header' `table' `level' `eform'
		exit
	}

	version 7, missing
	if replay() {
		if _by() {
			error 190
		}
		if "`e(cmd)'" != "glm" {
			error 301
		}
		if "`e(predict)'" == "glm_p" {
			glm_6 `0'
			exit
		}
		Display `0'
		error `e(rc)'
		exit
	}
	if _caller() < 7 {
		mac drop SGLM_*
		capture noi `BY' glm_6 `0'
		mac drop SGLM_*
		exit _rc
	}
	mac drop SGLM_running SGLM_nonstan
	capture noisily `version' `BY' glm_7 `0'
	if _rc {
		if "$SGLM_running" ~= "" & "$SGLM_nonstan" ~= "" {
			di as err /*
		*/"data not suitable for non-standard family-link combination"
			mac drop SGLM_*
			exit 459
		}
		else {
			mac drop SGLM_*
			exit _rc
		}
	}
	if "$SGLM_running" == "" & "$SGLM_nonstan" == "" {
		mac drop SGLM_*
	}
end

program glm_7, eclass byable(recall) sort
	version 7, missing

/* Parse. */

	if _caller() < 9 {
		local iteropt ITERate(integer 50)
	}
	else	local iteropt ITERate(integer `c(maxiter)')

	local dltol -1

	syntax varlist(numeric ts) [fw aw pw iw]	///
		[if] [in] [,				///
		EForm					///
		FROM(string)				///
		LEvel(cilevel)				///
		OFFset(varname numeric)			///
        	LNOFFset(varname numeric)		///
		noCONstant				///
		Robust					///
		CLuster(varname)			///
        	Family(string)				///
		Link(string)				///
		SEARCH					///
		IRLS					///
        	LTOLerance(real `dltol')		///
        	SCore(string)				///
		noLOg					///
		noTABLE					///
		noDISPLAY				///
		FISHER(integer -1)			///
		noDOTS					///
		`iteropt'				///
		TRAce					///
		SCAle(string)				///
		MU(varname)				///
        	INIt(varname)				///
        	T(string)				///
		noHEADer				///
		DISP(real 1)				///
		VFactor(real 1)				///
        	EXPosure(varname numeric)		///
		*					///
	]

	// NOTE: GetVCE creates one or more of the following local macros:
	// 	bstrap
	// 	brep
	// 	jknife
	// 	jknife1
	// 	hac
	//	oim
	//	opg
	//	unbiased
	//	vce
	// 	options
	GetVCE, `irls' `options'

	if "`family'"!= "" {
		local argfam `"family(`family')"'
	}
	if "`link'" != "" {
		local arglink `"link(`link')"'
	}

	if "`init'" != "" & "`mu'" != "" & "`init'" != "`mu'" {
		noi di as err "mu() and init() should use the same varname"
		exit 198
	}

	if "`init'" != "" {
		local init "`init'"
	}
	else	local init "`mu'"

	if "`dots'" == "nodots" {
		local dots ""
	}
	else    local dots "dots"

	if "`unbiased'" != "" {
		local robust "robust"
	}

	if `vfactor' <= 0.0 {
		di as err "vfactor() must be positive"
		exit 198
	}

	if `iterate' < 0 {
		di as err "iterate() must be positive"
		exit 198
	}

	mlopts mlopts, `options'
	if `ltolerance' == -1 {
		local ltol 1e-6
	}
	else {
		local ltol `ltolerance'
		if "`irls'" == "" {
			local mlopts `"`mlopts' ltol(`ltol')"'
		}
	}
	if `ltol'>=1 | `ltol'<0 {
		di as err "ltolerance() must be in [0,1)"
		exit 198
	}

/* Check syntax. */

	if `"`score'"'!="" {
		if "`irls'" != "" {
			di as err "deviance scores are calculated " /*
			*/ "postestimation via predict, score"
			exit 198
		}
		confirm new variable `score'
		local nword : word count `score'
		if `nword' > 1 {
			di as err "score() must contain the name of only " /*
			*/ "one new variable"
			exit 198
		}
		tempvar scvar
		local scopt score(`scvar')
	}
	if `"`cluster'"'!="" {
		if "`jknife'`bstrap'" == "" {
			local robust "robust"
		}
		local clopt cluster(`cluster')
	}

	if "`bstrap'" == "" {
		if !inlist("`brep'","","-1") {
			di as err "brep() only valid with bstrap"
			exit 198
		}
		else	local brep -1
	}

	if "`bstrap'" != "" {
		if `brep' == -1 {
			local brep = 199
		}
		if `brep' <= 40 {
			di as err "bootstrap replications must be > 40"
			exit 198
		}
	}

	local tsops = index("`varlist'", ".")
	if (`tsops' & "`bstrap'`jknife'`jknife1'" != "") {
		di as error "cannot use time-series operators with " _c
		di as error "bstrap, jknife, or jknife1"
		exit 101
	}

	if `fisher' != -1 {
		if "`irls'" != "" {
			di as err /*
		*/ "fisher() only valid with Newton-Raphson optimization"
			exit 198
		}
		if `fisher' <= 0 {
			di as err "fisher() must be positive"
			exit 198
		}
	}

	if `disp' < 0 {
		di as err "disp() must be positive"
		exit 198
	}

        if `"`scale'"'!="" {
		if "`robust'`opg'`hac'`jknife1'`jknife'`bstrap'" != "" {
			di as err "cannot use scale() with alternate variances"
			exit 198
		}
                if `"`scale'"'=="x2" {
			local scale 0
		}
                else if `"`scale'"'=="dev" {
			local scale -1
		}
                else {
                        capture confirm number `scale'
                        if _rc {
                                di as err "invalid scale()"
                                exit 198
                        }
			if `scale' <= 0 {
				di as err "scale(#) must be positive"
				exit 198
			}
                }
        }

	if "`scale'" == "-1" & "`irls'" == "" {
		di as err "scale(dev) only allowed with irls option"
		exit 198
	}

	if `disp' != 1 & "`irls'" == "" {
		di as err "disp() only allowed with irls option"
		exit 198
	}

	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)'

	if "`link'" == "" | "`family'" == "" {
		di as err "incomplete specification of family() and link()"
		exit 198
	}

	if "`offset'" != "" & "`lnoffset'" != "" {
		di as err ///
		"only one of offset() or lnoffset() can be specified"
		exit 198
	}
	if "`offset'" != "" & "`exposur'" != "" {
		di as err ///
		"only one of offset() or exposure() can be specified"
		exit 198
	}
	if "`lnoffset'" != "" & "`exposur'" != "" {
		di as err ///
		"only one of lnoffset() or exposure() can be specified"
		exit 198
	}
	if "`exposur'" != "" {
		local lnoffset `exposur'
		local exposur
	}
	if "`constant'" != "" {
		local nvar : word count `varlist'"
		if `nvar' == 1 {
			di as err ///
"independent variables required with noconstant option"
			exit 100
		}
	}

/* Mark sample except for offset/exposure. */

	marksample touse
	if "`cluster'" != "" {
		markout `touse' `cluster', strok
	}

	if `"`hac'"' == "" {
		if `"`t'"' ~= "" {
			di as err "t() only valid with vce(hac ...)"
			exit 198
		}
	}

	if "`hac'" != "" {	/* jknife1 and jknife removed -- rgg */
		xt_tis `t'
		local tvar `"`s(timevar)'"'
		markout `touse' `tvar'
		sort `touse' `tvar'
		cap assert `tvar'[_n-1] != `tvar' if `touse' & /*
		*/ (`touse'[_n-1]==1)
		if _rc {
			di as err "repeated time values in sample"
			exit 451
		}
	}

	if "`clopt'" != "" & "`vce'" != "" {
		opts_exclusive "cluster() vce()"
	}

/* Process offset/exposure. */

	if "`lnoffset'" != "" {
		capture assert `lnoffset' > 0 if `touse'
		if _rc {
			di as err "exposure() must be greater than zero"
			exit 459
		}
		tempvar offset
		qui gen double `offset' = ln(`lnoffset')
		local offvar "ln(`lnoffset')"
	}

	if "`offset'" != "" {
		markout `touse' `offset'
		local offopt "offset(`offset')"
		if "`offvar'" == "" {
			local offvar "`offset'"
		}
	}

/* Count obs and check for negative values of `y'. */

	gettoken y xvars : varlist
	tsunab y : `y'
	local yname : subinstr local y "." "_"

	summarize `y' if `touse', meanonly
	if r(N) == 0 { error 2000 }
	if r(N) == 1 { error 2001 }
	local nobs = round(r(N),1)

	MapHAC `nobs' `hac'
	local hacnam `"`r(hacnam)'"'
	local haclag `"`r(haclag)'"'

	local fwt "1"
	local awt "1"
	tempvar wt
	if "`weight'" != "" {
		qui gen double `wt' `exp' if `touse'
		local awt "`wt'"
		if `"`weight'"'=="aweight" {
			summ `wt' , meanonly
			qui replace `wt' = `wt'/r(mean)
		}
		else if `"`weight'"'=="fweight" {
			if "`bstrap'`hac'" != "" {
				di as err "fweights not allowed"
				exit 198
			}
			if "`jknife'" != ""  & "`cluster'" != "" {
				sort `cluster' `tvar' `touse'
				cap by `cluster' : assert `wt' == `wt'[_n-1]/*
				*/ if `touse' & (`touse'[_n-1]==1)

				if _rc {
					noi di as err "weight must be " /*
						*/ "constant within `cluster'"
					exit 198
				}
			}
			summ `touse' [fw=`wt'] if `touse' , meanonly
			local nobs = round(r(N),1)
			local fwt "`wt'"
			local awt "1"
		}
		if "`weight'" == "pweight" {
			local robust "robust"
		}
	}
	else {
		qui gen byte `wt' = `touse' if `touse'
	}

	if "`robust'" != "" {
		if "`hacnam'`opg'`jknife'`jknife1'" != "" {
			di as err "robust only allowed with hessian matrices"
			exit 198
		}
	}

	if "`weight'" != "" {
		if "`weight'" != "fweight" {
			if "`hacnam'`opg'`jknife'`jknife1'" != "" {
				di as err ///
"only fweights (if any) are available for non-hessian, non-bootstrap matrices"
				exit 198
			}
		}
	}

/* Remove collinearity. */

	if "`display'" != "" & "`weight'" == "aweight" {
		qui _rmcoll `xvars' [iw`exp'] if `touse', `constant'
			/* aweights produce "sum of wgt" message,
			   which is not wanted for -nodisplay-
			*/
	}
	else	_rmcoll `xvars' [`weight'`exp'] if `touse', `constant'
	local xvars `r(varlist)'

/* Get initial values. */

	global SGLM_V	`"`family'"'		/* V(mu) program */
	global SGLM_L	`"`link'"'		/* g(mu) g^(-1)(eta) program */
	global SGLM_A	`"`ancilla'"'		/* Ancillary params program */
	global SGLM_y	`"`y'"'			/* dep varname */
	global SGLM_m	`"`m'"'			/* Binomial denominator */
	global SGLM_a	`"`k'"'			/* Nbinom alpha */
	global SGLM_p	`"`pow'"'		/* Power, or argument */
						/* for user-defined link */
	global SGLM_f   `"`fisher'"'		/* Number of EIM iterations */
	global SGLM_mu                          /* Set in $SGLM_V -1 call */
	global SGLM_s1  `"`scale1'"'		/* Phi=1 indicator */
	/* Set titles, set range restrictions for mu, check sensibility of y */

	$SGLM_V -1 `touse'
	$SGLM_L -1 `touse'

	global SGLM_running 1

        if "$SGLM_V" == "glim_v1" | "$SGLM_V" == "glim_v3" | /*
	*/ "$SGLM_V" == "glim_v4" | "$SGLM_V" == "glim_v5" {
		if "$SGLM_L" == "glim_l02" | "$SGLM_L" == "glim_l04" | /*
		*/ "$SGLM_L" == "glim_l05" | "$SGLM_L" == "glim_l06" | /*
		*/ "$SGLM_L" == "glim_l07" | "$SGLM_L" == "glim_l08" | /*
		*/ "$SGLM_L" == "glim_l12" {
			global SGLM_nonstan 1
		}
	}
	if "$SGLM_V" == "glim_v2" {
		if "$SGLM_L" == "glim_l04" {
			global SGLM_nonstan 1
		}
	}
	if "$SGLM_V" == "glim_v6" {
		if "$SGLM_L" == "glim_l02" |  /*
                */ "$SGLM_L" == "glim_l05" | "$SGLM_L" == "glim_l06" | /*
                */ "$SGLM_L" == "glim_l07" | "$SGLM_L" == "glim_l08" | /*
                */ "$SGLM_L" == "glim_l12" {
                        global SGLM_nonstan 1
                }
	}

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

	if "`offset'" != "" {
		local moffset = "-`offset'"
	}

	local reg = "$SGLM_V" == "glim_v1" & "$SGLM_L" == "glim_l01"
	/* Newton-Raphson optimization using ml commands */
	if "`irls'" == "" {
		tempvar  eta mu dmu v z W
		tempname Wscale b0
		if "`vce'" != "eim" {
			local oim "oim"
		}

		if "`from'" == "" & "`search'" == "" {
			quietly {
				if `"`init'"'!="" {
					gen double `mu' = `init' if `touse'
				}
				else {
					sum `y' [aw=`wt'] if `touse'
					if "$SGLM_V" != "glim_v2" {
						gen double `mu' = /*
							*/ (`y'+r(mean))/(`m'+1)
					}
					else {
						gen double `mu' = /*
							*/ `m'*(`y'+.5)/(`m'+1)
					}
				}
				capture drop `eta'
				$SGLM_L 0 `eta' `mu'
				$SGLM_V 1 `eta' `mu' `v'
				$SGLM_L 2 `eta' `mu' `dmu'
				gen double `z' = `eta' + /*
					*/ (`y'-`mu')/`dmu' `moffset'
				gen double `W' = `dmu'*`dmu'/`v'
				summ `W' [aw=`wt'], meanonly
				scalar `Wscale' = r(mean)
				cap noisily quietly regress `z' `xvars' /*
				*/  [iw=`W'*`wt'/`Wscale'], mse1 `constant'
				if e(N)>=. {exit 459}
				if _rc {exit _rc}
/*
		mat `b0' = e(b)
		local k = colsof(`b0')
		regress `z' [iw=`W'*`wt'/`Wscale'], mse1 `constant'

		capture glm `y', fam(`family') link(`link') irls iter(2)
*/

				capture drop `eta'
				capture drop `mu'
				predict double `eta', xb
				$SGLM_L 1 `eta' `mu'