/****/
/* PROGRAM:  Example SAS Analysis Programs.SAS 	 						  */
/* 																		  */
/* PURPOSE:  THIS PROGRAM PROVIDES EXAMPLE ANALYSIS CODE FOR PRODUCING 	  */
/*			 ESTIMATES FROM THE 2023 NIS-CHILD PUF USING SAS.	 		  */
/*			 THE PROGRAM IS DIVIDED INTO FIVE SECTIONS: 		 		  */
/* 																		  */
/* EXAMPLE 1:   PRODUCE ESTIMATION AREA ESTIMATES AND STANDARD ERRORS 	  */
/*				FOR P_UTD431H314_ROUT_S			 						  */
/* 																		  */
/* EXAMPLE 2:   PRODUCE STATE ESTIMATES AND STANDARD ERRORS 			 */
/*				FOR P_UTD431H314_ROUT_S			 						  */
/* 																		  */
/* EXAMPLE 3:   PRODUCE A TABLE OF HAD_CPOX BY STATE FOR ALL HOUSEHOLD	  */
/*				COMPLETES USING RDDWT_C 								  */
/* 																		  */
/* EXAMPLE 4:   PRODUCE A TABLE OF P_UTD431H314_ROUT_S BY INCPOV1 		  */
/*				BY RACE_K 									 			  */
/* 																		  */
/* EXAMPLE 5:   PRODUCE A BAR CHART OF P_UTD431H314_ROUT_S BY 			  */
/*				INCPOV1 BY RACE_K 								 		  */
/****/

options ps = 78 ls = 90 obs = max;

*--- ASSIGN SAS LIBRARIES AND FILE NAMES ---*; 
libname dd 'c:\nispuf23'; *--- Specify path to sas dataset ---*;
libname library 'c:\nispuf23'; *--- If dataset was created with formats stored ---*;
%let out = 'c:\nispuf23'; *--- Specify output location ---*; 
libname out &out.;

*--- PERMANENTLY SPECIFY PATH TO LIBRARY ---*;
*--- OTHERWISE COMMENT THIS STATEMENT OUT ---*;
%let yy = 23; *--- YEAR OF DATA COLLECTION ---*;
%let in_file = dd.nispuf&yy.; *--- NAME OF SAS DATASET ---*;
%let ESTIAP = ESTIAP&yy.; * --- ESTIMATION AREA VARIABLE TO USE ---*;
%let RDDWT = RDDWT_C; * --- Weight to use for household variables (RDDWT_C is the single-frame cellular phone weight excluding territories. Use RDDWT_C_TERR to include territories) ---*;
%let PROVWT = PROVWT_C; * --- Weight to use for provider variables (PROVWT_C is the single-frame cellular phone weight excluding territories. Use PROVWT_C_TERR to include territories) ---*;
%let STRAT = STRATUM; * --- STRATUM variable to use for variance estimation (use STRATUM for all estimation) ---*;



