gnbreg.ado

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

*! version 3.3.7  21jun2005
program define gnbreg, eclass sort prop(svyb svyj svyr irr ml_score)
	_vce_parserun gnbreg : `0'
	if "`s(exit)'" != "" {
		exit
	}

	version 6, missing
	if _caller() <= 5 {
		gnbreg_5 `0'
		exit
	}
	if replay() {
		if "`e(cmd)'"!="gnbreg" { error 301 }
		Display `0'
		error `e(rc)'
		exit
	}

/* Parse. */
	local awopt = cond(_caller()<9, "aw", "")
	syntax varlist(numeric) [fw `awopt' pw iw] [if] [in] [, IRr /*
	*/ FROM(string) Level(cilevel) OFFset(varname numeric) 		/*
	*/ Exposure(varname numeric) noCONstant Robust CLuster(varname) /*
	*/ SCore(string) noLOg LNAlpha(varlist numeric)			/*
	*/ CRITTYPE(passthru) * ]

	mlopts mlopts, `options'
	local cns `s(constraints)'
	local mlopts `mlopts' `crittype'
	local lrtest nolrtest  /* rgg: I'm leaving in the code for the */
			       /* lrtest for when the dist. becomes known */

/* Check syntax. */

	if "`score'" != "" { 
		local n : word count `score'
		if `n'==1 & substr("`score'",-1,1)=="*" { 
			local score = /*
			*/ substr("`score'",1,length("`score'")-1)
			local score `score'1 `score'2 
			local n 2
		}
		confirm new variable `score'
		if `n' != 2 {
			di in red "score() must contain the names of " /*
			*/ "two new variables"
			exit 198
		}
		local scname1 : word 1 of `score'
		local scname2 : word 2 of `score'
		tempvar scvar1 scvar2
		local scopt "score(`scvar1' `scvar2')"
	}
	if "`offset'"!="" & "`exposur'"!="" {
		di in red "only one of offset() or exposure() can be specified"
		exit 198
	}
	if "`constan'"!="" {
		local nvar : word count `varlist'
		if `nvar' == 1 {
			di in red "independent variables required with " /*
			*/ "noconstant option"
			exit 100
		}
	}

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

	marksample touse

	if `"`cluster'"'!="" {
		markout `touse' `cluster', strok
		local clopt cluster(`cluster')
		local robust robust
	}
	if `"`weight'"' == "pweight" {
		local robust robust
	}

	markout `touse' `lnalpha'

/* Process offset/exposure. */

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

	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

	if "`weight'"!="" {
		if "`weight'"=="fweight" {
			local wt `"[fw`exp']"'
		}
		else	local wt `"[aw`exp']"'
	}

	summarize `y' `wt' if `touse', meanonly

	if r(N) == 0 { error 2000 }
	if r(N) == 1 { error 2001 }

	if r(min) < 0 {
		di in red "`y' must be greater than or equal to zero"
		exit 459
	}
	if r(min) == r(max) & r(min) == 0 {
		di in red "`y' is zero for all observations"
		exit 498
	}

	tempname mean
	scalar `mean' = r(mean)

/* Check whether `y' is integer-valued. */

	if "`display'"=="" {
		capture assert `y' == int(`y') if `touse'
		if _rc {
			di in blu "note: you are responsible for " /*
			*/ "interpretation of non-count dep. variable"
		}
	}

/* Print out aweight message. */

	if "`weight'"=="aweight" {
		di in blu "note: you are responsible for interpretation " /*
		*/ "of analytic weights"
	}

/* Remove collinearity. */

	_rmcoll `xvars' [`weight'`exp'] if `touse', `constan'
	local xvars `r(varlist)'

	if "`weight'"=="aweight" {
		_rmcoll `lnalpha' [iw`exp'] if `touse'
			/* aweights produce "sum of wgt" message,
			   which should not be repeated
			*/
	}
	else	_rmcoll `lnalpha' [`weight'`exp'] if `touse'
	local lnalpha `r(varlist)'

