xtgee.ado

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

		tempvar ei t cc lagw

		sort `ivar' `tvar' `obs'
		while `iterate' < `iter' & `diff' > `tol' {

			mat `b0' = `b1'
			scalar `phi' = 0.0

			* switch(link) page 4

			cap drop `mui'
			mat score double `mui' = `b1'
			`addoff' /* replace `mui' = `mui' + `offset' */
			xtgee_plink "`link'" `mui' `binom'

			* switch (var_mean_rel) page 4

			cap drop `ei'
			gen double `ei' = $S_X_depv - `mui'
			xtgee_elink "`family'" `mui' `ei' `binom'

			cap drop `cc'

			gen double `cc' = `ww'*`ei'*`ei'
			cap drop `t'

/* Use Liang Zeger 1986 moment estimators instead of Karim estimators */
			summ `cc'
			if r(N) != `nobs' {
				noi di in red "estimates diverging " /*
					*/ "(missing predictions)"
				exit 430
			}

/* N-p changed to N */
			scalar `phi' = (`sumw'-`nmp'*`p') / r(sum)

			by $S_X_ivar : gen `t' = /*
				*/ .5*`ww'[_N]*_N*(_N-1) if _n==_N
			summ `t'
			tempname cF
/* N-p changed to N */
			scalar `cF' = r(sum) - `nmp'*`p'

			cap drop `lagw'
			local bw1 = `band'+1
			by $S_X_ivar : gen double `lagw' = _N/(_N-_n+1) /*
			*/ if _n <= `bw1'

			_GEEUC `ei' `ww', a(`alpha') ty(`corr') n(`maxni') /*				*/ b(`band') l(`lagw') `topt' by($S_X_ivar) /*
				*/ fixsing

			if "`trace'" != "" {
                                tempname traa aa
                                mat `traa' = `b1'`addme'
                                mat coleq `traa' = beta
                        }

			if "`corr'" == "exchan" {
				scalar `alpha' = `alpha'*`phi'/`cF'

                                if "`trace'" != "" {
                                        local aaa = `alpha'
                                        mat `aa'  = (`aaa')
                                        mat colnames `aa' = alpha
                                        mat `traa' = `traa',`aa'
                                }

				if ((`alpha'<-1) | (`alpha'>1)) {
					noi di in red "estimates diverging " /*
					*/ "(absolute correlation > 1)"
					exit 430

				}
			}
			else if "`corr'" != "fixed" & "`corr'" != "indep" {
				if "`corr'"=="statm" | "`corr'"=="arm" {
					local kk = colsof(`alpha')
					tempname ttt
					local i 1

/*                        noi mat list `alpha'
			noi di in red `phi' */

					while `i' <= `kk' {
						scalar `ttt' = `phi' / /*
					*/ (`sumw'+`nclust'-`nclust'*`i'-/*
							*/ `nmp'*`p') /*
						p out with nmp option
							*/
						if `alpha'[1,`i'] != 0 {
						mat `alpha'[1,`i'] = /*
					*/ `alpha'[1,`i'] * `ttt'
						}
						local i = `i'+1
					}

/*                        noi mat list `alpha'
			noi di in red `ttt' */
                                        if "`trace'" != "" {
                                                mat `aa' = `alpha'[1,1..`bw1']
                                                mat colnames `aa' = `aanames'
                                                mat coleq `aa' = alpha
                                                mat `traa' = `traa',`aa'
                                        }
				}
				else if "`corr'" == "nonst" {
					tempname ttt
					scalar `ttt' = `phi' / (`nclust')
					mat `alpha' = `ttt' * `alpha'
				}
				else {
					/* unstructured */
					tempname ttt
					scalar `ttt' = `phi' / (`nclust')
					mat `alpha' = `ttt' * `alpha'
				}
			}

                        if "`trace'" != "" {
                                mat rownames `traa' = " "
                                noi mat list `traa', noheader
                        }

			* switch (corstruct) page 6

			if "`corr'" == "arm" {
				rarm `R' `maxni' `alpha' `bw1'
			}
			else _GEERC , c(`corr') r(`R') n(`maxni') /*
			*/ a(`alpha') b(`bw1')

			* Bottom of page 6

			cap drop `mui'
			mat score double `mui' = `b1'
			`addoff' /* replace `mui' = `mui' + `offset' */

			_GEEBT $S_X_idep, mui(`mui') N(`binom') by(`ivar') /*
			*/ `warg' xy(`S1') xx(`S2') R(`R') /*
			*/ b1(`b1') depv($S_X_depv) n(`maxni') `topt'

			* get S1 and S2 page 8
			mat `b1' = syminv(`S2')
			mat `b1' = `b1'*`S1'
			mat `b1' = `b1''

			mat colnames `b1' = `sna'

			* now get max diff between oldbeta and newbeta
			local diff = mreldif(`b1', `b0')
			`prel' in gr "Iteration " `iterate' /*
			*/ ": tolerance = " in ye `diff'
			local iterate = `iterate'+1
		}
		if `iterate' >= `iter' {
			local finalrc 430
		}
		else	local finalrc 0

		* Page 9 loop
		cap drop `zi'

		scalar `phi' = 0.0
		cap drop `mui'
		mat score double `mui' = `b1'
		`addoff' /* replace `mui' = `mui' + `offset' */

		if "`score'" != "" {
			qui gen double `score' = .
			local sarg "score(`score')"
		}

		/* Get final variance estimates */

		if "`weight'" == "fweight" {
			local fwt "fwt"
		}
		_GEEBT $S_X_idep, mui(`mui') N(`binom') by(`ivar') /*
		*/ `warg' phi(`phi') `sarg' `robust' xy(`S1') xx(`S2') /*
		*/ R(`R') b1(`b1') depv($S_X_depv) n(`maxni') `topt' `fwt'

		scalar `phi' = `avgni'*`phi' / (`sumw' -`nmp'*`p') /*  N-p */

		global S_E_dof = `sumw'-`nmp'*`p' /* N-p */
		local mydof $S_E_dof
		global S_E_scp "$S_X_scp"

		local wtex  "[aweight=`ww']"
		FixScale $S_X_depv `mui' `binom' "`wtex'"
		if "$S_X_fix" != "" {
			scalar `phi' = $S_1
		}

		mat `S2' = syminv(`S2')

		tempname nse rse

		mat `nse' = `phi'*`S2'
		mat `rse' = `S2'*`S1'
		mat `rse' = `rse'*`S2'
		mat drop `S1'
		mat drop `S2'
	}
	mat colnames `nse' = `sna'
	mat rownames `nse' = `sna'
	mat colnames `rse' = `sna'
	mat rownames `rse' = `sna'

	if "`weight'" == "fweight" {
		local nobs = `sumw'
	}
	else if "`weight'" == "aweight" {
		local nclust `keep'
	}

	if "`robust'" != "" & "`weight'" !="fweight" {
		local nclust `keep'
	}

	global S_X_9 = "$S_E_dis"
	global S_X_A = "$S_E_scp"
	if "`robust'" != "" {
		local df = `nclust'-1
		local factor = (`nclust')/(`df')
		/* The usual adjustment for regression includes this as well*/
		/* Making this adjustment prevents Robust Var. estimate from
		   being invariant to scale of weights. Only make it when
		   asked via RGF option */
		if "`rgf'" != "" {
 			if "`family'"=="gauss" {
 				local factor =`factor'*(`sumw'-1)/(`sumw'-`p')
 			}
		}
		mat `rse' = `rse' * `factor'
		mat `rse' = (`rse' + `rse'')/2
		est post `b1' `rse', depname(`depname') obs(`nobs')
		global S_X_vce "Semi-robust"
	}
	else {
		est post `b1' `nse', depname(`depname') obs(`nobs')
	}

	if "`corr'" == "arm" {
		rarm `R' `maxni' `alpha' `bw1'
	}
	else _GEERC , c(`corr') r(`R') n(`maxni') a(`alpha') b(`bw1')
	if _caller() < 6 {
		mat S_E_R = `R'
	}
	est matrix R `R'

	capture qui test $S_X_HAS, min
	est local chi2type "Wald"
	if _rc == 0 {
		est scalar chi2 = r(chi2)
		est scalar df_m    = r(df)
	}
	else {
		est scalar chi2 = .
		est scalar df_m = 0
	}
	global S_E_chi2 = e(chi2)	/* double saves */
	global S_E_chdf = e(df_m)	/* double saves */
	global S_E_prob = chiprob(e(df_m),e(chi2))	/* double saves */

	if "$S_LOG" != "" & "`display'" == "" {
		noi di
	}

	/* Save results */

	SaveMac `mydof' `nobs' `nclust' `rmaxni' `minni' `avgni' `tol' `diff' /*
		*/ `phi' `family' `link' `corr' "`tvar'" "`band'" `level'
	est local ivar "`rivar'"
	if "`nmp'" == "1" {
		est local nmp "nmp"
	}
	est local rgf "`rgf'"
	global S_E_ivar "`rivar'"  /* double saves */
	fixwords
	if "`robust'" != "" {
		global S_E_cvn "`rivar'"
	}

	est local offset "`offstr'"
	est local denom "`oarg'"

	est local estat_cmd "xtgee_estat"
	est local predict "xtgee_p"
	est local cmd "xtgee"
	global S_E_cmd "`e(cmd)'"		/* double saves */

	global S_LOG

	if "`score'" != "" {
		qui replace `score' = `score'/sqrt(`phi')
		keep `merge' `score'
		sort `merge'
		tempfile tmpf
		qui save `"`tmpf'"'
		restore
		cap confirm variable _merge
		if _rc==0 {
			tempvar _merge1
			rename _merge `_merge1'
		}
		sort `merge'
		cap merge `merge' using `"`tmpf'"'
		if _rc {
			if _rc==1 { exit 1 }
			noi di in red "Unable to generate score `score'"
			exit 198
		}
		drop _merge
		cap rename `_merge1' _merge
		est local scorevars `score'
	}
	else {
		restore
	}
	est scalar rc = `finalrc'
	mat repost, esample(`touse')

end

program define geerpt
	global S_X_maxn
	global S_X_corr
	global S_X_mvar
	global S_X_link
	global S_X_idep
	global S_X_depv
	global S_X_con
	global S_X_9
	syntax [, LEvel(real $X_level) EForm]

	if `"`e(cmd)'"'!="xtgee" { error 301 }

	global S_E_ef
	if "`eform'" != "" {
		if substr("`e(family)'",1,3) == "neg" {
			local eform "eform(IRR)"
		}
		else if "`e(family)'"=="Poisson" & "`e(link)'" != "probit" & /*
		*/ "`e(link)'" != "logit" {
			local eform "eform(IRR)"
		}
		else if "`e(link)'" == "probit" & "`e(family)'" != "Poisson" {
			local eform "eform(ExpB)"
		}
		else if "`e(link)'" == "logit" & "`e(family)'" != "Poisson" {
			local eform "eform(Odds Ratio)"
		}
		else if "`e(link)'" == "opower" {
			local eform "eform(Odds Ratio)"
		}
		else {
			local eform "eform(e^coef)"
		}
		global S_E_ef "`eform'"		/* double save */
	}

	local q "_col(68)"
	local l = 71 + 8 - length("`e(ivar)'")
	noi di in green "GEE population-averaged model"  /*
	*/ in gr _col(49) "Number of obs" `q' "=" in ye _col(70) %9.0g e(N)

	local ivar = abbrev("`e(ivar)'", 12)
	if "`e(tvar)'"!="" {
		local tvar = abbrev("`e(tvar)'", 12)
		local skip = max(43-length("`ivar' `tvar'"),1)

		di in gr "Group and time vars:" in ye _col(`skip') /*
		*/ "`ivar' `tvar'" /*
		*/ in gr _col(49) "Number of groups" `q' "=" /*
		*/ in ye _col(70) %9.0g e(N_g)
	}
	else {
		local skip = max(43-length("`ivar'"),1)
		di in gr "Group variable:" in ye _col(`skip')  /*
		*/ "`ivar'" /*
		*/ in gr _col(49) "Number of groups" `q' "=" /*
		*/ in ye _col(70) %9.0g e(N_g)
	}

	local skip = max(43 - length("`e(link)'"),1)
	noi di in green "Link:" in ye _col(`skip') "`e(link)'" /*
	*/ in gr _col(49) "Obs per group: min" `q' "=" in ye _col(70) /*
	*/ %9.0g e(g_min)

	local skip = max(43 - length("`e(family)'"),1)
	noi di in green "Family:" in ye _col(`skip') "`e(family)'" /*
	*/ in gr _col(64) "avg" `q' "=" in ye _col(70) /*
	*/ %9.1f e(g_avg)

	local skip = max(43 - length("`e(corr)'"),1)
	noi di in green "Correlation:" in ye _col(`skip') /*
	*/ "`e(corr)'" in gr _col(64) "max" `q' "=" /*
	*/ in ye _col(70) %9.0g e(g_max)

	if e(chi2)>999999 {
		local cfmt "%9.0g"
	}
	else	local cfmt "%9.2f"

	if !missing(e(df_r)) {
		local model as txt _col(49) "F(" ///
			as res %4.0f e(df_m) as txt "," ///
			as res %7.0f e(df_r) as txt ")" ///
			`q' "=" _col(70) as res %9.2f abs(e(F))
		local pvalue _col(49) as txt "Prob > F" `q' "=" ///
			as res _col(73) %6.4f Ftail(e(df_m),e(df_r),e(F))
	}
	else {
		local model as txt _col(49) "Wald chi2(" ///
			as res e(df_m) as txt ")" ///
			`q' "=" _col(70) as res `cfmt' e(chi2)
		local pvalue _col(49) as txt "Prob > chi2" `q' "=" ///
			as res _col(73) %6.4f chiprob(e(df_m),abs(e(chi2)))
	}

	noi di `model'
	noi di in green /*
	*/ "Scale parameter: "  _col(34) in ye %9.0g e(phi) `pvalue'

	if "`e(deviance)'" != "" & "`e(corr)'" == "independent" {
		if e(deviance)>999999 {
			local cfmt "%10.0g"
		}
		else	local cfmt "%10.2f"
		noi di in gr _n "Pearson chi2(" in ye e(df_pear) in gr "):" /*
		*/ in ye _col(33) `cfmt' e(chi2_dev) /*
		*/ in gr _col(49) "Deviance" `q' "=" in ye _col(69) /*
		*/ `cfmt' e(deviance)

		noi di in gr "Dispersion (Pearson):" in ye _col(34) /*
		*/ %9.0g e(chi2_dis) /*
		*/ in gr _col(49) "Dispersion" `q' "=" in ye _col(70) /*
		*/ %9.0g e(dispers)
	}

	di
	est di, level(`level') `eform'
	if "`e(disp)'" == "`e(phi)'" {
		cap confirm number `e(scale)'
		local j = _rc
		if "`e(scale)'" == "dev" {
			local mm "deviance based dispersion"
		}
		else if "`e(scale)'" == "x2" {
			local mm "Pearson X2-based dispersion"
		}
		else {
			local mm "`e(scale)'"
		}
		if /*
*/ (("`e(family)'"=="Gaussian" | "`e(family)'"=="gamma") & /*
*/ "`e(scale)'"!="x2") | /*
*/ (("`e(family)'"=="Poisson" | "`e(family)'"=="binomial") & /*
*/ "`e(scale)'"!="1") {
			noi di in gr  /*
			*/ "(Standard errors scaled using square root of `mm')"
		}
	}
	_prefix_footnote
	if e(rc) & !missing(e(rc)) {
		error e(rc)
	}
end

program define rarm
	args R n a b

	mat `a'[1,1] = 1
	tempname tmp wt al
	local b2 = `b'-1
	mat `tmp' = J(1,`n',0)
	local i 1
	while `i' <= `n' {
		mat `tmp'[1,`i'] = `a'[1,`i']
		local i = `i'+1
	}
	mat `wt' = J(`b2',`b2',0)
	local i 1
	while `i' <= `b2' {
		local j = `i'
		while `j' <= `b2' {
			local k = `j'-`i'+1