*! version 3.0.1 PR 15-Feb-94.  STB-18: sg22
* Based on J. Hilbe, STB-16  sg16.3, 11/93.
program define glmr
	version 3.1
	if "`*'"=="" | substr("`1'",1,1)=="," { 
		_glmrpt `*'
		exit
	}

	local varlist "req ex"
	local options "LEvel(int $S_level) EForm"
	#delimit ;
	local options "`options' LTol(real -1) ITerate(int 0) Family(string)
Link(string) K(real 1) DISP(real 1) Scale(string) FRVars Offset(string)
LNOffset(string) noCONStant INIt(string) noLOg" ;
	#delimit cr
	local in "opt"
	local if "opt"
	local weight "aweight fweight"
	parse "`*'"
	if "`init'"!="" { conf var `init' }

	parse "`varlist'", parse(" ")

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

	if "`constan'"!="" { local cons "nocons" }
	local small 1e-6
	if `level'<10 | `level'>99 { local level 95 }
	if `ltol'<0 { local ltol 1e-8 }
	if `iterate'<=0 { local iterate 50 }
/*
	Validate family and set up link indicator
*/
	_jprglld "`family'" "`link'" `k'
	local family "$S_1"
	local link "$S_2"
	local pow $S_3
	local scale1 $S_4 /* 1 for families with default scale parameter 1 */
	local m $S_5
	local mfixed $S_6
	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'"!="" {
		_jprglef `family' "`link'"
		local eform "$S_1"
	}
 	tempvar mu eta W z dev chi2 wt touse tmp
	quietly {
		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
			_crcunab `lnoffse'
			local offstr "ln($S_1)"
			gen `offvar' = `offstr'
		}
		else if "`offset'"!="" {
			_crcunab `offset'
			local offvar "$S_1"
			local offstr "$S_1"
		}
		if "`offstr'"!="" {
			local offset "$S_1"
			local offmiss "| `offvar'==."
			local offopt ", offset `offstr'"
		}

/*
TO DELETE
		gen byte `touse' = 1 `if' `in'
		replace `touse' = 0 if `y'==. | `touse'==. `mmiss' `offmiss'
		if "`exp'"!="" {
			gen `wt' `exp' if `touse'
			replace `wt' = . if `wt'<=0
			replace `touse' = 0 if `wt'==.
			if "`weight'"=="aweight" {
				summ `wt'
				replace `wt' = `wt' / _result(3)
			}
		}
		else gen `wt' = `touse' if `touse'
END OF TODELETE
*/
		mark `touse' [`weight' `exp'] `if' `in'
		markout `touse' `y' `*' 
		replace `touse'=0 if 0 `mmiss' `offmiss'
		if "`exp'"!="" { 
			gen `wt' `exp' if `touse'
			replace `touse'=0 if `wt'==.  /* bug in mark ! */
			if "`weight'"=="aweight" {
				summ `wt'
				replace `wt' = `wt'/_result(3)
				local wgt "[aw=`wt']"
			}
			else	local wgt "[fw=`wt']"
		}
		else gen `wt' = `touse' if `touse'

		sum `y' if `touse'		/* purposefully not weighted */
		if _result(1)==0 { noisily error 2000 } 
		if _result(1)==1 { noisily error 2001 }
		if _result(5)==_result(6) {
			noi di _n in bl "outcome does not vary."
			exit 2002
		}

				/* 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
*/
		regress `y' `*' `wgt' if `touse', `cons'
		local nobs = _result(1)
		if "`init'"!="" {
			gen `mu' = `init' if `touse'
		}
		else {
			sum `y' `wgt' if `touse'
			gen `mu' = (`y'+_result(3))/(`m'+1)
		}
		gen `W' = 0 if `touse' 
		gen `z' = 0 if `touse'  
		gen `dev' = 0 if `touse'
		gen `tmp' = .
/*
	Initialise linear predictor
*/
		if "`link'"!="pow" {
			if "`family'"=="bin" {
				if "`init'"=="" {
					replace `mu' = `m'*(`y'+.5)/(`m'+1) /*
					 */ if `touse'
				}
				if "`link'"=="l" {
					gen `eta' = ln(`mu'/(`m'-`mu'))
				}
				else if "`link'"=="p" {
					gen `eta' = invnorm(`mu'/`m')
				}
				else if "`link'"=="c" {
					gen `eta' = ln(-ln(1-`mu'/`m'))
				}
				if "`link'"=="opo" {
					gen `eta' = cond(abs(`pow')<`small', /*
					 */ ln(`mu'/(`m'-`mu')), /*
					 */ ((`mu'/(`m'-`mu'))^`pow'-1)/`pow')
				}
			}
			if "`family'"=="gau" {
				gen `eta' = `mu'
			}
			else if "`family'"=="poi" {
				gen `eta' = ln(`mu')
			}
			else if "`family'"=="gam" {
				gen `eta' = 1/`mu'
			}
			else if "`family'"=="ivg" {
				gen `eta' = 1/`mu'^2
			}
			else if "`family'"=="nb" {
				gen `eta' = ln(`mu'/(`mu'+`k'))
			}
		}
		else {
			if abs(`pow')<`small' {
				gen `eta' =  ln(`mu')
			}
			else gen `eta' = `mu'^`pow'
		}
