*! brsq -- An improved R-squared                        STB-14: sg18
*! version 1.0.0        June 1993       Patrick Royston, Richard Goldstein
program define brsq
	version 3.0
	local varlist "req ex min(2)"
	local if "opt"
	local in "opt"
	local weight "aweight fweight" 
	local options "Boxcox(string) Tvars(string) Reweight Leave *"
	local small 1e-6
	parse "`*'"
	parse "`varlist'", parse(" ")
	if "`boxcox'"=="" {
		local boxcox 0
	}
	else {
		cap confirm var `boxcox'
		if _rc {
			confirm number `boxcox'
			local b `boxcox'
			local boxcox
			tempvar boxcox
			gen `boxcox' = `b'
		}
	}
	parse "`varlist'", parse(" ")
	tempvar y yfit s logwt
	local lhs "`1'"
	mac shift
	local rhs1 "`*'"
	if "`tvars'"=="" { local rhs2 "`rhs1'" }
	else { local rhs2 "`tvars'" }
	quietly {
		gen `y' = `lhs' `if' `in'
		gen `s' = 1 if `y'!=.
		local userwt "[aweight=`s']"
		local mnlnwt 0
		if "`exp'"!="" {
			replace `s' `exp' if `y'!=.
			replace `y' = . if `s'<=0 | `s'==.
			if "`weight'"=="fweight" { 
				capture assert `s'==int(`s') if `y'!=.
				if _rc { noisily error 401 } 
			}
			sum `s'
			replace `s' = `s'/_result(3)
/*
	Mean log normalized weights
*/
			gen `logwt' = log(`s')
			sum `logwt'
			local mnlnwt = _result(3)
			drop `logwt'
		}
/*
	Deviances and df, untransformed data.
	Preserve fitted values.
*/
		reg `y' `userwt'
		local obs = _result(1)
		local ssr0u = _result(4)
		local dev0u = `obs'*(1+log(2*_pi*`ssr0u'/`obs')-`mnlnwt')
		reg `y' `rhs1' `userwt',`options'
		if _result(5)<1 { error 2001 } /* need at least 1 residual df */
		pred `yfit'
		local dfmu = _result(3)
		local ssr = _result(4)
		local dfru = _result(5)
		local rsdu = _result(9)
		local devu = `obs'*(1+log(2*_pi*`ssr'/`obs')-`mnlnwt')
/*
	For transformation, find the geometric mean of y
	and transform y to standardized form.
	First get deviances for unstandardized transformation.
*/
		tempvar tempy
		gen `tempy' = log(`y')
		sum `tempy'
		local gmean = exp(_result(3))
		replace `tempy' = -2*(`boxcox'-1)*log(`y'/`gmean')
		sum `tempy'
		local j = _result(3)*`obs'
		replace `tempy' = cond(abs(`boxcox')<`small', /*
		 */ log(`y'), /*
		 */ (`y'^`boxcox'-1)/`boxcox')
		if "`reweigh'"!="" {
			replace `s' = `s'*exp((2-2*`boxcox')*log(abs(`yfit')))
			sum `s'
			replace `s' = `s'/_result(3)
			gen `logwt' = log(`s')
			sum `logwt'
			local mnlnwt = _result(3)
		}
		reg `tempy' `userwt'
		local ssr = _result(4)
		local dev00t = `obs'*(1+log(2*_pi*`ssr'/`obs')-`mnlnwt')
		replace `y' = cond(abs(`boxcox')<`small', /*
		 */ `gmean'*`tempy', /*
		 */ `tempy'/`gmean'^(`boxcox'-1))
		drop `tempy'
		if "`reweigh'"!="" {
			replace `s' = `s'*exp((2*`boxcox'-2)*log(`gmean'))
			sum `s'
			replace `s' = `s'/_result(3)
			replace `logwt' = log(`s')
			sum `logwt'
			local mnlnwt = _result(3)
			drop `logwt'
		}
