*! version 1.0.5  16jun1999 STB-50 sbe29
program define binreg, eclass
	version 6.0
	if replay() {
		Playback `0'
		exit
	}

	syntax varlist [if] [in] [aw fw iw] [, LEvel(int $S_level) N(string) /*
		*/ LTol(real -1) ITerate(int 0) DISP(real 1) Scale(string) /*
		*/ FRVars Offset(varname numeric) LNOffset(varname numeric) /*
		*/ noCONStant INIt(varname numeric) noLOg OR RR RD HR COEFF ]

	local nc : word count `or' `rr' `rd' `hr'
	if `nc' > 1 {
		di in red "May choose only one of OR, RR, RD, HR"
		exit 198
	}
	local eopt 0 
	if `nc' != 0 & "`coeff'" == "" { local eopt 1 } 
	if `nc' == 0 { 
		local link "l" 
		local or "or"
	}

	if "`rd'" != "" {
		if "`coeff'" == "" { local eopt 0 }
		else               { local eopt 1 } 
	}

	marksample touse 		/* must be done early */

	if "`n'" == ""			{ local m 1		}
	else				{ local m `"`n'"'	}
	if `"`log'"'!=""		{ local iflog `"*"'	}
	else				{ local iflog "noi"	}
	if `"`constan'"'!=""		{ local cons "nocons"	}
	if `level'<10 | `level'>99	{ local level 95	}
	if `ltol'<0			{ local ltol 1e-6	} 
	if `iterate'<=0			{ local iterate 50	}

	local small 1e-6

	local family	`"bin"'
	local pow	1
	local k		0
	local scale1	1

	capture confirm integer number `m'
	local mfixed = (_rc==0)		/* Fixed if not a variable */

	if "`or'" != ""		{ 
		local link "l"				/* logistic	*/
		local efopt "eform(Odds Ratio)"
		local cmsg  "Odds ratio coefficients"
	}
	else if "`rr'" != ""	{ 
		local link "g"				/* log		*/
		local efopt "eform(Risk Ratio)"
		local cmsg  "Risk ratio coefficients"
	}
	else if "`rd'" != ""	{ 
		local link "i"				/* identity	*/
		/* local efopt "eform(Risk Diff.)" */
		local cmsg  "Risk difference coefficients"
	}
	else if "`hr'" != ""	{ 
		local link "c"				/* log compl.	*/
		local efopt "eform(Health Ratio)"
		local cmsg  "Health ratio coefficients"
	}

	if !`eopt' { local efopt } 

	local bernoul = `mfixed' & `"`m'"'=="1"		/* Bernoulli dist */

/*
	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
			}
		}
	}

 	tempvar mu eta W z V dev wt
	quietly {
		tokenize `varlist'
		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
			local offstr `"ln(`lnoffse')"'
			gen double `offvar' = `offstr'
		}
		else if `"`offset'"'!="" {
			local offvar `"`offset'"'
			local offstr `"`offset'"'
		}
		if `"`offstr'"'!="" {
			local offopt `", offset `offstr'"'
		}

		markout `touse' `offvar' `mmiss'

		sum `y' if `touse'	/* purposely not weighted */
		if r(N)==0 { noisily error 2000 } 
		if r(N)==1 { noisily error 2001 }
		if r(min)==r(max) & `"`offstr'"'=="" {
			noi di _n in bl "outcome does not vary"
			exit 2002
		}

		local nobs = r(N) /* nobs reset below when fweights */
		if `"`weight'"'!="" { 
			gen double `wt' `exp' if `touse'
			if `"`weight'"'=="aweight" {
				summ `wt'
				replace `wt' = `wt'/r(mean)
			}
			else if `"`weight'"'=="fweight" {
				summ `touse' [fw=`wt'] if `touse'
				local nobs = round(r(N),1)
			}
		}
		else gen byte `wt' = `touse' if `touse'

/* check for sensibility of dependent variable */

		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 */
		if `"`init'"'!="" {
			gen double `mu' = `init' if `touse'
		}
		else {
			sum `y' [aw=`wt'] if `touse'
			gen double `mu' = (`y'+r(mean))/(`m'+1)
		}
		gen double `W' = . in 1
		gen double `z' = . in 1
		gen double `V' = . in 1
		gen double `dev' = . in 1

