SAS/STAT Examples

data inputdata;
   input los @@;
   datalines;                                                                         
1671 1300   722   586   552   525  364                                                
 359  321   302   272   248   226  216
 208  182   141   141   132   120  117
 115  114   113   104   103   101   99                                          
  96   94    93    92    88    84   83
  81   79    74    70    63    62   62
  61   56    55    53    53    51   51
  50   49    36    33    33    33   29
  28   26    24    19    16    16   15
    14 9 9 5 5 3 3
     2 2 1 24028 23871 22852 22464
22317 21525 21264 21131 21077 20944 20608
20133 19583 19038 18821 18761 18742 18479
18449 18241 18128 17830 17800 17746 17285
16763 16444 15922 15226 15224 15178 14885
14744 14715 14601 14482 14174 14149 13295
13274 12850 12661 12660 12405 12300 12248
11821 11769 11577 11159 10680 10414 10393
10253 10095 9926 9722 9670 9426 9389
 9109 9004 9003 8959 8927 8583 8465
 8464 8317 8176 8156 8114 8096 7884
 7639 6552 6387 6283 6243 6155 5917
 5397 5370 5165 5007 4899 4408 4395
 4290 4282 4208 4182 4153 3855 3810
 3735 3719 3676 3467 3225 3149 3115
 2989 2933 2927 2918 2802 2746 2701
 2542 2430 2421 2380 2340 2331 2318
 2235 2209 2053 1997 1978 1923 1856
 1812 1804 1792 1757 1717 1670 1666
 1661 1615 1610 1609 1525 1481 1424
 1388 1376 1349 1302 1159 1134 1057
 1006 937 756 733 723 588 533
  522 348 212 439 188 181 141
  131 116 96 83 56 40 40
   37 36 26 23 19 14 8
20809 20319 18918 15058 13540 13446 10161
 9163 6474 5420 4584 4341 3805 3640
 3640 3578 3453 3375 3186 2756 2724
 2658 2599 2583 2416 2261 2232 2127
 2126 2035 2015 1986 1903 1897 1874
 1780 1776 1759 1752 1749 1730 1624
 1580 1552 1549 1546 1493 1472 1469
 1456 1441 1433 1428 1399 1391 1378
 1363 1359 1317 1295 1279 1262 1247
 1201 1199 1195 1182 1175 1141 1135
 1106 1049 1037 1031 1027 1006 994
  972 952 944 938 938 936 933
  897 887 880 880 863 842 806
  792 791 789 764 762 726 712
  707 706 703 697 691 687 677
   649 635 613 607 602 602 601
   594 590 576 569 552 544 544
   539 532 530 525 511 509 502
   497 477 462 453 446 443 434
   426 414 411 407 399 397 391
   375 374 370 366 365 359 357
   351 350 348 345 343 341 337
   336 336 333 330 330 310 306
   300 289 285 273 260 252 238
   231 223 219 215 197 197 197
   187 187 161 126 125 123 107
   100 93 83 58 48 37 37
  35   28    20    19     9     5    2
   2    2 22317 14006 11549 11006 8981
8402 7947  7266  6693  4408  4010 4003
1970 1857  1849  1833  1770  1769 1514
1217  956   924   611   386   280   93
;


/*-----------------------------------------------------------------
  Example: Fitting a Mixture of Exponential Distributions for
           Patients Length of Stay 
  Requires: SAS/STAT
  Version: 9.2                                                                       
  ------------------------------------------------------------------*/

** Calculates the Average number of days in hospital stays **;
proc means data = inputdata;
   ods output Summary=summ;
run;

** Creates a macro variable for the mean LOS **;
data _null_;
   set summ;
   call symputx("los_Mean",los_Mean);
run;

** Calculates the Scaled density of Exponential **;
data a;
   do x = 10 to 30000 by 1; 
      ex1 = 300*pdf("expo", x, &los_Mean); 
      output;
   end;
run;

** Merges inputdata and scaled density data set **;
data a_inputdata;
   set a inputdata;
run; 

** Template for displaying Histogram, and both scaled density plots **;
proc template;
   define statgraph sasuser.hist;
      begingraph;
      entrytitle "Distribution and Scaled Posterior Densities of Length of Stays"; 
         layout overlay / border=FALSE 
            xaxisopts=(label="Length of Stay (days)")
            yaxisopts=(label="Scaled Density");   
            seriesplot x = x y = ex1 / primary = true 
               lineattrs=(thickness = 2 color=red) 
               name = '1' legendlabel = "Exponential Density";  
            histogram los / binwidth=300 scale=proportion 
               fillattrs=(color=lightgray transparency = .6); 
            discretelegend '1'/across=2 autoalign=(center)
               location=outside;
         endlayout;
      endgraph;
   end;
run;
   
ods graphics on;   
** Create graph **;
proc sgrender data=a_inputdata template=sasuser.hist;
run;     
ods graphics off;                                             

ods graphics on;
proc mcmc data=inputdata seed=1010 nmc=50000 thin=5 
   propcov=quanew;
   ods output PostSummaries = post_summ;
   parms B 100 D 6000 pi 0.5;
   prior B D ~ igamma(3/10, scale = 10/3);
   prior pi ~ uniform(0,1);
   if (B < D) then
      llike = log(pi*pdf("expo", los, B) + (1-pi)*pdf("expo", los, D));
   else
      llike = .;
   model general(llike);
run;    
ods graphics off;  

** Creates several macro variables for calculating posterior density **;
data _null_;
   set post_summ;
   if(Parameter = "B") then do;
      b_temp = Mean;
      call symputx("Bmean",b_temp);
      end;   
   if(Parameter = "D") then do;
      d_temp = Mean;
      call symputx("Dmean",d_temp);
      end;
   if(Parameter = "pi") then do;
      pi_temp = Mean;
      call symputx("pimean",pi_temp);
      end;
run;

** Creates ex2 variable, scaled density values of mixture **;
data b;
   do x2 = 200 to 30000 by 1; 
      ex2 = 300*(&pimean*pdf("expo", x2, &Bmean) + 
              (1- &pimean)*pdf("expo", x2, &Dmean));  
      output;
   end;
run;

** Creates combined output data sets **;
data b_inputdata;
   set a b inputdata;
run;

** Template for displaying Histogram, and scaled density plots **;
proc template;                                           
   define statgraph sasuser.hist2;
      begingraph;
      entrytitle "Distribution and Scaled Posterior Densities of Length of Stays"; 
         layout overlay / border=FALSE 
            xaxisopts=(label="Length of Stay (days)")
            yaxisopts=(label="Scaled Density");   
         seriesplot x = x y = ex1 /primary = true 
            lineattrs=(thickness = 2 color=red) 
            name = '1' legendlabel = "Exponential Density";  
         seriesplot x = x2 y = ex2 /primary = true 
            lineattrs=(thickness = 2 color=blue)
            name = '2' legendlabel = "Mixture Density";
         histogram los/binwidth=300 scale=proportion 
            fillattrs=(color=lightgray transparency = .6); 
         discretelegend '1' '2'/across=2 autoalign=(center) 
            location=outside;
         endlayout;
      endgraph;
   end;
run;

ods graphics on;   
** Create graph **;
proc sgrender data=b_inputdata template=sasuser.hist2;                        
run;     
ods graphics off;