*! 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