/*
	Deviances
*/
		reg `y' `userwt'
		local ssr = _result(4)
		local dev0t = `obs'*(1+log(2*_pi*`ssr'/`obs')-`mnlnwt')+`j'
		reg `y' `rhs2' `userwt',`options'
		local dfmt = _result(3)
		local ssr = _result(4)
		local dfrt = _result(5)
/*
	For some reason, the following expression for the residual sd
	with transformed model is incorrect with reweighting.
	Rsd calculated instead via the scaled R^2.

		local rsdt = sqrt(`ssr'/_result(5))*exp(-`j'/`obs')
*/
		local adevt = `obs'*(1+log(2*_pi*`ssr'/`obs')-`mnlnwt')+`j'
		local adevu `devu'
		local prRsqu = 1-exp(-(`dev0u'-`adevu')/`obs')
		local aprRsqu = 1-((`obs'-1)/`dfru')*(1-`prRsqu')
		local prRsqt = 1-exp(-(`dev0u'-`adevt')/`obs')
		local aprRsqt = 1-((`obs'-1)/`dfrt')*(1-`prRsqt')
		local rsdt = sqrt((1-`prRsqt')*`ssr0u'/`dfrt')
		local dev0t = `dev00t'-(`dev0t'-`adevt')
		local fu = (`dfru'/`dfmu')*`prRsqu'/(1-`prRsqu')
		local ft = (`dfrt'/`dfmt')*`prRsqt'/(1-`prRsqt')
		if "`leave'"!="" {
			cap drop _fy
			cap drop _fyf
			cap drop _sfy
			gen _fy = `y'
			pred _fyf
			gen _sfy = `rsdt'/sqrt(`s')
		}
	}
	di in gr _n "An improved R-squared" _n _dup(22) "-"
	di in gr "Model: Y-var = " in ye "`lhs'" in gr /*
		*/ ", X-vars = " in ye "`rhs1'" _c
	if "`tvars'"!="" {
		di in gr ", T-vars = " in ye "`tvars'" _c
	}
	#delimit ;
	di in gr _n(2) _col(18) "Untransformed Y-var"
		_col(41) "Transformed Y-var" ;
	di in gr _col(18) _dup(19) "-" _col(41) _dup(17) "-" ;
	di in gr "Log-likelihood"
		_col(19) in ye %12.3f `devu'/-2 _col(40) %12.3f `dev0t'/-2 ;
	di in gr "Deviance-0"
		_col(19) in ye %12.3f `dev0u'   _col(40) %12.3f `dev00t' ;
	di in gr "Deviance"
		_col(19) in ye %12.3f `devu'    _col(40) %12.3f `dev0t' ;
	di in gr "Deviance-Scaled"
		_col(19) in ye %12.3f `adevu'   _col(40) %12.3f `adevt' ;
	di in gr _n "R-squared-Scaled"
		_col(19) in ye %12.4f `prRsqu'  _col(40) %12.4f `prRsqt' ;
	di in gr "Adj-R-sq-Scaled"
		_col(19) in ye %12.4f `aprRsqu' _col(40) %12.4f `aprRsqt' ;
	di in gr "Residual-SD-Scaled"
		_col(19) in ye %12.4f `rsdu'  _col(40) %12.4f `rsdt' ;
	#delimit cr
	di in gr _n "Note: " _c
	if "`reweigh'"=="" {
		di "same" in gr /*
		 */ " weights used in untransformed and transformed analyses"
	}
	else {
		di "transformed" in gr " weights used in transformed analysis"
	}
	mac def S_1 `adevu'
	mac def S_2 `adevt'
	mac def S_3 `prRsqu'
	mac def S_4 `prRsqt'
	mac def S_5 `aprRsqu'
	mac def S_6 `aprRsqt'
	mac def S_7 `rsdu'
	mac def S_8 `rsdt'
	mac def S_9 `fu'
	mac def S_10 `ft'
end