*--- DEFINE FORMATS USED FOR ANALYSIS ---*;
proc format;

	value putd431h314f
		1 = "4:3:1:3:3:1:4 Up-to-Date"
		0 = "Not 4:3:1:3:3:1:4 Up-to-Date"
		;

	value estiapf
		0 = "U.S. Total"
		1 = "CT"
		2 = "MA"
		4 = "ME"
		5 = "NH"
		6 = "RI"
		7 = "VT"
		8 = "NJ"
		10 = "NY-Rest of State"
		11 = "NY-City of New York"
		12 = "DC"
		13 = "DE"
		14 = "MD"
		16 = "PA-Rest of State"
		17 = "PA-Philadelphia County"
		18 = "VA"
		19 = "WV"
		20 = "AL"
		22 = "FL"
		25 = "GA"
		27 = "KY"
		28 = "MS"
		29 = "NC"
		30 = "SC"
		31 = "TN"
		34 = "IL-Rest of State"
		35 = "IL-City of Chicago"
		36 = "IN"
		38 = "MI"
		40 = "MN"
		41 = "OH"
		44 = "WI"
		46 = "AR"
		47 = "LA"
		49 = "NM"
		50 = "OK"
		51 = "TX-Rest of State"
		54 = "TX-City of Houston"
		55 = "TX-Bexar County"
		56 = "IA"
		57 = "KS"
		58 = "MO"
		59 = "NE"
		60 = "CO"
		61 = "MT"
		62 = "ND"
		63 = "SD"
		64 = "UT"
		65 = "WY"
		66 = "AZ"
		68 = "CA"
		72 = "HI"
		73 = "NV"
		74 = "AK"
		75 = "ID"
		76 = "OR"
		77 = "WA"
		106 = "Puerto Rico"
		;
		
	value statef
		0 = "U.S. Total"
		1 = "Alabama"
		2 = "Alaska"
		4 = "Arizona"
		5 = "Arkansas"
		6 = "California"
		8 = "Colorado"
		9 = "Connecticut"
		10 = "Delaware"
		11 = "District of Columbia"
		12 = "Florida"
		13 = "Georgia"
		15 = "Hawaii"
		16 = "Idaho"
		17 = "Illinois"
		18 = "Indiana"
		19 = "Iowa"
		20 = "Kansas"
		21 = "Kentucky"
		22 = "Louisiana"
		23 = "Maine"
		24 = "Maryland"
		25 = "Massachusetts"
		26 = "Michigan"
		27 = "Minnesota"
		28 = "Mississippi"
		29 = "Missouri"
		30 = "Montana"
		31 = "Nebraska"
		32 = "Nevada"
		33 = "New Hampshire"
		34 = "New Jersey"
		35 = "New Mexico"
		36 = "New York"
		37 = "North Carolina"
		38 = "North Dakota"
		39 = "Ohio"
		40 = "Oklahoma"
		41 = "Oregon"
		42 = "Pennsylvania"
		44 = "Rhode Island"
		45 = "South Carolina"
		46 = "South Dakota"
		47 = "Tennessee"
		48 = "Texas"
		49 = "Utah"
		50 = "Vermont"
		51 = "Virginia"
		53 = "Washington"
		54 = "West Virginia"
		55 = "Wisconsin"
		56 = "Wyoming"
		72 = "Puerto Rico"
		;
	
	value hadcpoxf
		1 = "Yes"
		2 = "No"
		;

	value incpvr2f
		1 = "Above Poverty, > $75,000"
		2 = "Above Poverty, < = $75,000"
		3 = "Below Poverty"
		4 = "Unknown"
		;
		
	value race_kf
		1 = "White Only"
		2 = "Black Only"
		3 = "Other/Mult Race"
		;
run;

data nispuf;
	set &in_file. (keep = SEQNUMHH SEQNUMC &ESTIAP. &RDDWT. &PROVWT. &STRAT. STATE P_UTD431H314_ROUT_S HAD_CPOX RACE_K INCPOV1);
	format &ESTIAP. ESTIAPf. 
		   STATE statef.
		   P_UTD431H314_ROUT_S putd431h314f.	 
		   HAD_CPOX hadcpoxf.
		   INCPOV1 incpvr2f.
		   RACE_K race_kf.;
run;





**
EXAMPLE 1:	PRODUCE ESTIMATION AREA ESTIMATES AND STANDARD ERRORS FOR 
			P_UTD431H314_ROUT_S
			VACCINATIONS (P_UTD431H314_ROUT_S)
**;

proc sort data = nispuf;by &estiap.;run;

ods exclude all;
ods output Statistics = sas_est;
proc surveymeans data = nispuf nobs sum mean stderr;
	stratum &STRAT.;
	cluster SEQNUMHH;
	weight &PROVWT.;
	class P_UTD431H314_ROUT_S;
	var P_UTD431H314_ROUT_S;
	by &ESTIAP.;
run;
ods exclude none;

data sas_est;
	set sas_est;
	mean = mean*100; *Convert to percent estimates;
	stderr = stderr*100;
run;

proc print data = sas_est (where = (varlevel = '4:3:1:3:3:1:4 Up-to-Date')) noobs label;
	format mean stderr 5.2;
	var &estiap. mean stderr;
	label mean = 'Percent 4:3:1:3:3:1:4 Up-to-Date'
		  stderr = 'Standard Error';
	title "4:3:1:3:3:1:4 Estimates by Estimation Area";
run;





EXAMPLE 2:	PRODUCE STATE ESTIMATES AND STANDARD ERRORS FOR 4+ HIB OF ANY 
			TYPE OR 2 HIB MERCK+1 HIB OF ANY TYPE VACCINATIONS 
			(P_UTD431H314_ROUT_S) 					
;

proc sort data = nispuf;by state;run;

ods exclude all;
ods output Statistics = sas_est2;
proc surveymeans data = nispuf nobs sum mean stderr;
	stratum &STRAT.;
	cluster SEQNUMHH;
	weight &PROVWT.;
	class P_UTD431H314_ROUT_S;
	var P_UTD431H314_ROUT_S;
	by STATE;
run;
ods exclude none;

data sas_est2;
	set sas_est2;
	mean = mean*100; *Convert to percent estimates;
	stderr = stderr*100;
