*! version 2.0.0 PR 27oct2014
program define ptrend, rclass /* chisq for trend */
	version 8.0
	syntax varlist(min=3 max=4 numeric) [if] [in]
	tokenize `varlist'
	local r `1'
	local nr `2'
	local x `3'
	local p `4'
	tempvar n pr touse
	quietly {
		marksample touse
		count if `touse'
		local rows = r(N)
		if (`rows' < 2) error 2001
		local r1 = `rows'-1
		gen long `n' = `r'+`nr' if `touse'
		if "`p'"=="" {
			local p "_prop"
			cap drop `p'
			gen `p' = `r'/`n'
			format `p' %8.3f
		}
/*
	Do regression for trend using binomial weights n for p = r/n.
	R = sum(r), N = sum(n) = sum(weights), P = R/N.

	Formula for chisq for trend adapted from that in TREND.MTB, needing
	some careful algebra and calculation of the weighted SSQ of x.
*/
		sum `r' if `touse'
		local R = r(mean)*`rows'
		sum `n'
		local N = r(mean)*`rows'
		local P = `R'/`N'
		sum `x' if `touse' [weight=`n']
		local ssx = `N'*r(Var)*(`r1')/`rows'
		regress `p' `x' [weight=`n']
		local b2 = _b[`x']^2
		local chitr = `b2'*`ssx'/(`P'*(1-`P'))
		local se = sqrt(`b2'/`chitr')
/*
	Ordinary chisquare
*/
		gen `pr' = sum(`p'*`r')
		local chisq = (`pr'[_N]-`R'*`P')/(`P'*(1-`P'))
		local chidep = `chisq' - `chitr'
		local dfdep = `rows'-2
	}
	list `r' `nr' `p' `x' if `touse'
	#delimit ;
	di as txt _n "Trend analysis for proportions" _n "{hline 30}" _n ;
	di as txt "Regression of p = " as res "`r'" as txt "/("
	 as res "`r'" as txt "+"
	 as res "`nr'" as txt ") on " as res "`x'" as txt ":" _n ;
	di as txt "Slope = " as res %7.0g _b[`x']
	 as txt ", std. error = " as res %7.0g `se'
	 as txt ", Z = " as res %7.3f sqrt(`chitr') ;
	di as txt _n "Overall chi2(" as res `rows'-1 as txt ") = " _col(27)
	 %7.3f as res `chisq' as txt ",  pr>chi2 = "
	 %6.4f as res chiprob(`r1', `chisq');
	di as txt "Chi2(" as res 1 as txt ") for trend = " _col(27)
	 %7.3f as res `chitr' as txt ",  pr>chi2 = "
	 %6.4f as res chiprob(1, `chitr') ;
	di as txt "Chi2(" as res `dfdep' as txt ") for departure = " _col(27)
	 %7.3f as res `chidep' as txt ",  pr>chi2 = "
	 %6.4f as res chiprob(`dfdep', `chidep') ;
	#delimit cr
	return scalar slope = _b[`x']
	return scalar se = `se'
	return scalar chi2trend = `chitr'
	return scalar chi2overall = `chisq'
	return scalar chi2dep = `chidep'
	return scalar dfdep = `dfdep'
	return scalar Pdep = chiprob(`dfdep', `chidep')
end