/*** This program provides the data and analysis for the main model reported in: ***/ /*** Walter et al. (2011) Accident Analysis and Prevention 43: 2064-2071. ***/ /*** These data are reproduced by permission of the NSW Ministry of Health ***/ /*** Custom Formats ***/ proc format; value injury 0 = 'Arm' 1 = 'Head'; value law 0 = 'Pre' 1 = 'Post'; run; /*** Data ***/ *Seasonally adjusted injury counts; *Note: seasonal adjustment via PROC X11 gives non-integer counts; data mhl; input time law injury count population cluster; LogPop=log(population); format injury injury. law law.; datalines; -17.5 0 0 64.47487894 5781579.792 1 -16.5 0 0 64.77414938 5786876.584 2 -15.5 0 0 59.00485663 5792173.376 3 -14.5 0 0 67.38284827 5797470.168 4 -13.5 0 0 75.28160988 5802766.96 5 -12.5 0 0 58.11457042 5808063.752 6 -11.5 0 0 73.7252308 5809961.67 7 -10.5 0 0 55.32709198 5814772.909 8 -9.5 0 0 62.85498662 5819584.147 9 -8.5 0 0 70.6898182 5824395.386 10 -7.5 0 0 78.88780608 5829206.624 11 -6.5 0 0 55.29506437 5834017.863 12 -5.5 0 0 55.61051436 5838832.761 13 -4.5 0 0 55.10430419 5843644.523 14 -3.5 0 0 59.76705742 5848456.284 15 -2.5 0 0 66.1524221 5853268.046 16 -1.5 0 0 63.84205226 5858079.807 17 -0.5 0 0 61.19292606 5862891.569 18 0.5 1 0 46.08706292 5870484.281 19 1.5 1 0 84.99399453 5875693.321 20 2.5 1 0 43.16521244 5880902.362 21 3.5 1 0 44.58563132 5886111.402 22 4.5 1 0 57.18844033 5891320.442 23 5.5 1 0 65.33160457 5896529.482 24 6.5 1 0 59.58535777 5904306.96 25 7.5 1 0 62.65281244 5909882.92 26 8.5 1 0 58.88563162 5915458.88 27 9.5 1 0 48.95443806 5921034.839 28 10.5 1 0 55.05635145 5926610.799 29 11.5 1 0 62.52178732 5932186.759 30 12.5 1 0 60.20458046 5938525.238 31 13.5 1 0 68.04549506 5944210.129 32 14.5 1 0 79.9502655 5949895.021 33 15.5 1 0 68.91860007 5955579.912 34 16.5 1 0 60.62523487 5961264.803 35 17.5 1 0 58.81606924 5966949.694 36 -17.5 0 1 62.50293693 5781579.792 1 -16.5 0 1 71.97713927 5786876.584 2 -15.5 0 1 66.96582169 5792173.376 3 -14.5 0 1 81.7965882 5797470.168 4 -13.5 0 1 96.57650578 5802766.96 5 -12.5 0 1 63.79850571 5808063.752 6 -11.5 0 1 70.74545503 5809961.67 7 -10.5 0 1 72.23872835 5814772.909 8 -9.5 0 1 73.65161838 5819584.147 9 -8.5 0 1 80.89882985 5824395.386 10 -7.5 0 1 84.17992221 5829206.624 11 -6.5 0 1 61.34604578 5834017.863 12 -5.5 0 1 67.72691578 5838832.761 13 -4.5 0 1 41.53723028 5843644.523 14 -3.5 0 1 65.24200163 5848456.284 15 -2.5 0 1 71.94773142 5853268.046 16 -1.5 0 1 58.27290466 5858079.807 17 -0.5 0 1 78.61231212 5862891.569 18 0.5 1 1 37.04651811 5870484.281 19 1.5 1 1 61.83736123 5875693.321 20 2.5 1 1 37.93629673 5880902.362 21 3.5 1 1 43.83379648 5886111.402 22 4.5 1 1 42.20046243 5891320.442 23 5.5 1 1 43.65799765 5896529.482 24 6.5 1 1 44.15371916 5904306.96 25 7.5 1 1 55.23137179 5909882.92 26 8.5 1 1 46.21871206 5915458.88 27 9.5 1 1 44.52648229 5921034.839 28 10.5 1 1 40.76181681 5926610.799 29 11.5 1 1 52.47785469 5932186.759 30 12.5 1 1 72.77560565 5938525.238 31 13.5 1 1 57.26175012 5944210.129 32 14.5 1 1 65.9085301 5949895.021 33 15.5 1 1 50.71240368 5955579.912 34 16.5 1 1 38.09350906 5961264.803 35 17.5 1 1 61.99779471 5966949.694 36 ; /*** Models ***/ *Negative binomial model to produce estimates reported by Walter et al. (2011); proc genmod data=mhl; class injury(ref=first) law(ref=first)/param=ref order=internal; model count=time|injury|law / dist=nb link=log offset=LogPop ; output out=res resraw=raw reschi=pearson; run; *Negative binomial GEE model to adjust for within-month correlation between head and arm injury counts; proc genmod data=mhl; class injury(ref=first) law(ref=first) cluster/param=ref order=internal; model count=time|injury|law / dist=nb link=log offset=LogPop ; repeated subject=cluster; run; /*** Residual plots ***/ proc gplot data=res; by injury; plot raw*time; run; /*** Autocorrelation and cross-correltion ***/ data head(keep=time count law pearson raw) arm(keep=time count law pearson raw); set res; if injury=1 then output head; if injury=0 then output arm; run; data arm; set arm; rename count=count2 pearson=pearson2 raw=raw2; run; data res; merge head arm; by time; run; *Autocorrelation for arm injuries; proc arima data=res; identify var=raw2; run; *Autocorrelation for head injuries and cross-correlation between head and arm injuries; proc arima data=res; identify var=raw crosscorr=raw2; run; /*** End of program ***/