************************************************************************
* Program:  C:\NHANES\agest_prev_sas_poptot.sas                              *
* Proposal: Generate population estimates in SAS/SAS Survey            * 
************************************************************************;
libname nh 'C:\NHANES\data\your_dir';
OPTIONS NODATE NOCENTER;

proc format;                                       
                                                   
VALUE SEXFMT   1 = 'Male'                          
               2 = 'Female'                        
                ;                                  
  

VALUE RACEFMT   1 = 'NH-White'
                2 = 'NH-Black'
			    3 = 'Mex-Am'
			    4 = 'Other'
			      ;
                       
VALUE AGEFMT                             
                1 = '20-39'                               
                2 = '40-59'                               
				3 = '60+'       
				;         


value hbpfmt   1 = 'high BP'
               2 = 'low BP'
				  ;
                ; 
VALUE hbp2fmt  100 = 'High BP'
                0 =  'No high BP'
                ;
value selfmt   1 = 'Age ge 20'
               2 = 'Age lt 20'
                ; 
run; 
DATA ANALYSIS_DATA;
SET NH.ANALYSIS_DATA;
  
if ridstatr = 2;  ***examined  ;

if ridageyr ge 20 then sel=1;
else sel=2;
                      
age = .;
if 20 LE ridageyr LE 39 then age=1;
if 40 LE ridageyr LE 59 then age=2;   
if ridageyr GE 60 then age=3;  

race=.;
if ridreth1=3 then race=1;
if ridreth1=4 then race=2; 
if ridreth1=1 then race=3;
if ridreth1=2 or ridreth1=5 then race=4;

*** code to define mean blood pressure measures and HBP ****;
n_sbp = n(of bpxsy1-bpxsy4);
      n_dbp = n(of bpxdi1-bpxdi4);
      *Setting DBP values of 0 as missing for calculating average;
      array _DBP bpxdi1-bpxdi4;
      do over _DBP;
      if (_DBP = 0) then _DBP = .;
      end;  

      mean_sbp = mean(of bpxsy1-bpxsy4);

      mean_dbp = mean(of bpxdi1-bpxdi4); 
      if BPQ050a=1 then HBP_trt=1;

      else if BPQ020 in (1,2) and BPQ050a < 7 then HBP_trt=0;

      if n_sbp>0 and n_dbp>0 then do;
            if mean_sbp>=140 then SBP140=1;
            else SBP140=0;
            if mean_dbp>=90 then DBP90=1;
            else DBP90=0; 
      end;
if HBP_trt>=0 and SBP140>=0 and DBP90>=0 then do;
            if HBP_trt=1 or SBP140=1 or DBP90=1 then HBP=1; 
            else HBP=2;
            end;
if hbp=1 then hbpx=100;
if hbp=2 then hbpx=0;
 
LABEL                           
       age = 'AGE GROUP'    
	   race = 'Race Ethnicity'	
	   riagendr = 'Gender'
       sel = 'Subpop: 1=age ge 20 2=age lt 20'
                     ;                                                     
                                                             
                            
                                                                
RUN;

title;

proc surveymeans data=analysis_data df nobs mean stderr clm cv ;  
strata sdmvstra;
cluster sdmvpsu;
weight wtmec4yr;
class riagendr race;
var hbpx;
domain sel sel*riagendr sel*race sel*riagendr*race;
ods OUTPUT domain(match_all)=unadj; /*dataset output for each domain listed on the domain statement*/
run;

/* append output datasets into one dataset*/
data bp_stats; 
set unadj unadj1 unadj2 unadj3;
if sel=1;                                    /*keep if age ge 20*/
if race=. then race=0; 
if riagendr=. then riagendr=0; 

*Wald 95% confidence interval;
ll=round(lowerclmean,.01);
ul=round(upperclmean,.01);

percent=round(mean,.01);
sepercent=round(stderr,.01);
run;


/* read in population totals from supplied SAS data set cpstot9902 .
 Create age, race, gender subgroups of interest */
/*rename variables*/
data cpstot9902;
	set nh.cpstot9902;
	rename ctutage=age ctutpopt=poptot ctutgndr=riagendr ctutrace=race;
run;

Proc means data=cpstot9902 ; where age >= 20;
var poptot;
output out=d1 n=n sum=sum;
run; 

proc sort data=cpstot9902; by riagendr; run;
proc means data=cpstot9902; where age >= 20;
var poptot;
by riagendr;
output out=d2 n=n sum=sum;
run; 

proc sort data=cpstot9902; by race; run;
proc means data=cpstot9902; where age >= 20;
var poptot;
by race;
output out=d3 n=n sum=sum;
run; 
proc sort data=cpstot9902; by riagendr race; run;
proc means data=cpstot9902; where age >= 20;
var poptot;
by riagendr race;
output out=d4 n=n sum=sum;
run; 

/* append output datasets into one dataset*/
data saspt9902;
set d1 d2 d3 d4;
if race=. then race=0;
if riagendr=. then riagendr=0;
run;

/* sort by subdomains being reported */

proc sort data=bp_stats; by  riagendr race ; run;
proc sort data=saspt9902  ; by riagendr race ; run;

/* merge prevalence and 95% CI with population totals */
data comb;
merge bp_stats(in=a) saspt9902(drop=n);
by riagendr race ;
if a ;

/* calculate and round the population counts by applying the population totals to the % pos 
 and the 95% CI on the % positive */
popmean=(percent/100)*sum;
popl=ll/100*sum;
popu=ul/100*sum;

/* round and format estimates to the nearest thousand */
   
poplr=round(popl,1000);
popur=round(popu,1000);
popmeanr=round(popmean,1000);
totalr=round(sum,1000) ;

format popmean popl popu  11.0 ;
run;

/* formatted printout of results   */

proc print noobs split='/' double;
var  riagendr race percent sepercent ll ul N totalr popmeanr poplr popur ;
format race racefmt. riagendr sexfmt.  n 5.0 percent 5.2 sepercent 5.2
ll 4.2 ul 4.2  df 2.0;
label   
percent='%'/'with'/'high'/'bp'
n='Num'/'bp'/'status'
sepercent='Std'/'error'
ll='Lower'/'95 %'/'Wald'/'CI'
ul='Upper'/'95 %'/'Wald'/'CI'
popmeanr='Pop'/'Est'/'US'/'with'/'high'/'bp'
totalr='Pop'/'total'/'US'
poplr='Pop Est'/'Lower'/'95 %'/'WALD'/'CI'
popur='Pop Est'/'Upper'/'95 %'/'WALD'/'CI' ;
;
title1 'Prevalence of persons with high Bp - US, 1999-2002';
title2 'Percent and population estimates of number with high Bp-Wald CI';
run;