/* Run comparison Poisson model. */

	if "`log'"!="" { local nohead "*" }

	if "`from'"=="" {
		if "`lrtest'"!="" {
				    /* we do not do comparison LR test,
		                       but we get starting values from
				       poisson, iterate(0)
				    */
			if "`weight'"!="" {
				local wt `"[iw`exp']"'
			}
			quietly poisson `y' `xvars' `wt' /*
			*/ if `touse', iterate(0) nodisplay /*
			*/ `offopt' `constan'
		}
		else {
			if "`robust'`cns'`cluster'"!="" | /*
			*/ "`weight'"=="pweight" {
				local lrtest "nolrtest"
			}
			`nohead' di in gr _n "Fitting Poisson model:"

			poisson `y' `xvars' [`weight'`exp'] if `touse', /*
			*/ nodisplay `offopt' `constan' `log' `mlopts'
		}

		tempname bp
		mat `bp' = e(b)

		if "`lrtest'"=="" {
			tempname llp
			scalar `llp' = e(ll)
			local kp `e(k)'
		}
	}

/* Fit constant-only model. */

	tempname lna

	if "`constan'"=="" & "`from'"=="" {

	/* Get starting values for constant-only model. */

		if "`offset'"!="" {
			SolveC `y' `offset' [`weight'`exp'] if `touse', /*
			*/ mean(`mean')
			local c = r(_cons)
		}
		else	local c = ln(`mean')

		tempname b0 ll0
		if `c'< . {
			mat `b0' = (`c', 0)
		}
		else	mat `b0' = (0, 0)

		mat colnames `b0' = `y':_cons lnalpha:_cons

		if "`weight'"=="aweight" {
			local wt `"[aw`exp']"'
		}
		else if "`weight'"!="" {
			local wt `"[iw`exp']"'
		}

		`nohead' di in gr _n "Fitting constant-only model:"

		ml model lf gnbre_lf (`y': `y'=, `constan' `offopt') /*
		*/ (lnalpha: `lnalpha') if `touse' `wt', /*
		*/ collinear missing max nooutput nopreserve wald(0) /*
		*/ init(`b0') search(off) `mlopts' `log' `robust'

		local lf0 "lf0(`e(k)' `e(ll)')"
		mat `b0' = e(b)
		mat `lna' = `b0'[1,"lnalpha:"]
		scalar `ll0' = e(ll)
	}
	else {
		mat `lna' = (0)
		mat colnames `lna' = lnalpha:_cons
	}

	if "`from'"=="" {
		mat `bp' = `bp' , `lna'

		if "`constan'"=="" {
			LLnbreg `y' `bp' [`weight'`exp'] if `touse', `offopt'
			if `ll0' > r(lnf) | r(lnf)>=. {
				local initopt "init(`b0')"
			}
		}
		if "`initopt'"=="" {
			local initopt "init(`bp')"
		}

		if "`constan'"=="" | "`lrtest'"=="" {
			`nohead' di in green _n "Fitting full model:"
		}
	}
	else    local initopt `"init(`from')"'

	ml model lf gnbre_lf (`y': `y'=`xvars', `constan' `offopt') /*
	*/ (lnalpha: `lnalpha') if `touse' [`weight'`exp'], collinear /*
	*/ missing max nooutput nopreserve `initopt' search(off) `mlopts' /*
	*/ `log' `scopt' `robust' `clopt' `lf0' /*
	*/ title(Generalized negative binomial regression)

	est local cmd

        if "`score'"!="" {
		label var `scvar1' "Score index for x*b from gnbreg"
		rename `scvar1' `scname1'
		label var `scvar2' "Score index for lnalpha from gnbreg"
		rename `scvar2' `scname2'
		est local scorevars `score'
	}

	if "`llp'"!="" {
		est local chi2_ct "LR"
		est scalar ll_c = `llp'
		est scalar k_c = `kp'
		if e(ll) < e(ll_c) & reldif(e(ll),e(ll_c)) < 1e-5 {
			est scalar chi2_c = 0
				/* otherwise, let it be negative when
				   it does not converge
				*/
		}
		else	est scalar chi2_c = 2*(e(ll)-e(ll_c))
	}

	if "`cluster'"=="" & "`weight'"!="pweight" {
		est scalar r2_p = 1 - e(ll)/e(ll_0)
	}

	est local offset1 "`offvar'"
	est local predict "gnbreg_p"
        est local cmd     "gnbreg"


	Display, `irr' level(`level')
	error `e(rc)'
end

program define Display
	syntax [, Level(cilevel) IRr]
	if "`irr'"!="" {
		version 9: ///
		ml di, level(`level') eform(IRR) first plus nofootnote

		di in smcl in gr /*
		*/ "     lnalpha {c |}  (type gnbreg to see ln(alpha) " /*
		*/ "coefficient estimates)" _n  /*
		*/ "{hline 13}{c BT}{hline 64}"
		ml_footnote
	}
	else	version 9: ml di, level(`level') nofootnote

	DispLr  /* Note: This is not a chi-square, chibar-square, or */
		/* anything known to man.*/
	_prefix_footnote
end

program define DispLr
	if "`e(chi2_ct)'"!="LR" { exit }
	if e(chi2_c) < 1e5 {
		local fmt "%9.2f"
	}
	else 	local fmt "%9.2e"
        di in gr "Likelihood-ratio test of alpha=0:   " /*
        */ in gr "chi2(" in ye e(k)-e(k_c) in gr ") =" in ye `fmt' /*
        */ e(chi2_c) in gr _col(59) "Prob > chi2 = " in ye %6.4f /*
        */ chiprob(e(k)-e(k_c), e(chi2_c))
end

program define SolveC, rclass /* modified from poisson.ado */
	gettoken y  0 : 0
	gettoken xb 0 : 0
	syntax [fw aw pw iw] [if] , Mean(string)
	if "`weight'"=="pweight" | "`weight'"=="iweight" {
		local weight "aweight"
	}
	summarize `xb' `if', meanonly
	if r(max) - r(min) > 2*709 { /* unavoidable exp() over/underflow */
		exit /* r(_cons) >= . */
	}
	if r(max) > 709 | r(min) < -709  {
		tempname shift
		if r(max) > 709 { scalar `shift' =  709 - r(max) }
		else		  scalar `shift' = -709 - r(min)
		local shift "+`shift'"
	}
	tempvar expoff
	qui gen double `expoff' = exp(`xb'`shift') `if'
	summarize `expoff' [`weight'`exp'], meanonly
	return scalar _cons = ln(`mean')-ln(r(mean))`shift'
end

program define LLnbreg, rclass
	gettoken y 0 : 0
	gettoken b 0 : 0
	syntax [fw aw pw iw] [if] [, offset(string)]
	if "`weight'"!="" {
		if "`weight'"=="fweight" {
			local wt `"[fw`exp']"'
		}
		else	local wt `"[aw`exp']"'
	}
	tempvar z lnalpha m
	tempname m b1 b2
	mat `b1' = `b'[1,"`y':"]
	mat `b2' = `b'[1,"lnalpha:"]

	quietly {
		mat score double `z' = `b1' `if'
		summarize `z' `wt' `if', meanonly
		local nobs `r(N)'

		if "`offset'"!="" {
			local zo "(`z'+`offset')"
		}
		else	local zo "`z'"

		mat score double `lnalpha' = `b2' `if'
		local bound -30
		replace `lnalpha' = `bound' if `lnalpha'<`bound' & `lnalpha'< .

		gen double `m' = exp(-`lnalpha')

		replace `z' = lngamma(`m'+`y') - lngamma(`y'+1) /*
		*/ - lngamma(`m') - `m'*ln(1+exp(`zo'+`lnalpha')) /*
		*/ - `y'*ln(1+exp(-`zo'-`lnalpha')) `if'

		summarize `z' `wt' `if', meanonly
		if r(N) != `nobs' { exit }
		if "`weight'"=="aweight" {
			ret scalar lnf = r(N)*r(mean)
				/* weights not normalized in r(sum) */
		}
		else	ret scalar lnf = r(sum)
	}
end