/*
	Local scoring algorithm: IRLS
*/
		local oldev 0
		local newdev 1  
		local i 1
		`iflog' di
		while abs(`newdev'-`oldev')>`ltol' & `i'<=`iterate' {
			local 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.
*/
			_jprglwz `touse' `y' `family' "`link'" `pow' `m' /*
			 */ `disp' `k' `eta' `mu' `W' `z' "`offvar'"
/*
	Weighted regression
*/
			regress `z' `*' [iw=`W'*`wt'], /*
			 */ mse1 dof(100000) dep(`y') `cons'
/*
	Get linear predictor
*/
			cap drop `eta'
			predict `eta' if `touse' 
			if "`offstr'"!="" { replace `eta' = `eta'+`offvar' }
/*
	Get inverse link
*/
			global S_E_fam "`family'" 
			global S_E_link "`link'"
			_jprglil `eta' `mu' `pow' `k'
			if "`family'"=="bin" & !`bernoul' {
				replace `mu' = `mu'*`m'
			}
*noi di "W, z, mu, eta"
*noi sum `W' `z' `mu' `eta'
/*
	Deviance
*/
			if "`family'"=="gau" {
				replace `dev' = (`y'-`mu')^2
			}
			else if "`family'"=="bin" & `bernoul' {
				replace `dev'= cond(`y', -2*ln(`mu'), /*
				 */ -2*ln(1-`mu'))
			}
			else if "`family'"=="bin" & !`bernoul' {
				replace `dev' = cond(`y'>0 & `y'<`m', /*
	 */ 2*`y'*ln(`y'/`mu') + 2*(`m'-`y')*ln((`m'-`y')/(`m'-`mu')), /*
	 */ cond(`y'==0, 2*`m'*ln(`m'/(`m'-`mu')), 2*`y'*ln(`y'/`mu')) )
*				replace `dev' = 2*(`m'-`y')*ln((`m'-`y')/ /*
*				 */ (`m'-`mu')) if `y'==0
*				replace `dev' = 2*`y'*ln(`y'/`mu') /*
*				 */ if `y'==`m'
*				replace `dev' = 2*`y'*ln(`y'/`mu')+ /*
*				 */ 2*(`m'-`y')*ln((`m'-`y')/ /*
*				 */ (`m'-`mu')) if `y'!=0 & `y'!=`m'
			}
			else if "`family'" == "poi" { 
				replace `dev' = cond(`y'==0, 2*`mu', /*
				 */ 2*(`y'*ln(`y'/`mu')-(`y'-`mu')))
			}
			else if "`family'" == "gam" {
				replace `dev' = -2*(ln(`y'/`mu')- /*
				 */ (`y'-`mu')/`mu')
			}
			else if "`family'" == "ivg" {
				replace `dev' = (`y'-`mu')^2/(`mu'^2*`y')
			}
			else if "`family'" == "nb" {
*				replace `dev' = 2*(`y'*ln(`y'/`mu')- /*
*				 */ (1+`k'*`y')/`k'*  /*
*				 */ ln((1+`k'*`y')/(1+`k'*`mu'))) if `y'>0
*				replace `dev' = 2*ln(1+`k'*`mu')/`k' if `y'==0
				replace `dev' = cond(`y'==0, /*
				 */ 2*ln(1+`k'*`mu')/`k', /*
				 */ 2*(`y'*ln(`y'/`mu')-(1+`k'*`y')/`k'*  /*
				 */ ln((1+`k'*`y')/(1+`k'*`mu'))))
			}
			replace `tmp' = sum(`dev'*`wt')
			local newdev = `tmp'[_N]/`disp'
			`iflog' di in gr "Iteration `i' : deviance = " /*
			 */ in ye %9.4f `newdev'
			local i = `i'+1
		}
/*
	Calculate chi2 and dispersions
*/
		if "`family'"=="gau" {
			gen `chi2' = `disp'*(`y'-`mu')^2
		}
		else if "`family'"=="bin" {
			gen `chi2' = (`y'-`mu')^2/(`mu'*(1-`mu'/`m')*`disp')
		}
		else if "`family'"=="poi" {
			gen `chi2' = (`y'-`mu')^2/(`mu'*`disp')
		}
		else if "`family'"=="gam" {
			gen `chi2' = (`y'-`mu')^2/(`mu'^2*`disp')
		}
		else if "`family'"=="ivg" {
			gen `chi2' = (`y'-`mu')^2/(`mu'^3*`disp')
		}
		else if "`family'"=="nb" {
			gen `chi2' = (`y'-`mu')^2/(`disp'*`mu'*(1+`k'*`mu'))
		}
		replace `chi2' = sum(`chi2'*`wt')
		local chisq = `chi2'[_N]
		local df = `nobs'-_result(3)-("`cons'"=="")
		local dispc = `chisq'/`df'
		local dispd = `newdev'/`df'
	}
/*
	Store results
*/

	_jprglfl "`family'" "`link'" `pow' `bernoul' `k'
	local o "$S_1, $S_2`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 qui regress `z' `*' [iw=`W'*`wt'/`delta'], /*
		 */ mse1 dof(`dof') dep(`y') `cons'
	} 

	tempname b V
	if `reg' {
		local dofopt "dof(`df')"
	}
	mat `b' = get(_b)
	mat `V' = get(VCE)
	if "`zapse'"=="yes" { 
		local i 1
		while `i'<=rowsof(`V') {
			mat `V'[`i',`i'] = 0 
			local i=`i'+1
		}
	}
	mat post `b' `V', depname(`y') obs(`nobs') `dofopt'

/* results can be saved now */

	if "`cd'"!="" {
		global S_E_msg1 "(Standard errors scaled using `cd')"
	}
	if `reg' {
		global S_E_msg2 "(Model is ordinary regression, use regress instead)"
	}
	if `disp'!=1 {
		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
*/
	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'

	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  "`offset'"		/* name of offset/lnoffse var */
	global S_E_lno  = "`lnoffse'"!=""	/* 1 if lnoffset, 0 otherwise */
	global S_E_ef   "`eform'"
	global S_E_k    "`k'"  
	global S_E_disp "`disp'"
	global S_E_del  "`delta'"
	global S_E_cmd  "glmr"			/* must be last */

	_glmrpt, level(`level') `efopt'
end