run;

proc print data = sas_est2 (where = (varlevel = '4:3:1:3:3:1:4 Up-to-Date')) noobs label;
	format mean stderr 5.2;
	var STATE mean stderr;
	label mean = 'Percent 4:3:1:3:3:1:4 Up-to-Date'
		  stderr = 'Standard Error';
	title "4:3:1:3:3:1:4 Estimates by State";
run;





EXAMPLE 3:	PRODUCE A TABLE OF HOUSEHOLD REPORT OF THE CHILD HAVING CHICKEN 
			POX BY STATE FOR ALL HOUSEHOLD COMPLETES USING RDDWT_C 
;

proc sort data = nispuf;by state;run;

ods exclude all;
ods output Statistics = sas_est3;
proc surveymeans data = nispuf nobs sum mean stderr;
	stratum &STRAT.;
	cluster SEQNUMHH;
	weight &RDDWT.;
	class HAD_CPOX;
	var HAD_CPOX;
	by STATE;
run;
ods exclude none;

data sas_est3;
	set sas_est3;
	mean = mean*100; *Convert to percent estimates;
	stderr = stderr*100;
run;

proc print data = sas_est3 (where = (varlevel = 'Yes')) noobs label;
	format mean stderr 5.2;
	var STATE mean stderr;
	label mean = 'Percent HAD_CPOX = Yes'
		  stderr = 'Standard Error';
	title "Child Had Chicken Pox by State";
run;





EXAMPLE 4:	PRODUCE A TABLE OF P_UTD431H314_ROUT_S BY INCPOV1 BY RACE_K
;

proc sort data = nispuf;by incpov1 race_k;run;

proc freq data = nispuf;
	where &PROVWT. ne .;
	tables P_UTD431H314_ROUT_S incpov1 race_k;
	title1 "Table 4A. Q1/20&yy. - Q4/20&yy.: Unweighted Frequencies";
run;

ods exclude all;
proc surveymeans data = nispuf nobs sum mean stderr;
	ods output Domain = sas_est4;
	stratum &STRAT.;
	cluster SEQNUMHH;
	weight &PROVWT.;
	class P_UTD431H314_ROUT_S;
	var P_UTD431H314_ROUT_S;
	domain INCPOV1*RACE_K;
run;
ods exclude none;

data sas_est4;
	set sas_est4;
	mean = mean*100; *Convert to percent estimates;
	stderr = stderr*100;
run;

proc print data = sas_est4 (where = (varlevel = '4:3:1:3:3:1:4 Up-to-Date')) noobs label;
	format mean stderr 5.2;
	var INCPOV1 RACE_K mean stderr;
	label mean = '4:3:1:3:3:1:4 Up-To-Date'
		  stderr = 'Standard Error';
	title1 "Table 4B. Q1/20&yy. - Q4/20&yy.: 4:3:1:3:3:1:4 Estimates and Standard Errors by INCPOV1 by RACE_K";
run;

data out.sas_est4;
	set sas_est4 (where = (varlevel = '4:3:1:3:3:1:4 Up-to-Date'));
	keep INCPOV1 RACE_K mean;
	label mean = '4:3:1:3:3:1:4 Up-to-Date';
	format mean 5.2; 
run;





EXAMPLE 5:	PRODUCE A BAR CHART OF P_UTD431H314_ROUT_S BY INCPOV1 BY RACE_K
;

data sas_est4;
	set out.sas_est4;
	format mean 3. 
		   INCPOV1 incpvr2f.
		   RACE_K race_kf.;
	label RACE_K = 'Race of Child'
		  INCPOV1 = 'Poverty Status';
	filename odsout &out.;
	ods listing close;
	
	*Set the graphics environment;
	goptions reset = global gunit = pct border ftext = swissb htitle = 4 htext = 1.5 device = gif;
	ods html body = 'graph_4.html' path = odsout;
run;

title1 h = 12pt "Percentage of Children Up-to-Date with 4:3:1:3:3:1:4 Vaccine Series";
title2 h = 12pt "by Race and Poverty Status, National Immunization Survey - Child, 20&yy.";
footnote j = r 'graph_4';

proc sgplot data = sas_est4;
	styleattrs datacolors = (wheat lightpink forestgreen) datacontrastcolors = (wheat lightpink forestgreen) datalinepatterns = (solid);
	vbar INCPOV1 / response = mean group = RACE_K groupdisplay = cluster stat = mean barwidth = 0.8;
	xaxis display = (noticks);
	yaxis grid;
	keylegend / across = 1 position = top;
run;

ods html close;
ods listing;