/* Initialize linear predictor */
		if `"`init'"'=="" {
			replace `mu' = `m'*(`y'+.5)/(`m'+1) /*
			 */ if `touse'
		}

		if `"`link'"'=="l" {
			gen double `eta' = ln(`mu'/(`m'-`mu'))
		}
		else if `"`link'"'=="g" {
			gen double `eta' = ln(`mu'/`m')
		}
		else if `"`link'"'=="c" {
			gen double `eta' = ln(1-`mu'/`m')
		}
		else {
			gen double `eta' = `mu'/`m'
		}

/* Local scoring algorithm: IRLS */
		tempname newdev oldev Wscale
		scalar `oldev' =  0
		scalar `newdev' =  1  

		local i 1
		`iflog' di
		while abs(`newdev'-`oldev')>`ltol' & `i'<=`iterate' {
			scalar `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.
*/
			CalcV `touse' `y' `family' `"`link'"' `pow' `m' /*
			 */ `disp' `k' `eta' `mu' `W' `z' `V' `"`offvar'"'

/* Weighted regression `Wscale' added by wms 10 Apr 1995 */
                        summarize `W' [aw=`wt']
                        scalar `Wscale' = r(mean)
                        regress `z' `*' [iw=`W'*`wt'/`Wscale'], mse1 `cons'

/* Get linear predictor */
			cap drop `eta'
			_predict double `eta' if `touse' 
			if `"`offstr'"'!="" { replace `eta' = `eta'+`offvar' }

/* Ensure predictions are within range and get inverse link */
			if "`link'" == "l" {
				replace `mu' = 1/(1+exp(-`eta'))
			}
			else if "`link'" == "g" {
				replace `eta' = min(-.0001,`eta')
				replace `mu' = exp(`eta')
			}
			else if "`link'" == "c" {
				replace `eta' = min(-.0001,`eta')
				replace `mu' = 1-exp(`eta')
			}
			else {
				replace `eta' = max(.0001,`eta')
				replace `eta' = min(.9999,`eta')
				replace `mu' = `eta'
			}
			if !`bernoul' {
				replace `mu' = `mu'*`m'
			}

			count if `mu'==. & `eta' != . & `touse'
			if r(N) != 0 {
				noi di in red "estimates diverging"
				exit 430
			}

/* Get weighted squared deviance contributions (residuals) */
			_crcgldv `y' `family' `bernoul' `m' `k' `mu' `wt' `dev'

/* Get deviance (note: analytic weights are already normalized). */
			replace `z' = sum(`dev')
			scalar `newdev' = `z'[_N]/`disp'
			`iflog' di in gr `"Iteration `i' : deviance = "' /*
			 */ in ye %9.4f `newdev'
			local i = `i'+1
		}
		local conrc = cond(abs(`newdev'-`oldev')>`ltol',430,0)
		if `conrc' { 
			noisily di in ye "(convergence not achieved)"
		}

/* Get weighted squared contributions (residuals) to Pearson X2 */
		replace `z' = sum(`wt'*(`y'-`mu')^2/`V')
		local chisq = `z'[_N]/`disp'
		local df = `nobs'-e(df_m)-(`"`cons'"'=="")
		local dispc = `chisq'/`df'
		local dispd = `newdev'/`df'
	}

/* Store results */
	DescFL `"`family'"' `"`link'"' `pow' `bernoul' `k' `m'
	local o `"`r(dist)', `r(link)'`offopt'"'
	local lll : word 1 of `r(link)'

	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 scalar `Wscale' = `Wscale'/`delta' /* scale to `delta' */
	} 

	tempname b V
	mat `b' = get(_b)
	mat `V' = get(VCE)
        scalar `Wscale' = 1/`Wscale'
        mat `V' = `Wscale'*`V' /* get rid of `Wscale' scaling */
	if `"`zapse'"'=="yes" { 
		local i 1
		while `i'<=rowsof(`V') {
			mat `V'[`i',`i'] = 0 
			local i=`i'+1
		}
	}
	tempvar mysamp
	qui gen byte `mysamp' = `touse'
	est post `b' `V', depname(`y') obs(`nobs') `dofopt' esample(`mysamp')

