/*---------------------------------------------- | Data from student in SAS forecasting | | class in Toronto. Hospital admissions in the | | years 2002 and 2003(SARS year) with period | | being pre-SARS (1), height of epidemic (2), | | post epidemic (3) each year. These are for | | something that could be postponed (I think | | it is asthma related visits) with the point | | being to see the effect of SARS on admissions | *----------------------------------------------*/ options symbolgen; %let days = d1 d2 d3 d4 d5 d6; %let start = "31mar2003"d; Title "Hospital elective visits"; title2 "Assigning &start as start of SARS scare"; Data SARS; array D(7); input DATE date9. VISITS ADMITS period; format date date7.; * Create some useful variables for analysis; dow=weekday(date); *day of the week; do i=1 to 7; D(i)=(dow=i); end; *day indicators ; julday = put(date, julday3.)+0.000; *day of the year; year=year(date); *year (2 or 3); yr_p = year + period/10; *year.period; yr_dow =year+ dow/10; *year.day_of_week; SARS=(date=&start); *SARS outbreak; XX = put(&start,julday3.)+0.000; call symput("SARS",xx); perm =(date>&start-1); *permanent effect; d_slope = (julday-xx)*perm; *slope change; datalines; 01MAR2002 178 11 1 02MAR2002 171 13 1 03MAR2002 176 25 1 04MAR2002 150 10 1 05MAR2002 125 12 1 06MAR2002 134 13 1 07MAR2002 150 13 1 08MAR2002 146 14 1 09MAR2002 157 11 1 10MAR2002 157 18 1 11MAR2002 151 19 1 12MAR2002 139 15 1 13MAR2002 159 19 1 14MAR2002 127 18 1 15MAR2002 159 14 1 16MAR2002 143 12 1 17MAR2002 175 13 1 18MAR2002 128 22 1 19MAR2002 142 15 1 20MAR2002 145 24 1 21MAR2002 115 4 1 22MAR2002 135 17 1 23MAR2002 146 16 1 24MAR2002 166 14 1 25MAR2002 147 12 1 26MAR2002 127 13 1 27MAR2002 128 9 1 01MAR2003 153 7 1 02MAR2003 169 21 1 03MAR2003 157 20 1 04MAR2003 141 18 1 05MAR2003 125 19 1 06MAR2003 127 11 1 07MAR2003 145 19 1 08MAR2003 143 13 1 09MAR2003 141 18 1 10MAR2003 125 20 1 11MAR2003 115 16 1 12MAR2003 117 12 1 13MAR2003 129 24 1 14MAR2003 137 15 1 15MAR2003 134 13 1 16MAR2003 158 8 1 17MAR2003 148 12 1 18MAR2003 160 12 1 19MAR2003 141 18 1 20MAR2003 119 17 1 21MAR2003 136 21 1 22MAR2003 158 18 1 23MAR2003 149 10 1 24MAR2003 144 19 1 25MAR2003 137 15 1 26MAR2003 130 13 1 27MAR2003 124 16 1 28MAR2002 138 13 2 29MAR2002 166 9 2 30MAR2002 182 15 2 31MAR2002 160 7 2 01APR2002 153 16 2 02APR2002 126 18 2 03APR2002 141 16 2 04APR2002 127 10 2 05APR2002 137 15 2 06APR2002 136 18 2 07APR2002 157 11 2 08APR2002 141 18 2 09APR2002 155 16 2 10APR2002 135 14 2 11APR2002 130 16 2 12APR2002 147 17 2 13APR2002 157 11 2 14APR2002 142 19 2 15APR2002 157 20 2 16APR2002 137 23 2 17APR2002 128 12 2 18APR2002 138 12 2 19APR2002 129 16 2 20APR2002 147 10 2 21APR2002 167 17 2 22APR2002 129 14 2 23APR2002 120 15 2 24APR2002 134 13 2 25APR2002 106 11 2 26APR2002 127 14 2 27APR2002 119 17 2 28APR2002 131 14 2 29APR2002 134 9 2 30APR2002 120 13 2 01MAY2002 113 17 2 02MAY2002 112 11 2 03MAY2002 142 16 2 04MAY2002 143 14 2 28MAR2003 118 6 2 29MAR2003 132 14 2 30MAR2003 130 11 2 31MAR2003 87 14 2 01APR2003 69 16 2 02APR2003 64 11 2 03APR2003 64 12 2 04APR2003 55 9 2 05APR2003 72 13 2 06APR2003 83 15 2 07APR2003 70 15 2 08APR2003 65 12 2 09APR2003 63 13 2 10APR2003 59 11 2 11APR2003 69 13 2 12APR2003 94 14 2 13APR2003 98 14 2 14APR2003 64 13 2 15APR2003 60 15 2 16APR2003 80 12 2 17APR2003 78 11 2 18APR2003 77 12 2 19APR2003 86 14 2 20APR2003 94 11 2 21APR2003 91 7 2 22APR2003 78 17 2 23APR2003 66 16 2 24APR2003 68 9 2 25APR2003 76 11 2 26APR2003 89 17 2 27APR2003 74 6 2 28APR2003 118 6 2 29APR2003 132 14 2 30APR2003 130 11 2 01MAY2003 89 11 2 02MAY2003 79 12 2 03MAY2003 114 14 2 04MAY2003 87 8 2 05MAY2002 150 11 3 06MAY2002 140 13 3 07MAY2002 125 19 3 08MAY2002 127 15 3 09MAY2002 115 12 3 10MAY2002 156 16 3 11MAY2002 131 12 3 12MAY2002 156 9 3 13MAY2002 130 16 3 14MAY2002 106 9 3 15MAY2002 118 14 3 16MAY2002 116 15 3 17MAY2002 135 18 3 18MAY2002 117 19 3 19MAY2002 147 15 3 20MAY2002 152 14 3 21MAY2002 129 23 3 22MAY2002 115 13 3 23MAY2002 128 20 3 24MAY2002 118 21 3 25MAY2002 123 9 3 26MAY2002 130 11 3 27MAY2002 119 12 3 28MAY2002 130 17 3 29MAY2002 143 18 3 30MAY2002 136 17 3 31MAY2002 131 17 3 01JUN2002 143 17 3 02JUN2002 139 19 3 03JUN2002 139 20 3 04JUN2002 125 17 3 05JUN2002 107 10 3 06JUN2002 107 10 3 07JUN2002 125 12 3 05MAY2003 106 11 3 06MAY2003 91 17 3 07MAY2003 90 14 3 08MAY2003 72 18 3 09MAY2003 90 7 3 10MAY2003 94 13 3 11MAY2003 99 12 3 12MAY2003 101 16 3 13MAY2003 84 12 3 14MAY2003 100 15 3 15MAY2003 98 20 3 16MAY2003 80 7 3 17MAY2003 75 11 3 18MAY2003 125 17 3 19MAY2003 119 20 3 20MAY2003 127 11 3 21MAY2003 85 14 3 22MAY2003 119 17 3 23MAY2003 108 16 3 24MAY2003 107 19 3 25MAY2003 124 12 3 26MAY2003 102 14 3 27MAY2003 88 11 3 28MAY2003 109 15 3 29MAY2003 110 9 3 30MAY2003 95 15 3 31MAY2003 113 5 3 01JUN2003 118 18 3 02JUN2003 110 19 3 03JUN2003 110 11 3 04JUN2003 87 14 3 05JUN2003 89 7 3 06JUN2003 93 12 3 07JUN2003 112 6 3 07JUN2002 . . 3 08JUN2002 . . 3 09JUN2002 . . 3 10JUN2002 . . 3 11JUN2002 . . 3 12JUN2002 . . 3 13JUN2002 . . 3 14JUN2002 . . 3 15JUN2002 . . 3 16JUN2002 . . 3 17JUN2002 . . 3 18JUN2002 . . 3 19JUN2002 . . 3 20JUN2002 . . 3 21JUN2002 . . 3 22JUN2002 . . 3 23JUN2002 . . 3 24JUN2002 . . 3 25JUN2002 . . 3 26JUN2002 . . 3 27JUN2002 . . 3 28JUN2002 . . 3 29JUN2002 . . 3 30JUN2002 . . 3 01JUL2002 . . 3 07JUN2003 . . 3 08JUN2003 . . 3 09JUN2003 . . 3 10JUN2003 . . 3 11JUN2003 . . 3 12JUN2003 . . 3 13JUN2003 . . 3 14JUN2003 . . 3 15JUN2003 . . 3 16JUN2003 . . 3 17JUN2003 . . 3 18JUN2003 . . 3 19JUN2003 . . 3 20JUN2003 . . 3 21JUN2003 . . 3 22JUN2003 . . 3 23JUN2003 . . 3 24JUN2003 . . 3 25JUN2003 . . 3 26JUN2003 . . 3 27JUN2003 . . 3 28JUN2003 . . 3 29JUN2003 . . 3 30JUN2003 . . 3 01JUL2003 . . 3 ; proc sort data=SARS; by date; * Create the trend function for graphing; data SARS; set SARS; EXP = SARS + .98084*EXP; retain EXP; if EXP=. then EXP=0; Yhat3 = 162.48914 -0.19969*julday -74.65255*EXP +7.87668*D1+0.12568*D2 -8.86646*D3-10.70244*D4-15.18572*D5-7.23870*D6; Yhat2=Yhat3*(year=2);Yhat3=Yhat3*(year=3); if Yhat2=0 then Yhat2=.; if Yhat3=0 then Yhat3=.; If year=2002 then Y2=visits; else Y2=.; Label Y2 = "Visits"; if year=2003 then Y3=visits; else Y3=.; proc print data=SARS; run; *** GRAPHICS *******; goptions reset=all; proc gplot data=SARS; plot visits*julday=yr_p/href=&SARS; symbol1 v=dot i=smooth30; where julday<161; run; proc gplot data=SARS; plot (Y2 Y3 yhat2 yhat3)*julday/href=&SARS overlay; symbol1 v=dot i=none c=green; symbol2 v=dot i=none c=orange; symbol3 v=none i=join c=black w=2; symbol4 v=none i=join c=red w=2; run; goptions reset=all; goptions colors=(black blue red green magenta cyan gray); proc gplot data=sars; plot visits*julday=yr_dow/href=&SARS; symbol1 v=dot i=join w=1; symbol2 v=circle h=1.5 i=join; where julday<161; run; ***** ANALYSIS *******; data SARS; set SARS; if VISITS>.; proc arima data=SARS; * where year=3; identify var=visits crosscor = (d_slope julday SARS perm &days year) noprint; * Permanent effct and change in slope, transitory effects not significant; estimate input = (julday perm d_slope /(1)SARS &days) p=1 ml; run; * Permanent effect and weekday effects only; estimate input = (perm &days) p=2 ml; run; * Standard exponential decay to "business as usual"; estimate input = (julday /(1)SARS &days) p=1 ml; forecast lead=0 id=date interval=1 out=out2; run; **** Residual Graphics *******; proc gplot data=out2; plot residual*date=part; symbol1 v=dot i=join; data out2; set out2; year=year(date); proc arima data=out2; where year=2002; identify var=residual; run; proc arima data=out2; where year=2003; identify var=residual; run; proc arima data=SARS; where year=2002; identify var=visits crosscor = (julday &days) noprint; estimate input = (julday &days) p=1 ml; run; proc format; value pred 3002.0 - 3002.99 = "yr 2002" 3003.0 - 3003.99 = "yr 2003"; data all; merge out2 SARS; by date; keep visits year dow julday date part yrdow EXP; visits = 162.48914 -0.19969*julday -74.65255*EXP +7.87668*D1+0.12568*D2 -8.86646*D3-10.70244*D4 -15.18572*D5-7.23870*D6; dow= year-2000+6; part=2; yrdow = year+1000+dow/10; output; set SARS; part=1; yrdow=year+dow/10; output; proc sort; by part date; proc print; run; goptions colors=(black blue red green magenta cyan gray); proc gplot; plot visits*julday=yrdow; symbol1 v=dot h=1.5 i=none; symbol2 v=circle h=1.8 i=none; symbol3 v=none i=join w=2; format yrdow pred.; run;