*! version 5.4.4  03oct2002
program aglm, eclass byable(onecall)
	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 _by() { 
		local by "by `_byvars'`_byrc0':"
	}
	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 `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_running SGLM_nonstan
                        exit 459
                }
                else {
                        mac drop SGLM_running SGLM_nonstan
                        exit _rc
                }
        }
end

program glm_7, eclass byable(recall) sort

/* Parse. */

	syntax varlist(numeric) [fw aw pw iw] [if] [in] [, EForm        /*
	*/ FROM(string) LEvel(integer $S_level) OFFset(varname numeric) /*
	*/ LNOFFset(varname numeric) noCONstant Robust CLuster(varname) /*
	*/ Family(string) Link(string) SEARCH IRLS LTOLerance(real 1e-6) /*
	*/ SCore(string) noLOg OIM noDISPLAY FISHER(integer -1) noDOTS    /*
	*/ ITERate(integer 50) TRAce SCAle(string) MU(varname)        /*
	*/ NWEST(string) BSTRAP BREP(int -1) JKNIFE JKNIFE1 UNBiased    /*
	*/ T(string) noHEADer BHHH OPG DISP(real 1) VFactor(real 1) * ]

	if "`bhhh'" != "" {
		local opg opg
	}

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

	local ltol `"`ltolerance'"'
	
	if `ltol'>=1 | `ltol'<0 {
		noi di as err "ltolerance() must be in [0,1)"
		exit 198
	}

	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'

/* Check syntax. */

	if `level' < 10 | `level' > 99 {
		di as err "level() must be between 10 and 99"
		exit 198
	}


	if `"`score'"'!="" {
		if "`irls'"!="" {
			di as err "deviance scores are calculated " /*
			*/ "post-estimation 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 `brep' ~= -1 {
			di as err "brep() only valid with bstrap"
			exit 198
		}
	}

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

	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'`nwest'`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
	}

	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 "`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 `"`nwest'"' == "" {
		if `"`t'"' ~= "" {
			di as err "t() only valid with nwest()"
			exit 198
		}
	}

	if "`nwest'" != "" {	/* 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'" != "" {
		if "`nwest'`jknife1'" != "" {
			di as err "cannot specify cluster() with " /*
				*/ "this variance estimate"
			exit 198
		}
	}

