* copy and paste 'srs','est_srs', and 'est_strs' SAS macros into a new program editor;
* then run them before using this program;

* url for data in INFILE;
filename uspop url 'http://webpages.uidaho.edu/~renaes/Data/uspop.csv';

data uspop;
length State $ 15;
* first row to read in is row 3 since row1 has headers and row2 is the whole US;
infile uspop firstobs=3 dlm=',';
input State $ TotalPop Section Pop18_24 PopGE18 Pop15_44 PopGE65 PopGE85 PctPov;
if State = 'Delaware' then totalpop = 807385;
if State = 'Delaware' then pop18_24 = 81501;
run;

proc means sum data = uspop;
var pop18_24;
run;

proc sort data=uspop;
by totalpop;
proc print;
run;

data stratusa ;
input group $ n ;
do i = 1 to n ;
input obs State $ TotalPop Section Pop18_24 PopGE18 Pop15_44 PopGE65 PopGE85 PctPov ;
output ;
end ;
drop i n ;
cards ;
small 5
1 Wyoming         498703 4   54248   376359   210398   59222   7273  8.8000
4 Alaska          643786 4   58738   451358   282205   39200   3073  8.7000
5 SouthDakota     761063 2   82635   565438   319741  108322  17021 10.0000
6 Delaware        807385 3   81501    11821    11821   11821  11821  7.9000
11 Maine          1294464 1  118126  1015406   539991  186383  25025 11.9000
medium 5
23 Oklahoma       3493714 3  377256  2620154  1491234  460459  58325 14.6000
26 SouthCarolina  4107183 3  429425  3128020  1794151  503256  55259 14.7000
33 Maryland       5458137 3  488911  4078212  2380371  616699  73543  7.3000
36 Washington     6068996 4  593628  4555636  2678937  677532  93072 10.8000
41 Georgia        8560310 3  868937  6291833  3956561  813652  95660 12.1000
large 5
43 Michigan      10050446 2  970466  7480182  4303040 1231920 155891 10.5000
44 Ohio          11421267 2 1098431  8541340  4811220 1513372 190926 10.1000
46 Illinois      12600620 2 1228541  9346097  5529191 1499249 206861 11.5000
47 Florida       16713149 3 1403624 12830878  6664700 2854838 360332 12.6000
49 Texas         21779893 3 2287194 15677577  9857869 2152896 255611 15.3000
;
proc print;
run;

proc boxplot data=stratusa;
plot pop18_24*group;
run;

data stratinfo;
input group $ popsize;
cards;
small 21
medium 21
large  8 
;
run;

* Estimates for StRS;
%est_strs(sample=stratusa,strata=group,setup=stratinfo,param=total,response=pop18_24,npop=popsize); 

* Compare to SRS;
%srs(frame=uspop,npop=50,n=16,sample=uspop_srs,seed=020509);

proc print data=uspop_srs;
run;

%est_srs(sample=uspop_srs,npop=50,response=Pop18_24,param=total);

* checking stratum-specific estimates;

data uspopgroup;
set uspop;
popgroup = .;
if totalpop le 3000000 then popgroup = 1;
if totalpop le 10000000 and totalpop gt 3000000 then popgroup = 2;
if totalpop gt 10000000 then popgroup = 3;
run;

proc means mean sum std;
class popgroup; 
var pop18_24;
run;

proc means data=stratusa;
class group;
var pop18_24;
run;