*!version 1.0.2 05/04/05 (STB-56: sg149) /* Mark S. Pearce, University of Newcastle upon Tyne */ /* Tests for seasonality with a variable population at risk */ /* Follows the tests given in Walter & Elwood (1975), Edwards (1961) */ /* Default is Walter & Elwood test using equal month lengths */ /* 05/04/05 Corrections to fitted (expected) frequencies and angle of maximum rate */ /* We have become aware of a problem when the program calculates the fitted, or expected, frequencies. In some instances, these frequencies are actually flipped around the horizontal axis, simply because the calculation of the angle has also been inverted. The angle flips is because tan(theta) = tan(theta+pi) = tan(theta-pi). Therefore, when the atan function is calculated, it picks an angle restricted to +- pi/2 and cannot produce any angles outside this range. In other words, no maximum rates could be predicted to lie outside October to March. A correction has been made by adding pi to thstar if xbar-xhat <0 (i.e. when the peak is April to September). We are grateful to Darren Greenwood for helpful discussions with this problem.*/ program define seast version 6.0 local options "EDwards EXAct NOTab Generate(string) SECtor(string) LENgth(string)" local varlist "req ex min(1) max(2)" parse "`*'" parse "`varlist'", parse(" ") local obs="`1'" if "`2'"=="" { tempvar pop local nopop=1 gen `pop'=1 } else { local pop ="`2'" local nopop=0 } if "`generate'"!="" { confirm new var `generate' } if "`sector'"=="" { local sector="month" } tempvar theta weight xmean ymean srtm cosbit sinbit sumcos sumsin nmi totb totc sumnmi ci days x ex1 ey1 ex2 ey2 sumci exp ti sumti oe2 gfit suml rate rank quietly { summ `sector' local k=r(max) egen `totb'=sum(`pop') summ `totb' local M=r(max) egen `totc'=sum(`obs') summ `totc' local N=r(max) local pi=3.1415927 if "`length'"!="" { egen `suml'=sum(`length') summ `suml' gen `x'=(360/`suml')*`length'*(1/57.296) gen `theta'=`pi'/`k' local stheta=`pi'/`k' local i=2 while `i'<=`k' { summ `x' if `sector'==`i'-1 replace `theta'=`stheta'+r(max) if `sector'==`i' local stheta=`stheta'+r(max) local i=`i'+1 } } else { if "`exact'"=="" { gen `theta' =`pi'/`k'*((2*`sector')-1) } else { if `k'==12 { gen `days'=28+2/9 if `sector'==2 for num 1 3 5 7 8 10 12: replace `days'=31 if `sector'==X replace `days'=30 if `days'==. gen `x'=(360/365.25)*`days'*(1/57.296) gen `theta'=`pi'/12 local stheta=`pi'/12 local i=2 while `i'<=12 { summ `x' if `sector'==`i'-1 replace `theta'=`stheta'+r(max) if `sector'==`i' local stheta=`stheta'+r(max) local i=`i'+1 } } else { if "`length'"=="" { di in red "Can only use exact with monthly data or by using the length option" e 100 } } } } egen `weight'=sum(sqrt(`obs')) summ `weight' local W=r(max) *Calcluate observed centre point egen `xmean'=sum((sqrt(`obs'))*(cos(`theta'))) summ `xmean' local xbar=r(max)/`W' egen `ymean'=sum((sqrt(`obs'))*(sin(`theta'))) summ `ymean' local ybar=r(max)/`W' *Getting expected centre point if "`edwards'"=="" { egen `srtm' = sum(sqrt(`pop')) summ `srtm' local srmi=r(max) gen `ex1'= (sqrt(`pop'))*(cos(`theta')) egen `ex2'=sum(`ex1') gen `ey1' =(sqrt(`pop'))*(sin(`theta')) egen `ey2'=sum(`ey1') summ `ex2' local xhat=r(max)/`srmi' summ `ey2' local yhat=r(max)/`srmi' } else { local xhat=0 local yhat =0 } *Variances gen `cosbit'=0.25*(cos(`theta'))*(cos(`theta')) gen `sinbit'=0.25*(sin(`theta'))*(sin(`theta')) egen `sumcos'=sum(`cosbit') summ `sumcos' local vx1=r(max) egen `sumsin'=sum(`sinbit') summ `sumsin' local vy1=r(max) if "`edwards'"=="" { gen `nmi'=sqrt(`N'*`pop'/`M') egen `sumnmi'=sum(`nmi') summ `sumnmi' local v2=r(max)*r(max) local varx= `vx1'/`v2' local vary=`vy1'/`v2' } local dsq=((`xbar'-`xhat')^2)+((`ybar'-`yhat')^2) local d=sqrt(`dsq') local thstar=atan((`ybar'-`yhat')/(`xbar'-`xhat')) if (`xbar'-`xhat'<0) { local thstar=_pi+`thstar' } local alpha = 4*`d' if "`edwards'"!="" { gen `ci'=(1+`alpha'*(cos(`theta'-`thstar'))) egen `sumci'=sum(`ci') summ `sumci' gen `exp'=`N'*`ci'/r(max) local test = 0.5*`alpha'*`alpha'*`N' } else { gen `ci'=`pop'*(1+`alpha'*(cos(`theta'-`thstar'))) egen `sumci'=sum(`ci') summ `sumci' gen `exp'=`N'*`ci'/r(max) local test = (((`xbar'-`xhat')/sqrt(`varx'))^2) + (((`ybar'-`yhat')/sqrt(`vary'))^2) } gen `oe2'=((`obs'-`exp')^2)/`exp' egen `gfit'=sum(`oe2') summ `gfit' local fit=r(max) } *end of quietly local thstar = `thstar'*57.296 if "`notab'"=="" { di in green _dup(50) "-" if `k'==12 { di in green "Month" _col(20) "|" _col(22) "Observed" _col(35) "Expected" } else { di in green "Sector" _col(20) "|" _col(22) "Observed" _col(35) "Expected" } di in green _dup(19) "-" _col(20) "+" _dup(30) "-" local i=1 while `i'<=`k' { quietly summ `obs' if `sector'==`i' local obsi=r(max) quietly summ `exp' if `sector'==`i' local expi=r(max) if `k'!=12 { di in gr `i' _col(20) "|" in ye _col(22) `obsi' _col(35) `expi' } else { if `i'==1 { di in gr "January" _col(20) "|" in ye _col(22) `obsi' _col(35) `expi' } else if `i'==2 { di in gr "February" _col(20) "|" in ye _col(22) `obsi' _col(35) `expi' } else if `i'==3 { di in gr "March" _col(20) "|" in ye _col(22) `obsi' _col(35) `expi' } else if `i'==4 { di in gr "April" _col(20) "|" in ye _col(22) `obsi' _col(35) `expi' } else if `i'==5 { di in gr "May" _col(20) "|" in ye _col(22) `obsi' _col(35) `expi' } else if `i'==6 { di in gr "June" _col(20) "|" in ye _col(22) `obsi' _col(35) `expi' } else if `i'==7 { di in gr "July" _col(20) "|" in ye _col(22) `obsi' _col(35) `expi' } else if `i'==8 { di in gr "August" _col(20) "|" in ye _col(22) `obsi' _col(35) `expi' } else if `i'==9 { di in gr "September" _col(20) "|" in ye _col(22) `obsi' _col(35) `expi' } else if `i'==10 { di in gr "October" _col(20) "|" in ye _col(22) `obsi' _col(35) `expi' } else if `i'==11 { di in gr "November" _col(20) "|" in ye _col(22) `obsi' _col(35) `expi' } else if `i'==12 { di in gr "December" _col(20) "|" in ye _col(22) `obsi' _col(35) `expi' } } local i=`i'+1 } di in green _dup(50) "-" } di _n if "`edwards'"!="" { di in gr "Edwards Test" di _n } else { if "`exact'"!="" { di in green "Walter & Elwood Test (using exact month lengths)" } else { di in gr "Walter & Elwood Test" } di _n } quietly { tempvar sumobs gen `sumobs'=sum(`obs') summ `sumobs' local totobs=r(max) } di in green "Total Number of Cases = " in yel `totobs' di _n di in green "Seasonality Test" _col(50) "Goodness of fit test" di in green "chi2(" in yel "2" in gr ")" _col(15) "= " in yel %9.4f `test' _col(50) in gr "chi2(" in yel "2" in gr ")" _col(65) "=" in yel %9.4f `fit' di in green "Prob > chi2" _col(15) "= " in yel %9.4f chiprob(2,`test') _col(50) in gr "Prob > chi2" _col(65) "=" in yel %9.4f chiprob(`k'-1,`fit') di _n di in green _dup(50) "-" di in green "Parameter" _col(38) "|"_col(40) "Estimate" di in green _dup(37) "-" _col(38) "+" _dup(12) "-" di in green "Amplitude of cyclic variation" _col(38) "|" _col(40) in yellow `alpha' di in green "Angle of maximum rate" _col(38) "|" _col(40) in yellow `thstar' di in green _dup(50) "-" if "`generate'"!="" { gen `generate'=`exp' } end