/* Process offset/exposure. */

	if "`lnoffset'"!="" {
		capture assert `lnoffset' > 0 if `touse'
		if _rc {
			di as err "lnoffset() 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

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

	MapNW `nobs' `nwest'
	local nwnam `"`r(nwnam)'"'
	local nwlag `"`r(nwlag)'"'

	local fwt "1"
	local awt "1"
	tempvar wt
	if "`weight'"!="" {
		qui gen double `wt' `exp' if `touse'
		local awt "`wt'"
		if `"`weight'"'=="aweight" {
			qui summ `wt' 
			qui replace `wt' = `wt'/r(mean)
		}
		else if `"`weight'"'=="fweight" {
			if "`bstrap'`nwest'" != "" {
				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
				}
			}
			qui summ `touse' [fw=`wt'] if `touse' 
			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 "`nwnam'`opg'`jknife'`jknife1'" != "" {
			di as err "robust only allowed with hessian matrices"
			exit 198
		}
	}

	if "`weight'" != "" {
		if "`weight'" != "fweight" {
			if "`nwnam'`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
		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']                        
				scalar `Wscale' = r(mean)
				cap 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' 
				$SGLM_mu `mu'
				mat `b0' = e(b)
				local k = colsof(`b0')
				local from "`b0'"
				replace `z' = sum(`wt'*(`y'-`mu')^2/`v')
				global SGLM_ph = `z'[_N]/(`nobs'-`k')
			}
		}

		if "`from'"!="" & "`search'" == "" {
			local initopt `"init(`from') search(off)"'
		}
		else	local initopt `"search(on) maxfeas(100)"'
		if `reg' & `"`offopt'"'=="" { 
			local iterate 0 
			local nowarn "nowarning" 
		}
		local initopt "`initopt' iter(`iterate') `nowarn' "

		if `"`scale'"'=="" {    /* default */
			local scale `scale1'
			local cd
		}
		else if `scale'==0 {    /* Pearson X2 scaling */
			if `scale1' {
				local cd "square root of Pearson"
				local cd "`cd' X2-based dispersion"
			}
		}
		else {                  /* user's scale parameter */
			if !`scale1' | (`scale1' & `scale'!=1) {
				local cd "dispersion equal to square"
				local cd "`cd' root of `scale'"
			}
		}

		if `scale1' == 1 & `scale' != `scale1' { 
			global SGLM_s1 = `scale1' 
			local fixme 1
		}
		else {
			global SGLM_s1 = `scale' 
			local fixme 0 
		}

		/* 
			        Need to handle k==N case separately 
				This must be handled by ml because the 
				variance matrix is missing. The irls 
				works fine in this case.  Maybe the right
				thing to do is to just pass it on to the
				irls method since the variance matrix is 
				missing anyway.
		*/
		tempname newdev

		if "`opg'`unbiased'" != "" {
			if "`scopt'" == "" {
				tempvar scvar
				local scopt score(`scvar')
			}
			if "`unbiased'" != "" {
				local cll "`clopt'"
				local robust
				local clopt
			}
		}

		if (!`iterate') { `iflog' di }
 
		ml model d2 glim_lf /* 
		*/ (`y': `y' = `xvars', `constant' `offopt') /*
		*/ [`weight'`exp'] if `touse', collinear missing max nooutput /*
		*/ nopreserve `mlopts' title(Generalized Linear Models) /*
		*/ `scopt' `robust' `clopt' `log' `initopt' `trace' 
		local rc = cond(e(rc) == 430 & `reg', 0, e(rc))

		if `fixme' {
			tempname V fix
			mat `V' = e(V)
			if `scale' { scalar `fix' = `scale' }
			else         scalar `fix' = $SGLM_ph 
			mat `V' = `fix'*`V'
			estimates repost V = `V'
		}
			

		local setype "OIM"

		if "`unbiased'" != "" {
			tempvar hat eta mu v dmu dv d2mu z W
			tempname  ll Wscale

			qui _predict double `eta', xb   /* sic -- no nooffset */
			qui $SGLM_L 1 `eta' `mu'
			qui $SGLM_V 1 `eta' `mu' `v'
			qui $SGLM_L 2 `eta' `mu' `dmu'

			qui gen double `z' = `eta' + /*
				*/ (`y'-`mu')/`dmu' `moffset' if `touse'

			qui $SGLM_V 2 `eta' `mu' `dv'
			qui $SGLM_L 3 `eta' `mu' `d2mu'
			qui gen double `W' = `dmu'*`dmu'/`v' - /*
				*/ (`mu'-`y')*(`dmu'*`dmu'*`dv'/(`v'*`v') - /*
				*/ `d2mu'/`v') if `touse'
			nobreak {
				estimates hold `ll'
				qui reg `z' `xvars' [iw=`W'*`wt'], /*
					*/ mse1 `constant'
				qui _predict double `hat' if `touse', /*
					*/ stdp `offset'
				qui replace `hat' = `hat'*`hat'*`v'
				estimates unhold `ll'
			}

			qui replace `scvar' = `scvar' / sqrt(1-`hat')
			qui _robust `scvar' [fweight=`fwt'] if `touse', /*
				*/ `cll' minus(0)
			local setype "Sandwich" 
			if "`cll'" != "" {
				local setype "Modified `setype'"
			}
			local setype "Unbiased `setype'"
			local vce    "Robust"
		}

		if "`opg'" != "" {
			tempname V1
			mat `V1' = e(V)
			local nms : colnames `V1'
			local ncol = colsof(`V1')
			mat `V1' = I(`ncol')
			mat colnames `V1' = `nms'
			mat rownames `V1' = `nms'
			_robust `scvar' [fweight=`fwt'] if `touse', var(`V1')
			mat `V1' = syminv(`V1')
			estimates repost V = `V1'
			local setype "Outer product of the gradient"
			local vce    "OPG"
		}

		if "`robust'" != "" { 
			local setype "Sandwich" 
			if "`clopt'" != "" {
				local setype "Modified `setype'"
			}
			if "`unbiased'" != "" {
				local setype "Unbiased `setype'"
			}
		}

		if "`nwnam'" != "" {
			tempname V
			mat `V' = e(V)
			local p = colsof(`V')-1
			local vce "Newey-West"
			capture drop `eta'
			capture drop `mu'
			capture drop `v'
			capture drop `dmu'
			capture drop `z'
			qui _predict `eta', xb
			qui $SGLM_L 1 `eta' `mu' 
			qui $SGLM_V 1 `eta' `mu' `v'
			qui $SGLM_L 2 `eta' `mu' `dmu'
			qui gen double `z' = `dmu'*(`y'-`mu')/`v'
			sort `touse' `tvar'
			qui NWest "`nwnam'" "`nwlag'" "`wt'" "`z'" /*
				*/ "`touse'" `"`constant'"' `"`xvars'"'
			tempname vm
			mat `vm' = r(V)
			local pm1 = `p'+1
			`nwnam' `nwlag' 0 `nobs' `pm1'  
			local delta = 1
			if `scale1' == 0 { local delta = $SGLM_ph } 
			tempname ss
			scalar `ss' = 1/(`delta'*`delta')
			mat `V' = `ss'*`V'*`vm'*`V''
			local setype "Weighted sandwich (`r(setype)')"
			local sewtype "`r(sewtype)'"
			estimates repost V = `V'
		}
		if "`jknife1'" != "" {
			tempname b V
			mat `b' = e(b)
			mat `V' = e(V)
			qui JK1nr `b' `V' `y' "`xvars'" `eta' `mu' `v' `dmu' /*
				*/ "`tvar'" `z' "`moffset'" "`constant'" /*
				*/ `wt' `touse' "`weight'"
			local vce    "1-step JKnife"
			local setype "One-step Jackknife"
			local p = colsof(`b')
			tempname ss
			scalar `ss' = (`nobs'-`p')/`nobs'
			mat `V' = `ss'*`V'
			estimates repost V = `V'
		}
		if "`jknife'" != "" {
			tempname b V
			mat `b' = e(b)
			mat `V' = e(V)
			noi JKnife `b' `V' `y' "`xvars'" "`tvar'" /*
				*/ "`offopt'" "`constant'" "`argfam'" "" /*
				*/ "`arglink'" "[`weight'`exp']" `touse' /*
				*/ "`cluster'" `dots' "`weight'" "`wt'"
			local mm = `nobs'
			local vce    "Jackknife"
			local setype "Jackknife"
			if "`cluster'" != "" { local mm = $SGLM_nc }
			else local mm = 1
			local p = colsof(`b')
			local nnn = `nobs'-$SGLM_bm
			tempname ss
			if "`cluster'" != "" {
				scalar `ss' = (`nnn'-`p')/`nnn'
			}
			else {
				scalar `ss' = (`nnn'-`p')/`nnn'
			}
			mat `V' = `ss'*`V'
			estimates repost V = `V'
		}
		if "`bstrap'" != "" {
			tempname b V
			mat `b' = e(b)
			mat `V' = e(V)
			noi Bstrap `b' `V' `y' "`xvars'" /*
				*/ "`offopt'" "`constant'" "`argfam'" "" /*
				*/ "`arglink'" "[`weight'`exp']" `touse' /*
				*/ `brep' "`clopt'" `dots'
			local vce    "Bootstrap"
			local setype "Bootstrap"
			local p = colsof(`b')
			tempname ss
			scalar `ss' = (`nobs'-`p')/(`nobs'*(`brep'-$SGLM_bm))
			mat `V' = `ss'*`V'
			estimates repost V = `V'
		}
		if `vfactor' != 1 {
			tempname V
			mat `V' = `vfactor'*e(V)
			estimates repost V = `V'
		}

		capture drop `eta'
		capture drop `mu'
		capture drop `v'
		qui predict double `eta', xb
		qui $SGLM_L 1 `eta' `mu' 
		qui $SGLM_mu `mu'
		qui $SGLM_V 1 `eta' `mu' `v'
		mat `b0' = e(b)
		local dfm = colsof(`b0')
		capture drop `z'
		qui gen double `z' = sum(`wt'*(`y'-`mu')^2/`v')
		local chi2  = `z'[_N]
		local df = e(N)-`dfm'

		if $SGLM_s1 { global SGLM_ph = $SGLM_s1 }

		capture drop `z'
		qui $SGLM_V 3 `eta' `mu' `z'
		qui replace `z' = sum(`wt'*`z')
		scalar `newdev' = `z'[_N]
		local disp  = $SGLM_ph
		est scalar vf = `vfactor'
		est scalar rc = `rc'
		est scalar aic = (-2*e(ll)+2*`dfm')/`nobs'

	}
	/* Classical IRLS */
	else {
		quietly {
			tempvar  eta mu dmu v W z
			tempname Wscale b0
			if "`offset'" != "" {
				local moffset = "-`offset'"
			}
			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)
				}
			}

			if `reg' { local iterate 1 }

			capture drop `eta'
			$SGLM_L 0 `eta' `mu'
			tempvar dev
			tempname newdev oldev 
			scalar `oldev'  =  0
			scalar `newdev' =  1
			local i 1
			`iflog' di
			while abs(`newdev'-`oldev')>`ltol' & `i'<=`iterate' {
				capture drop `v'
				capture drop `dmu'
				capture drop `W'
				capture drop `z'
				scalar `oldev' = `newdev'
				$SGLM_V 1 `eta' `mu' `v'
				$SGLM_L 2 `eta' `mu' `dmu'
				gen double `z' = `eta' + /*
					*/ (`y'-`mu')/`dmu' `moffset' if `touse'
				gen double `W' = `disp'*`dmu'*`dmu'/`v' /*
					*/ if `touse'
				summ `W' [aw=`wt']                        
				scalar `Wscale' = r(mean)

				cap areg `z' `xvars' /*
				*/   [aw=`W'*`wt'/`Wscale'], absorb(indiv) `constant'
				if e(N)>=. {exit 459}
				if _rc {exit _rc}
				capture drop `eta'
				capture drop `mu'
				capture drop `dev'
				predict double `eta' if `touse', xbd
				if "`offset'" != "" {
					replace `eta' = `eta'+`offset'
				}
				$SGLM_L 1 `eta' `mu'
				$SGLM_mu `mu'
				$SGLM_V 3 `eta' `mu' `dev'
				replace `dev' = sum(`wt'*`dev')
				scalar `newdev' = `dev'[_N]/`disp'
				`iflog' di as txt /*
					*/ `"Iteration `i' : deviance = "' /*
					 */ as res %9.0g `newdev'
				if "`trace'" != "" {
					tempname beta
					`iflog' mat list e(b), noblank noheader
					`iflog' di as txt "{hline 78}" 
					`iflog' di
				}
				local i = `i'+1
			}
			if "`unbiased'`jknife1'`jknife'" != "" {
				cap regress `z' `xvars' [iw=`W'*`wt'], /*
					*/ mse1 `constant'
				if e(N)>=. {exit 459}
				if _rc {exit _rc}
				tempvar h
				qui _predict double `h' if `touse', stdp 
				qui replace `h' = `h'*`h'*`v'
				local minus "minus(0)"
			}
			else    local h = 0
			local rc = cond(abs(`newdev'-`oldev')>`ltol' /* 
				*/ & !`reg',430,0)
			replace `z' = sum(`wt'*(`y'-`mu')^2/`v')
			local chi2  = `z'[_N]/`disp'
			local p     = e(df_m)
			local df    = `nobs'-`p'-(`"`constant'"'=="")
			local dispc = `chi2'/`df'
			local dispd = `newdev'/`df'

			global SGLM_ph = `dispc'

			/* Apply appropriate scaling */

			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"
					local cd "`cd' X2-based dispersion"
				}
			}
			else if `scale'==-1 {   /* deviance scaling */
				local delta `dispd'
				local cd "square root of deviance-based"
				local cd "`cd' dispersion"
			}
			else {                  /* user's scale parameter */
				local delta `scale'
				if !`scale1' | (`scale1' & `scale'!=1) {
					local cd "dispersion equal to square"
					local cd "`cd' 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'
				}
			}

			tempname b V df_a
			mat `b' = e(b)
			mat `V' = e(V)/(e(rss)/e(df_r))
			scalar `df_a'=e(df_a)

			local setype "EIM"
			local vce "EIM"
			if "`robust'`oim'`nwnam'`opg'" != ""  {
				capture drop `W'
				capture drop `z'
				capture drop `v'
				capture drop `dmu'

				tempvar xb
				_predict double `xb', xb /* sic--no nooffset */
				$SGLM_V 1 `eta' `mu' `v'
				$SGLM_L 2 `eta' `mu' `dmu'

				tempvar ys yi
				gen double `ys' = `dmu'*(`mu'-`y')/`v'
				gen double `yi' = `dmu'*`dmu'
				if "`oim'" != "" {
					local setype "OIM"
					tempvar dv d2m
					$SGLM_V 2 `eta' `mu' `dv'
					$SGLM_L 3 `eta' `mu' `d2m'
					replace `yi' = `yi' + /*
						*/ `d2m'*(`mu'-`y') - /*
						*/ `dmu'*`dv'*`ys'
				}
				replace `yi' = `awt' * `yi' / `v'
				replace `ys' = `awt' * `ys' / sqrt(1-`h')

				reg `xb' `xvars' if `touse' /*
					*/ [iweight=`fwt'*`yi'], `constant' mse1
				mat `V' = e(V)
				local vce "OIM"
				if "`robust'" != "" {
					if "`oim'" != "" {
						local vce "Robust"
					}
					else {
						local vce "Semi-Robust"
					}
					_robust `ys' [fweight=`fwt'] /*
						*/ if `touse', var(`V') /*
						*/ `clopt' `minus'
					local setype "Sandwich"
					if "`clopt'" != "" {
						local setype "Modified `setype'"
					}
					if "`unbiased'" != "" {
						local setype "Unbiased `setype'"
					}
				}
				else if "`opg'" != "" {
					local nms : colnames `V'
					local ncol = colsof(`V')
					mat `V' = I(`ncol')
					mat colnames `V' = `nms'
					mat rownames `V' = `nms'
					replace `ys' = `ys'/`disp'
					_robust `ys' [fweight=`fwt'] /*
						*/ if `touse', var(`V')
					mat `V'=`delta'*`delta'*syminv(`V')
					local setype /*
					*/ "Outer product of the gradient"
					local vce    "OPG"
				}
				else if "`nwnam'" == "" {
					tempname ss
					scalar `ss' = `delta'*`Wscale'
					mat `V' = `ss'*`V'
				}
			}

			if "`nwnam'" != "" {
				local vce "Newey-West"
				capture drop `v'
				capture drop `dmu'
				capture drop `z'
				$SGLM_V 1 `eta' `mu' `v'
				$SGLM_L 2 `eta' `mu' `dmu'
				gen double `z' = `dmu'*(`y'-`mu')/`v'
				sort `touse' `tvar'
				qui NWest "`nwnam'" "`nwlag'" "`wt'" "`z'" /*
					*/ "`touse'" `"`constant'"' `"`xvars'"'
				tempname vm
				mat `vm' = r(V)
				local pm1 = `p'+1
				`nwnam' `nwlag' 0 `nobs' `pm1'  
				local setype "Weighted sandwich (`r(setype)')"
				local sewtype "`r(sewtype)'"
				local ddd = 1
				tempname ss

				scalar `ss' = 1/(`ddd'*`ddd')
				mat `V' = `ss'*`V'*`vm'*`V''
			}

			if "`jknife1'" != "" {
				JK1irls `b' `V' `y' "`xvars'" `h' `eta' /*
					*/ `mu' `v' `dmu' "`tvar'" `z' /*
					*/ "`moffset'" "`constant'"  /*
					*/ `wt' `touse' "`weight'"
				local vce    "1-step JKnife"
				local setype "One-step Jackknife"
				local p = colsof(`b')
				tempname ss
				scalar `ss' = (`nobs'-`p')/`nobs'
				mat `V' = `ss'*`V'
			}

			if "`jknife'" != "" {
				noi JKnife `b' `V' `y' "`xvars'" "`tvar'" /*
					*/ "`offopt'" "`constant'" /*
					*/ "`argfam'" "irls" "`arglink'" /*
					*/ "[`weight'`exp']" `touse' /*
					*/ "`cluster'" `dots' "`weight'" "`wt'"
				local mm = `nobs'
				local vce    "Jackknife"
				local setype "Jackknife"
				if "`cluster'" != "" { local mm = $SGLM_nc }
				else local mm=1
				local p = colsof(`b')
				local nnn = `nobs'-$SGLM_bm
				tempname ss
				scalar `ss' = (`nnn'-`p')/(`nnn'*`mm') 
				mat `V' = `ss'*`V'
			}
			if "`bstrap'" != "" {
				tempname b V
				mat `b' = e(b)
				mat `V' = e(V)
				noi Bstrap `b' `V' `y' "`xvars'" "`offopt'" /*
					*/ "`constant'" "`argfam'" "irls" /*
					*/ "`arglink'" "[`weight'`exp']" /*
					*/ `touse' `brep'       /*
					*/ "`clopt'" `dots'
				local vce    "Bootstrap"
				local setype "Bootstrap"
				if "`cluster'" != "" { local mm = $SGLM_nc }
				else local mm=1
				local p = colsof(`b')
				tempname ss
				scalar `ss' = (`nobs'-`p')/(`nobs'*(`brep'-$SGLM_bm))
				mat `V' = `ss'*`V'
			}

					/* get rid of `Wscale' scaling */

			if "`opg'`nwnam'`robust'`jknife'`jknife1'`bstrap'" == ""  {
				scalar `Wscale' = 1/`Wscale'
				mat `V' = `Wscale'*`V' 
			}

			if `"`zapse'"'=="yes" {
				local i 1
				while `i'<=rowsof(`V') {
					mat `V'[`i',`i'] = 0
					local i=`i'+1
				}
			}

			if `vfactor' != 1 {
				mat `V' = `vfactor'*`V'
			}

			tempvar mysamp
			local k = colsof(`b')
			qui gen byte `mysamp' = `touse'
			est post `b' `V', depname(`y') obs(`nobs') /*
				*/ `dofopt' esample(`mysamp')
			est local depvar "`y'"
			est local wtype "`weight'"
			est local wexp  "`exp'"
			est scalar rc = `rc'
			est scalar disp = `disp'
			est scalar k = `k'
			est local crittype "deviance"
		}
	}

	mac drop SGLM_running SGLM_nonstan

	if "`score'" != "" {
		label var `scvar' "Score index from glm"
		rename `scvar' `score'
		est local scorevars `score'
	}

	local dfm : word count `xvars'
	local p  = `dfm' + (`"`constant'"'=="")
	local df = `nobs'-`p'-`df_a'

	est local link      "$SGLM_L"
	est local varfunc  "$SGLM_V"
	est local m	   "$SGLM_m"
	est local a	   "$SGLM_a"
	est local msg      "`cd'"
	est local cons	   "`constant'"
	est local oim	   "`oim'"
	est local clustvar "`cluster'"
 	capture confirm number $SGLM_p
	if !_rc {
		est scalar power = $SGLM_p
	}	
	else {
		est scalar power = 0
	}
	
	if "$SGLM_V" == "glim_v2" & "$SGLM_L" == "glim_l01" & "`eform'" == "" {
		est local msg2 "Coefficients are the risk differences"
	}
	if "`irls'" != "" {
		if "`oim'" == "" { 
			est local opt  "MQL Fisher scoring"
			est local opt2 "(IRLS EIM)" 
		}
		else  { 
			est local opt  "MQL Newton-Raphson"
			est local opt2 "(IRLS OIM)" 
		}
	}
	else {
		est local opt = "ML: Newton-Raphson"
		est local opt2 
	}
	if "$SGLM_bm" != "" {
		if $SGLM_bm > 0 { est scalar Nf = $SGLM_bm }
	}
	est scalar df_m = `dfm'
	est scalar df = `df'
	if "$SGLM_nc" != "" {
		est scalar N_clust = $SGLM_nc
	}
	if `brep' > 0 {
		est scalar N_brep  = `brep'
	}

	est scalar vf             = `vfactor'
	est scalar phi            = `disp'
	est scalar deviance       = abs(`newdev')
	est scalar dispers        = e(deviance)/e(df) 
	est scalar deviance_s     = e(deviance)/e(phi)
	est scalar dispers_s      = e(deviance_s)/e(df)
	est scalar deviance_p     = abs(`chi2')
	est scalar dispers_p      = e(deviance_p)/e(df)
	est scalar deviance_ps    = e(deviance_p)/e(phi)
	est scalar dispers_ps     = e(deviance_ps)/e(df)

	capture drop `v'
	qui gen double `v' = sum((`y'-`mu')^2) if `touse'
	est scalar bic = e(deviance) - `df'*log(`nobs')

	est local varfuncf "$SGLM_vf"
	est local varfunct "$SGLM_vt"
	est local linkf "$SGLM_lf"
	est local linkt "$SGLM_lt"
	est local vcetype "`vce'"

	est local se1  "`setype'"
	est local se2  "`sewtype'"

	est local offset  "`offvar'"
	est local offset1 /* erase; set by -ml- */
	est local predict "glim_p"
	est local cmd     "glm"

	capture mac drop SGLM_bm SGLM_nc

	if "`display'"=="" {
		Display, `eform' level(`level') `header'
	}
	error `e(rc)'
end

program Display
	di 
	syntax [, noHEADer Level(int $S_level) eform]
	if `level' < 10 | `level' > 99 {
		di as err "level() must be between 10 and 99"
		exit 198
	}
	
	if "`eform'" != "" {
		Eform 
		local eform  "`r(eform)'"
	}
	if "`eform'" != "" {
		local eopt "eform(`eform')"
	}

	if "`header'" == "" { Head }
	est di, level(`level') `eopt' first
	if "`e(msg)'" != "" {
		noi di as txt "(Standard errors scaled using `e(msg)')"
	}
	if "`e(msg2)'" != "" {
		noi di as txt "`e(msg2)'"
	}
	if "`e(Nf)'" != "" {
		if `e(Nf)' > 1 { local s "s" }
		noi di as txt "(" as res `e(Nf)' as txt /*
			*/ " model`s' failed to converge)"
	}
end

program Head
/*
----+----1----+----2----+----3----+----4----+----5----+----6----+----7----+---
Generalized linear models                          No. of obs      = ######### 
Optimization     :                                 Residual df     = #########  
                                                   Scale parameter = #########
Deviance         = #########                       (1/df) Deviance = #########
Pearson          = #########                       (1/df) Pearson  = #########

Variance function: V(u) = ########################    ########################
Link function    : g(u) = ########################    ########################
Standard errors  : ###############################
                   ###############################

Log likelihood   = #########                       AIC             = #########
BIC              = #########                       APC             = #########
*/
	di as txt "Generalized linear models" /*
		*/ _col(52) "No. of obs"        _col(68) "=" /*
		*/ _col(70) as res %9.0g e(N)
	di as txt "Optimization     : " as res "`e(opt)'" /*
		*/ as txt _col(52) "Residual df" _col(68) "=" /*
		*/ _col(70) as res %9.0g e(df)

	di as res _col(20) "`e(opt2)'" as txt _col(52) "Scale parameter"  /*
		*/ _col(68) "="  _col(70) as res %9.0g e(phi)

	di as txt "Deviance" _col(18) "=" as res _col(20) %12.0g e(deviance) /*
		*/ as txt _col(52) "(1/df) Deviance" /*
		*/ _col(68) "=" as res _col(70) %9.0g e(dispers)
	di as txt "Pearson" _col(18) "=" as res _col(20) %12.0g e(deviance_p) /*
		*/ as txt _col(52) "(1/df) Pearson" /*
		*/ _col(68) "=" as res _col(70) %9.0g e(dispers_p)

	di
	di as txt "Variance function: " as res "V(u) = " /*
		*/ as res _col(27) "`e(varfuncf)'" /*
		*/ _col(52) as txt "[" as res "`e(varfunct)'" as txt "]"
	di as txt "Link function    : " as res "g(u) = " /*
		*/ as res _col(27) "`e(linkf)'" /*
		*/ _col(52) as txt "[" as res "`e(linkt)'" as txt "]"

	di as txt "Standard errors  : " as res "`e(se1)'"
	if "`e(se2)'" != "" {
		di as res _col(20) "`e(se2)' weights"
	}
	di

	if "`e(ll)'" != "" {
		local crtype = upper(substr(`"`e(crittype)'"',1,1)) + /*
			*/ substr(`"`e(crittype)'"',2,.)
		local crlen = max(18,length(`"`crtype'"') + 2)
		di as txt "`crtype'" _col(`crlen') "= " /*
			*/ as res %12.0g e(ll)  /*
			*/ as txt _col(52) "AIC" _col(68) "=" /*
			*/ as res _col(70) %9.0g e(aic)
	}
	else	local crlen 18
	di as txt "BIC" _col(`crlen') "=" as res _col(20) %12.0g e(bic) 
	else if "`e(disp)'" != "" & "`e(disp)'" != "1" {
		di as txt "Quasi-likelihood model with dispersion: " /*
			*/ as res `e(disp)'
	}
	di
end

program MapNW, rclass /* f [lag] */
	args nobs f lag

        local f = lower(trim(`"`f'"'))
        local l = length(`"`f'"')

	local s1
	local s2

        if `"`f'"'==""                                   { local s1 "" }
        else if `"`f'"'==substr("nwest",1,max(`l',2))    { local s1 "glim_nw1" }
        else if `"`f'"'==substr("bartlett",1,max(`l',2)) { local s1 "glim_nw1" }
        else if `"`f'"'==substr("gallant",1,max(`l',2))  { local s1 "glim_nw2" }
        else if `"`f'"'==substr("parzen",1,max(`l',2))   { local s1 "glim_nw2" }
        else if `"`f'"'==substr("anderson",1,max(`l',2)) { local s1 "glim_nw3" }
	else                                             { local s1 "`f'" }

	if "`lag'" != "" {
		confirm integer number `lag'
		if `lag' < 0 {
			noi di as err "Newey-West lag must be positive"
			exit 198
		}
	}
	else local lag = `nobs'-2

	if `lag' >= `nobs'-2 {
		local lag = `nobs'-2
	}
	local s2 `lag'

	ret local nwnam "`s1'"
	ret local nwlag "`s2'"
end

program MapFL, rclass /* family link */
        args f ulink

        MapFam `f'                      /* map user-specified family    */
        local fam `"`r(famcode)'"'      /* store program in fam         */

        local mfixed 1
        local m 1
        local k 1
        if `"`fam'"'=="glim_v2" {           /* bin takes an optional argument */
                tokenize `"`f'"'
                if `"`2'"'!="" {
                        capture confirm number `2'
                        if _rc {
                                unabbrev `2'
                                local m `"`s(varlist)'"'
                                local mfixed 0
                        }
                        else {
                                if `2'>=. | `2'<1 {
                                        di as err /*
                                     */ `"`2' in family(binomial `2') invalid"'
                                        exit 198
                                }
                                local m `2'
                        }
                }
		if `"`3'"'!="" {
			di as err "family(`f') invalid"
			exit 198
		}
        }
        if `"`fam'"'=="glim_v6" {    /* nb takes an optional argument */
                tokenize `"`f'"'
                if `"`2'"' != "" {
                        confirm number `2'
                        if `2'<=0 {
                                di as err /*
                                */ `"`2' in family(nbinomial `2') invalid"'
                                exit 198
                        }
                        local k `2'
                }
		if `"`3'"'!="" {
			di as err "family(`f') invalid"
			exit 198
		}
        }
	if `"`fam'"'~="glim_v1" & `"`fam'"'~="glim_v2" & /*
	*/ `"`fam'"'~="glim_v3" & `"`fam'"'~="glim_v4" & /*
	*/ `"`fam'"'~="glim_v5" & `"`fam'"'~="glim_v6" {
		tokenize `"`f'"'
		if `"`2'"' != "" {
			global SGLM_fa `2' 
		}                 
		if `"`3'"'!="" {
			di as err "family(`f') invalid"
			exit 198
		}
	}        
	MapLink `ulink'
        local link `"`r(link)'"'
        local pow `"`r(power)'"'
        local scale1 = `"`fam'"'=="glim_v2" | /*
			*/ `"`fam'"'=="glim_v3" | `"`fam'"'=="glim_v6"

        if `"`link'"'=="" {             /* apply defaults */
                local link "glim_l11"
                if `"`fam'"'=="glim_v1"      { local link "glim_l01"
					       local pow 1           }
                else if `"`fam'"'=="glim_v2" { local link "glim_l02" 
					       local pow 0           }
                else if `"`fam'"'=="glim_v3" { local link "glim_l03"
					       local pow  0          }
                else if `"`fam'"'=="glim_v4" { local link "glim_l09"
					       local pow -1          }
                else if `"`fam'"'=="glim_v5" { local link "glim_l10"
					       local pow -2          }
                else if `"`fam'"'=="glim_v6" { local link "glim_l03"
					       local pow  0          }
                else local link
        }

        ret local family `"`fam'"'
        ret local link   `"`link'"'
        ret local power  `pow'
        ret local scale  `scale1'
        ret local m      `m'
        ret local mfixed `mfixed'
        ret local k      `k'
end


program MapFam, rclass /* <code> */
	args f

        local s1
        local f = lower(trim(`"`f'"'))
        local l = length(`"`f'"')

        if `"`f'"'==""                                    { local s1 "1" }
        else if `"`f'"'==substr("gaussian",1,max(`l',3))  { local s1 "1" }
        else if `"`f'"'==substr("normal",1,`l')           { local s1 "1" }
        else if `"`f'"'==substr("binomial",1,`l')         { local s1 "2" }
        else if `"`f'"'==substr("bernoulli",1,`l')        { local s1 "2" }
        else if `"`f'"'==substr("poisson",1,`l')          { local s1 "3" }
        else if `"`f'"'==substr("gamma",1,max(`l',3))     { local s1 "4" }
        else if `"`f'"'==substr("igaussian",1,max(`l',2)) { local s1 "5" }
        else if `"`f'"'==substr("inormal",1,max(`l',2))   { local s1 "5" }
        else if `"`f'"'=="ivg"                            { local s1 "5" }
        else if `"`f'"'==substr("nbinomial",1,max(2,`l')) { local s1 "6" }

	if "`s1'" != "" {
		local s1 "glim_v`s1'"
	}
	else {
		local s1 "`f'"
	}
        ret local famcode `s1'
end

program MapLink, rclass /* <code> <coderest> ... */
	args ulink upow rest
        if `"`rest'"'!="" {
                di as err "link(`ulink' `upow' `rest') invalid"
                exit 198
        }

        local ulink = lower(`"`ulink'"')
        local l     = length(`"`ulink'"')
        local s1      "11"                /* power(a) link */
        local s2      .                   /* a of power(a) */

        if `"`ulink'"'==""                              { local s1 "" }
        else if `"`ulink'"'==substr("identity",1,`l')   { local s1 "01" 
							  local s2 1    }
        else if `"`ulink'"'==substr("reciprocal",1,`l') { local s1 "09"
							  local s2 -1   }
        else if `"`ulink'"'=="log"                      { local s1 "03"
							  local s2 0    }
        else if `"`ulink'"'==substr("power",1,max(`l',3)) {
                capture confirm number `upow'
                if _rc {
                        di as err "invalid # in link(power #)"
                        exit 198
                }
		if `upow' == 0 {
			local s1 "03"
		}
		if `upow' == -1 {
			local s1 "09"
		}
		else if `upow' == -2 {
			local s1 "10"
		}
                local s2 `upow'
                local upow
        }
        else if `"`ulink'"'==substr("opower",1,max(`l',3)) {
                capture confirm number `upow'
                if _rc {
                        di as err "invalid # in link(opower #)"
                        exit 198
                }
                if `upow'==0 {
                        local s1 "02"	/* logit */
                }
                else {
                        local s1 "12"   /* odds-power */
                }
		local s2 `upow'
                local upow
        }
        else if `"`ulink'"'==substr("logit",1,`l')          { local s1 "02" 
							      local s2 0    }
        else if `"`ulink'"'==substr("nbinomial",1,max(`l',2)) { local s1 "04" }
        else if `"`ulink'"'==substr("logc",1,4)             { local s1 "05" }
        else if `"`ulink'"'==substr("loglog",1,max(`l',4))  { local s1 "06" }
        else if `"`ulink'"'==substr("cloglog",1,`l')        { local s1 "07" }
        else if `"`ulink'"'==substr("probit",1,`l')         { local s1 "08" }

        else {
		local s1
		local s2 "`upow'"
        }

	if `"`s1'"' != "" {
		local s1 "glim_l`s1'"
	}
	else {
		local s1 "`ulink'"
	}

        ret local link  `s1'
        ret local power `s2'
end

program Eform , rclass
        local var "`e(varfunc)'"
        local lnk "`e(link)'"

        local eform "ExpB"
        if ("`var'" == "glim_v3" | "`var'" == "glim_v6") {
                if "`lnk'" == "glim_l03" {
                        local eform "IRR"
                }
        }
        else if "`var'" == "glim_v2" {
                if "`lnk'" == "glim_l02" {
                        local eform "Odds Ratio"
                }
                else if "`lnk'" == "glim_l03" {
                        local eform "Risk Ratio"
                }
                else if "`lnk'" == "glim_l05" {
                        local eform "Hlth Ratio"
                }
                else local eform "ExpB"
        }
        return local eform "`eform'"
end

program NWest, rclass
	args nwnam nwlag wv res touse constant xvars

	local xv "`xvars'"
	if "`constant'" == "" {
		tempvar CONS
		gen byte `CONS' = 1
		local xn "`xv' _cons"
		local xv "`xv' `CONS'"
	}
	else {  local xn "`xv'" }
	local nx : word count `xv'

	qui summ `touse', meanonly
	local nobs = r(N)

	tempvar vt1
	tempname tt xtx tx

	gen double `vt1' = . 
	local j 0
	while `j' <= `nwlag' {
		capture mat drop `tt'
		local i 1
		while `i' <= `nx' {
			local x : word `i' of `xv'
			qui replace `vt1' = `x'[_n-`j']*`res'*`res'[_n-`j']* /* 
					*/ `wv'[_n-`j'] if `touse'
			mat vecaccum `tx' = `vt1' `xv' if `touse', nocons
			mat `tt' = nullmat(`tt') \ `tx'
			local i = `i'+1
		}
		`nwnam' `nwlag' `j' 
		if `r(wt)'>=. {
			noi di as err "Newey-West weights are missing"
			exit 199
		}
		mat `tt' = (`tt' + `tt'')*r(wt)
		if `j' > 0 {
			mat `xtx' = `xtx'+`tt'
		}
		else {
			mat `xtx' = `tt'*.5
		}
		local j = `j'+1
	}
	mat colnames `xtx' = `xn'
	mat rownames `xtx' = `xn'
	return matrix V `xtx'
end

program LoadXi
	args beta xi i 

	local p = colsof(`beta')
	local xv : colnames(`beta')
	local j 1
	while `j'<=`p' {
		local x : word `j' of `xv'
		if "`x'" != "_cons" {
			mat `xi'[1,`j'] = `x'[`i']
		}
		else    mat `xi'[1,`j'] = 1 
		local j = `j'+1
	}
end

program JK1nr
	args b V y xvars eta mu v dmu tvar z moffset constant wt touse weight 
	capture drop `eta'
	capture drop `mu'
	capture drop `v'
	capture drop `dmu'
	_predict `eta', xb
	$SGLM_L 1 `eta' `mu' 
	$SGLM_V 1 `eta' `mu' `v'
	$SGLM_L 2 `eta' `mu' `dmu'
	capture drop `z'
	tempvar W h
	tempname ll 
	gen double `z' = `eta' + (`y'-`mu')/`dmu' `moffset' if `touse'
	gen double `W' = `dmu'*`dmu'/`v' if `touse'
	estimates hold `ll'
	cap regress `z' `xvars' [iw=`W'*`wt'], mse1 `constant'
	if e(N)>=. {exit 459}
	if _rc {exit _rc}
	_predict double `h' if `touse', stdp 
	replace `h' = `h'*`h'*`v'
	estimates unhold `ll'

	capture drop `z'
	gen double `z' = (`mu'-`y')/(1-`h') 
	sort `touse' `tvar'
	capture drop `dmu'

	gen `dmu' = _n if `touse'
	local p = colsof(`b')

	tempname bi xtxi xi ss
	mat `xi' = J(1,`p',0)
	mat `V' = 0*`V'
	if "`weight'" == "fweight" {
		mat accum `xtxi' = `xvars' [iw=`v'*`wt'] if `touse', `constant'
	}
	else    mat accum `xtxi' = `xvars' [iw=`v'] if `touse', `constant'
	mat `xtxi' = syminv(`xtxi')
	summ `dmu'
	local factor 1
	local max = r(max)
	local i   = r(min)
	while `i' <= `max' {
		if !`touse'[`i'] {
			local i = `i' + 1
			continue
		}	
		LoadXi "`b'" "`xi'" `i'
		mat `bi' = `xi'*`xtxi'
		scalar `ss' = `z'[`i']
		mat `bi' = `ss'*`bi'
		if "`weight'" == "fweight" { local factor = `wt'[`i'] }
		mat `V' = `V' + `factor'*(`bi')'*`bi' 
		local i = `i'+1
	}
end

program JK1irls
	args b V y xvars h eta mu v dmu tvar z moffset constant wt touse weight
	capture drop `v'
	capture drop `dmu'
	capture drop `z'
	$SGLM_V 1 `eta' `mu' `v'
	gen double `z' = (`mu'-`y')/(1-`h')
	sort `touse' `tvar'
	gen `dmu' = _n if `touse'
	local p = colsof(`b')

	tempname bi xtxi xi ss
	mat `xi' = J(1,`p',0)
	mat `V' = 0*`V'
	mat accum `xtxi' = `xvars' [iw=`v'*`wt'] if `touse', `constant'
	mat `xtxi' = syminv(`xtxi')
	summ `dmu'
	local factor 1
	local max = r(max)
	local i = r(min)
	while `i' <= `max' {
		if !`touse'[`i'] {
			local i = `i' + 1
			continue
		}
		LoadXi "`b'" "`xi'" `i'
		mat `bi' = `xi'*`xtxi'
		scalar `ss' = `z'[`i']
		mat `bi' = `ss'*`bi'
		if "`weight'" == "fweight" { local factor = `wt'[`i'] }
		mat `V' = `V' + `factor'*(`bi')'*`bi'
		local i = `i'+1
	}
end

program JKnife
	args b V y xvars tvar offset constant family irls link /*
		*/ wtopt touse cluster dots weight wt

	if "`dots'" != "" { noi di } 
	sort `touse' `tvar'
	local p = colsof(`b')
	mat `V' = 0*`V'

	if "`irls'" == "" {
		local iopt "from(`b')"
	}
	else {
		tempvar xb mu
		qui _predict double `xb' if `touse', xb
		qui $SGLM_L 1 `xb' `mu'
		local iopt "mu(`mu')"
	}
	local reg = "$SGLM_V"=="glim_v1" & "$SGLM_L"=="glim_l01" & "`irls'"==""
	if `reg' { local iopt }
	tempname bi ll
	estimates hold `ll'
	tempvar dd
	if "`cluster'" != "" {
		qui egen `dd' = group(`cluster') if `touse'
		qui summ `dd'
		local nsize = r(max)
	}
	else    {
		qui gen  `dd' = _n if `touse'
		qui summ `dd'
		local nsize = r(N)
	}
	local max = r(max)
	local min = r(min) 
	local i   = `min'
	local nb = 0

	if "`dots'" != "" {
		di as txt "Jackknife iterations (" /* 
			*/ as res `nsize' as txt ")"
		di as txt "{hline 4}{c +}{hline 3} 1 "/*
			*/"{hline 3}{c +}{hline 3} 2 "/*
			*/"{hline 3}{c +}{hline 3} 3 "/*
			*/"{hline 3}{c +}{hline 3} 4 "/*
			*/"{hline 3}{c +}{hline 3} 5" 
	}

	if "`weight'" == "fweight" {
		local wtopt "[fweight=`wt']"
		local fixwt 1
	}
	else    local fixwt 0


	local nx : word count `xvars'

	while `i' <= `max' {
		if "`cluster'"=="" & !`touse'[`i'] {
			local i = `i'+1
			continue
		}
		if `fixwt' == 0 {
			qui _rmcoll `xvars' `wtopt' if `touse' & `dd'!=`i', `constant'
			local nxrm : word count `r(varlist)'
		}
		else    local nxrm = `nx'

		if `nxrm' == `nx' {
			capture {
				if `fixwt' { 
					qui summ `wt' if `dd'==`i', meanonly
					local factor = r(mean) 
					qui replace `wt'=`wt'-1 if `dd'==`i' 
					qui glm `y' `xvars' if `touse' /*
						*/ `wtopt' /*
						*/ , `family' /*
						*/ `link' `offset' /*
						*/ `constant' `irls' `iopt' 
				} 
				else {
					local factor 1
					qui glm `y' `xvars' if `touse' & /*
					*/ `dd'!=`i' `wtopt' /*
					*/ , `family' /*
					*/ `link' `offset' /*
					*/ `constant' `irls' `iopt' 
				}
				mat `bi' = e(b)
				mat `V' = `V' + `factor'*(`bi'-`b')'*(`bi'-`b')
			}
			if _rc {
				if _rc == 1 { exit _rc } 
				local err 1
				local nb = `nb'+`factor'
			}
			else    local err 0
			if `fixwt' {
				qui replace `wt'=`wt'+1 if `dd'==`i' 
			}
		}
		else    local err 1
		if "`dots'" != "" {
			local j = `i'-`min'+1
			if `err' {
				noi di as err "x" _c 
			}
			else {
				noi di as txt "." _c 
			}
			if int(`j'/50)*50 == `j' {
				noi di as res " " `j' 
			}
		}
		local i = `i'+1
	}
	estimates unhold `ll'
	global SGLM_nc = `nsize'
	global SGLM_bm = `nb'
	if "`dots'" != "" { noi di } 
end

program Bstrap
	args b V y xvars offset constant family irls link /*
		*/ wtopt touse brep clopt dots

	if "`dots'" != "" { noi di } 
	sort `touse' `tvar'
	local p = colsof(`b')
	mat `V' = 0*`V'

	if "`irls'" == "" {
		local iopt "from(`b')"
	}
	else {
		tempvar xb mu
		qui _predict double `xb' if `touse', xb
		qui $SGLM_L 1 `xb' `mu'
		local iopt "mu(`mu')"
	}

	local reg = "$SGLM_V"=="glim_v1" & "$SGLM_L"=="glim_l01" & "`irls'"==""
	if `reg' { local iopt }

	tempname bi ll
	preserve
	qui drop if `touse'==0
	tempfile tmp
	estimates hold `ll'
	qui save `tmp'
	local nb = 0 
	local i 1

	if "`dots'" != "" {
		di as txt "Bootstrap iterations (" /* 
			*/ as res `brep' as txt ")"
		di as txt "{hline 4}{c +}{hline 3} 1 "/*
                        */"{hline 3}{c +}{hline 3} 2 "/*
                        */"{hline 3}{c +}{hline 3} 3 "/*
                        */"{hline 3}{c +}{hline 3} 4 "/*
                        */"{hline 3}{c +}{hline 3} 5"
	}

	local nx : word count `xvars'

	while `i' <= `brep' {
		qui use `tmp', clear
		qui bsample , `clopt'

		qui _rmcoll `xvars' `wtopt', `constant'
		local nxrm : word count `r(varlist)'

		local err 0
		if `nxrm' == `nx' {
			capture {
				qui glm `y' `xvars' `wtopt' /*
					*/ , `family' /*
					*/ `link' `offset' /*
					*/ `constant' `irls' `iopt'
				mat `bi' = e(b)
				mat `V' = `V' + (`bi'-`b')'*(`bi'-`b')
			}
			if _rc {
				local err 1
				if _rc == 1 { exit _rc } 
				local nb = `nb'+1
			}
			else    local err 0
		}
		else local err 1
		if "`dots'" != "" {
			if `err' {
				noi di as err "x" _c
			}
			else {
				noi di as txt "." _c 
			}
			if int(`i'/50)*50 == `i' {
				noi di as res " " `i' 
			}
		}
		local i = `i'+1
	}
	restore
	estimates unhold `ll'
	global SGLM_bm = `nb'
	if "`dots'" != "" { noi di } 
end
exit