/* results can be saved now */

	if `"`cd'"'!="" {
		est local msg_1 `"(Standard errors scaled using `cd')"'
	}
	if `disp'!=1 {
		est local msg_2 /*
			*/ `"(Quasi-likelihood model with dispersion `disp')"'
	}
	if !`eopt' {
		est local msg_3 `"`cmsg'"'
	}
	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 */
	est scalar df_pear = `df'
	est scalar N = `nobs'
	est scalar chi2 = `chisq'
	est scalar deviance = `newdev'
	est scalar dispersp = `dispc'
	est scalar dispers = `dispd'
	est local  wtype `"`weight'"'
	est local  wexp `"`exp'"'
	est local  family `"`family'"'
	est local  link `"`lll'"'
	est local  title_fl `"`o'"'
	est scalar bernoul = `bernoul'
	est local  m `"`m'"'
	est local  depvar `"`y'"'
	est scalar power = `pow'
	est local  offset `"`offstr'"'	        /* offset/lnoffset */
	est local  eform `"`eform'"'
	est scalar k = `k'  
	est scalar disp = `disp'
	est scalar delta = `delta'
	est scalar rc = `conrc'
	est local  predict "binr_p"
	est local  cmd "binreg"

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


program define Playback
	syntax [, LEvel(int $S_level) EForm(passthru)]
	if `"`e(cmd)'"'!="binreg" { error 301 } 

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

	#delimit ;
	di _n in gr "Residual df  = " in ye %9.0g e(df_pear)
		_col(57) in gr "No. of obs = "  in ye %9.0g e(N) ;
	di in gr "Pearson X2   = " in ye %9.0g e(chi2) in gr 
		_col(57) "Deviance   = "  in ye %9.0g e(deviance) ;
	di in gr "Dispersion   = " in ye %9.0g e(dispersp) in gr
		_col(57) "Dispersion = "  in ye %9.0g e(dispers) ;
	#delimit cr

	di in gr _n `"`e(title_fl)'"' /* family, link, offset */
	if `"`e(msg_3)'"' != "" {
		di in gr "`e(msg_3)'"
	}
	est di, level(`level') `eform'
	local i 1
	while `i' <= 2 {
		if `"`e(msg_`i')'"' != "" { di in gr `"`e(msg_`i')'"' }
		local i = `i' + 1
	}
	error `e(rc)'
end


program define DescFL, rclass /* return description of family and link */
	args f l pow bernoul k m

	local s1 " distribution"
	local s2 " link"

	if `bernoul'	{ local s1 `"Bernoulli`s1'"' }
	else 		{ local s1 `"Binomial (N=`m')`s1'"' }

	if `"`l'"'=="l"		{ local s2 `"logit`s2'"' }
	else if `"`l'"'=="g"	{ local s2 `"log`s2'"' }
	else if `"`l'"'=="c"	{ local s2 `"log-complement`s2'"' }
	else if `"`l'"'=="i"	{ local s2 `"identity`s2'"' }

	ret local dist `"`s1'"'
	ret local link `"`s2'"'
end


* Update of _glmwgt (sg16.3), PR 22-Dec-93, 11-Apr-94, 11-Oct-94.
* Calculates var functn, iterative weights and working dependent variable.
program define CalcV
	args touse y fam link pow m disp k eta mu W z V moffset
	if `"`moffset'"'!="" {
		local moffset `"-`14'"'
	}

/* Calculate variance function, V, for each family */
	tempvar w
	replace `V' = `mu'*(1-`mu'/`m')

	if `"`link'"' == "l" {
		gen double `w' = `V'
	}
	else if `"`link'"' == "i" {
		gen double `w' = `m'    /*   `m'==N in paper */
	}
	else if `"`link'"' == "g" {	
		gen double `w' = `mu'
	}
	else if `"`link'"' == "c" {	
		gen double `w' = `m'*(`mu'/`m'-1)
	}

	replace `W' = `disp'*`w'^2/`V' /* `disp' could be in wrong place */
	replace `z' = `eta'+(`y'-`mu')/`w' `moffset' 
end
exit


Notes:
	This program fits the model described in:

	Wacholoder, Sholom 1986.
		Binomial regression in GLIM: estimating 
		risk ratios and risk differences
		American Journal of Epidemiology, 123:174--184.






James Hardin and Mario Cleves
Email correspondence to: mcleves@@stata.com