/* ======================================================*/
/*===========NYC HANES 2013-14 AND 2004 TUTORIAL===========*/
/* ======================================================*/
/* Created by NYC-HANES 2013-14 Team: 
	Claudia Chernov 					(cchernov@health.nyc.gov)
	Rania Kanchi 						(rania.kanchi@nyumc.org)
	Lorna Thorpe							(lorna.thorpe@nyumc.org)
	Sharon Perlman 					(sperlma1@health.nyc.gov)

Thanks to Jesica Rodriguez-Lopez for her contribution in creating this sample code
/* ======================================================*/

/*======================DISCLOSURE======================*/
/*All material included in this tutorial presentation is available for public use 
as an educational tool.
Results from this tutorial analysis should not be used as official estimates 
or references for research.
/*=======================================================*/

/*=======================================================*/
/*This tutorial Includes:	======================================
- Recoding of variables	======================================
- Weight adjustment due to non-response 	========================
- Age standardization		======================================
- Prevalence estimate calculations	==============================
- Stratified analysis (Prevalence estimates by gender)	================
- Comparison across years		==================================*/
/*=======================================================*/

/*=======================================================*/
/* DIABETES PREVALENCE ANALYSIS USSING 2004 AND 2013-14 DATA */
/* Objective:  Calculate the prevalence of Diabetes in NYC residents, total and 
	by gender, and compare diabetes prevalence across survey years======*/
/*=======================================================*/

/*======================== 2013-14 ========================*/
/*================= PREVALENCE ESTIMATES =================*/
/* Step 1:  Call in dataset  */
/*Download datasets and format code, and save on your own computer*/
/*For this tutorial, we will assume that datasets and format code are saved in C:\Data\2014*/

libname NYCH14 'C:\Data\2014';
%include 'C:\Data\2014\SAS variable formats_V2_0916.sas';
run;

data NYCHANES14;
	set NYCH14.NYCHANES_2013_14_V2_09212016;
run;

/* Step 2:  Recode/create outcome of interest  */

data dbts14;
set NYCHANES14;

data dbts14;
set NYCHANES14;

*Define diabetes using FPG or A1c or previous diagnosis of diabetes;

/*Diabetic: Fasted and have plasma glucose =>126 or A1C =>6.5 or self-reported diagnosis*/
if DBT_Fasting=1 and (Glucose >=126 or A1c >=6.5) or DIQ_1=1 then dbts=1; 

/*Diabetic: Did not fast and have A1c =>6.5 or self-reported diagnosis*/
else if DBT_Fasting=2 and (A1c =>6.5 or DIQ_1=1) then dbts=1; 

/*Non-Diabetic: Fasted and have plasma glucose<126 and A1c<6.5 and no self-reported diagnosis*/
else if DBT_Fasting=1 and (0<Glucose<126 or 0<A1c<6.5)  and DIQ_1 ne 1 then dbts=2;  

/*Non-Diabetic: Did not fast and have A1c<6.5 and no self-reported diagnosis*/
else if DBT_Fasting=2 and 0<A1c<6.5  and DIQ_1 ne 1 then dbts=2;  

label dbts = 'Diabetes based on FPG, A1c, or Self-Report';

/*	Step 4:  Adjust for non-response	*/
 /*1. Initial Diabetes weight. This is the first crucial step in adjusting for non-response. 
		Set the weight to zero for cases that will not be included in the analysis (that is, 
		cases with a missing value for the dbts variable).*/
if dbts=. then dbts_Wt=0; 
/*2. Set the weight to the identified appropriate weight for cases that will be included 
		in the analysis. (We will use CAPI_WT, because our definition is inclusive.)*/
else dbts_Wt =CAPI_WT; 

run;


/*Step 3:  Identify appropriate weight	*/
/*Number of cases to be included in the analysis*/
proc freq data=dbts14;
	tables dbts/list;
run;

/*	Step 4:  Adjust for non-response	*/
/* 2. Create a new dataset called SAMPLEWTS, in which you sum the initial diabetes weight 
	by age, sex, race/ethnicity	*/
proc summary data=dbts14 nway ;
	var dbts_Wt;
	class Agegrp4c gender race;
	output out=SAMPLEWTS sum=samplewt; /*3. Create a new variable called samplewt 
																		(this will be the sum of the initial diabetes 
																		weight for all observations by age,sex, and 
																		race/ethnicity)	*/
run;

proc print data=samplewts noobs;
run;
/*4. Call in the appropriate control total dataset	*/
data Control_capi; 
	set NYCH14.Control_capi;
run;
/*5. Sort the control total dataset and the samplewts dataset by the variables 
		on which they will merge:Agegrp4c gender race 	*/
proc sort data=samplewts; 
	by Agegrp4c gender race;
run;
proc sort data=Control_capi; 
	by Agegrp4c gender race;
run;

/*6. Merge the dataset that contains weight control totals with the dataset of 
	summed weights (from the analytic dataset), and create weight factor	*/

data wtfactor; 
	merge samplewts Control_capi; 
	by Agegrp4c gender race;
	wt_factor_dbts=control_totals/samplewt; /*weight factor*/
	keep Agegrp4c gender race wt_factor_dbts;
run;

proc print data=wtfactor noobs; 
run;
	
/*7. 	Sort the new weight factor dataset and your analytic dataset, and then 
		merge the two datasets. Create the weight for your analysis by multiplying 
		the weights of participants with non-missing data by the weight factor	*/

proc sort data=wtfactor; 
	by Agegrp4c gender race;
run;

proc sort data=dbts14; 
	by Agegrp4c gender race;
run;

data dbts14_2; 
merge dbts14 (in=a) wtfactor; 
	by Agegrp4c gender race;
	if a;
	/*Adjusted Diabetes Weight*/
	dbts_wt_new = dbts_Wt * wt_factor_dbts; 
	drop dbts_Wt wt_factor_dbts;
	/*Adjusted Diabetes Weight Indicator*/
	if dbts_wt_new = 0 then wt_dbts = 2;
		else wt_dbts = 1;
	label dbts_wt_new = 'Diabetes Weight'
			wt_dbts = 'Diabetes Weight Indicator';
	format wt_dbts TE_3F.;
run;


/*Step 5:  Test your new weight	*/
/* 1. Make sure your adjusted weight is defined for all the cases included 
	in the analysis: 1246	*/
proc freq data = dbts14_2;
	table wt_dbts;
run;

/* 2. Make sure total weighted frequency is the NYC total population 6,285,749, 
	and the gender distribution matches that of the total sample: 46.9% males and
	53.1% females */

proc freq data = dbts14_2;
	table gender;
	weight dbts_wt_new;
run;

/*Step 6:  Generate prevalence estimates 	*/
/*SUDAAN	*/
proc descript data=dbts14_2 filetype=sas design=wr NOTSORTED;
	nest Borostratum PSUnest HHnest; /*Sampling design variables*/
	weight dbts_wt_new; /*Weight*/
	var dbts dbts; /*Variable to be analyzed. State the variable as many times as the categories.*/
	catlevel 1 2; /*State the categories of the variable*/
	tables _one_ ; /*Prevalence estimate for entire sample*/
	subgroup dbts AGEGROUP; /*State all variables in your analysis*/
	levels  2 3; /*For each variable in subgroup, state the number of categories*/	
	stdvar AGEGROUP; /*Age standardization variable*/
	stdwgt  0.396579  0.371795  0.231626; /*Standardizing proportions*/
	/*The following concern the appearance of the results*/
	rtitle "Age-Standardized Weighted Prevalence of Diabetes, NYC HANES 2013-14 ";
	setenv COLWIDTH=16 COLSPCE=2; 
	print  NSUM="SAMSIZE" WSUM="POPSIZE" total="weights for %"  percent lowpct uppct sepercent / style = nchs nsumfmt=F8.0 wsumfmt=F8.0 totalfmt=F8.0 sepercentfmt=F3.1 percentfmt=F6.1 lowpctfmt=F6.1 uppctfmt=F6.1 ; 
run;

/*Generate prevalence estimates*/
/*SAS	*/
/*  In SAS you can use the proc surveyreg to get weighted, age-standardized prevalence 
	estimates.*/
/*  proc surveyreg is an analogous procedure to proc reg. Both will give you linear regression 
	estimates. */
/*  Since your outcome variable (dbts) is dichotomous, you need to recode so that a negative 
	outcome is coded as 0 and a positive outcome is coded as 100.  That will make it easier to 
	read and interpret the results from the proc surveyreg.*/

data dbts14_2;
      set dbts14_2;
      if dbts = 2 then dbts_10=0;
      else if dbts= 1 then dbts_10 = 100;
run;

proc surveyreg data=dbts14_2 nomcar;
stratum Borostratum ;/*Sampling design variable*/
cluster PSUnest HHNest;/*Sampling design variable*/
weight dbts_wt_new; /*weighting variable*/
class AGEGROUP;/*Age standardization variable*/
model dbts_10 = AGEGROUP/noint solution vadjust=none CLPARM; /*To get prevalence estimate*/
estimate AGEGROUP  0.396579 0.371795 0.231626;/*standardizing proportions*/
/*The following lines are all about appearance of the results*/
TITLE1 "Age-Standardized Weighted Total Prevalence of Diabetes ";
TITLE2 "NYC-HANES 2013-14";
run;

/*================= STRATIFIED ANALYSIS =================*/
/*SUDAAN*/

/* 1. Gender specific prevalence estimates	*/

proc descript data=dbts14_2 filetype=sas design=wr NOTSORTED;
nest Borostratum PSUnest HHnest; /*Sampling design variables*/
weight dbts_wt_new; /*Weight*/
var dbts; /*Variable to be analyzed*/
catlevel 1; /*State the category of the positive outcome*/
tables gender; /*Prevalence estimate by gender*/
subgroup dbts gender AGEGROUP; /*State all variables in your analysis*/
levels  2 2 3; /*State number of categories for each variable*/	
stdvar AGEGROUP; /*Adge standardization variable*/
stdwgt 0.396579 0.371795 0.231626; /*Standardizing proportions*/
/*The following concern the appearance of the results*/
rtitle "Age-Standardized Weighted Prevalence of Diabetes by Gender, NYC HANES 2013-14 ";
setenv COLWIDTH=16 COLSPCE=2; 
print  NSUM="SAMSIZE" WSUM="POPSIZE" total="weights for %"  percent lowpct uppct  sepercent/  style = nchs nsumfmt=F8.0 wsumfmt=F8.0 totalfmt=F8.0 sepercentfmt=F3.1 percentfmt=F6.1 lowpctfmt=F6.1 uppctfmt=F6.1; 
run;

/* 2. Statistical testing	*/
proc descript data=dbts14_2 filetype=sas design=wr NOTSORTED;
	nest Borostratum PSUnest HHNest; /*Sampling design variables*/
		weight dbts_wt_new; /*Weight*/
		var dbts; /*Variable to be analyzed*/
		catlevel 1; /*State the category of the positive outcome*/
		class gender AGEGROUP; /*State all variables in your analysis*/
		pairwise gender; /*Prevalence estimate by gender*/
		stdvar AGEGROUP; /*Age standardization variable*/
		stdwgt 0.396579 0.371795 0.231626; /*Standardizing proportions*/
		/*The following concern the appearance of the results*/
		rtitle "Age-Standardized Weighted Prevalence of Diabetes by Gender, NYC HANES 2013-14";
		setenv COLWIDTH=16 COLSPCE=2; 
		print  NSUM="SAMSIZE" WSUM="POPSIZE" total="weights for %"  percent 
		   lowpct uppct sepercent/ style = nchs nsumfmt=F8.0 wsumfmt=F8.0 totalfmt=F8.0 
		   sepercentfmt=F3.1 percentfmt=F6.1 lowpctfmt=F6.1 uppctfmt=F6.1; 
run;


/*SAS*/

/* 1. Gender specific prevalence estimates and statistical testing*/
proc surveyreg data=dbts14_2 nomcar;
	stratum Borostratum; /*Sampling design variable*/
	cluster PSUnest HHNest; /*Sampling design variable*/
	weight dbts_wt_new; /*Weighting variable*/
	class gender AGEGROUP; /*Stratification and age standardization variables*/
	model dbts_10 = Gender AGEGROUP gender*AGEGROUP/noint solution vadjust=none CLPARM; /*Age-standardized prevalence estimate by Gender*/
	estimate '1: Male' gender 1 0 AGEGROUP  0.396579 0.371795 0.231626
				gender*AGEGROUP 0.396579 0.371795 0.231626 0 0 0;
	estimate '1: Female' gender 0 1 AGEGROUP  0.396579 0.371795 0.231626
				gender*AGEGROUP 0 0 0 0.396579 0.371795 0.231626;
	/*The following concern the appearance of the results*/
	TITLE1 "Age-Standardized Weighted Total Prevalence of Diabetes by Gender";
	TITLE2 "NYC HANES 2013-14";
run;

/*========================== 2004 =========================*/
/*================= PREVALENCE ESTIMATES =================*/

/* Step 1:  Call in dataset  */
/*Download datasets and save*/
/*For this tutorial, we will assume that datasets are saved in  C:\Data\2004*/
libname NYCH04 'C:\Data\2004';
%include 'C:\Data\2014\SAS-variable-formats.sas';

/*Each dataset must be sorted on merging variable SP_ID*/
proc sort data = NYCH04.spfile out=spfile; by sp_id;
proc sort data = NYCH04.capi out=capi; by sp_id;
proc sort data = NYCH04.blood out=blood; by sp_id;
data NYCH04; 
	merge spfile (in=a) capi blood; 
	by sp_id;
	if a; /*Chooses observations that are in the SPfile */
run;


/* Step 2:  Recode/create outcome of interest  */
data dbts04;
set NYCH04; 

*Define diabetes using FPG or A1c or previous diagnosis of diabetes;

/*Diabetic: Fasted and have plasma glucose >=126 or A1C>=6.5 or self-reported diagnosis*/
if WTSF1F>0 and (lbxglu>=126 or lbxgh>=6.5)  or diq010=1 then dbts=1; 

/*Diabetic: Did not fast and have A1c>=6.5 or self-reported diagnosis*/
else if WTSF1F<=0 and lbxgh>=6.5 or diq010=1 then dbts=1; 

/*Non-Diabetic: Fasted and have plasma glucose<126 and A1C<6.5 and no self-reported diagnosis*/
else if WTSF1F>0 and (0<lbxglu<126 or 0<lbxgh<6.5)  and diq010 ne 1 then dbts=2;  

/* Non-Diabetic:  Did not fast and have A1C<6.5 and no self-reported diagnosis */
else if WTSF1F<=0 and 0<lbxgh<6.5  and diq010 ne 1 then dbts=2;  

label dbts = 'Diabetes based on FPG, A1c, or Self-Report';

/*	Step 4:  Adjust for non-response	*/
/*Initial Diabetes weight*/
	/*this is the first step in adjusting for non-response, and is KEY*/
	/*1. Set to zero the weight for cases who will not be included in the analysis (Cases missing dbts)*/
	if dbts = . then dbts_Wt=0; 
	/*2. Set to the identified appropriate weight (WTSF1CH) the weight for cases who will be included in the analysis*/
	else dbts_Wt = WTSF1CH;
run;


proc freq data=dbts04;
	tables dbts/list;
run;

/*	Step 4:  Adjust for non-response	*/
/* 2. Create a new dataset called SAMPLEWTS, 
	in which you sum the {clinic + home} weight by age, sex, race/ethnicity */

proc summary data=dbts04 nway ;
	var  dbts_Wt ; *use modified clinic + home weight created in previous step; 
	class agewt riagendr racewt;
	output out = samplewts sum = samplewt; /*3. Create a new variable called samplewt 
																		(this will be the sum of the {clinic + home} 
																		weight for all observations by age, sex, 
																		and race/ethnicity):	*/
run;

proc print data=samplewts noobs;
run;
	
/*4. Call in the appropriate control total dataset	*/
data int_con_tots; 
	set NYCH04.Clinic_home_totals; 
run;
/*5. Sort the control total dataset and the samplewts dataset by 
		the variables on which they will merge: Agegrp4c gender race	*/
proc sort data=int_con_tots; *clinic + home control totals;
by agewt riagendr racewt;
run;
proc sort data=samplewts;	 *analysis dataset;
by agewt riagendr racewt;
run;

/*6. Merge the dataset that contains weight control totals with the dataset of 
	summed weights (from the analytic dataset), and create weight factor	*/

data wtfactor; 
	merge samplewts int_con_tots; 
	by agewt riagendr racewt;
	wt_factor_dbts=clinic_home_control_total/samplewt;/*weight factor*/
	keep agewt riagendr racewt wt_factor_dbts;
run;

proc print data = wtfactor noobs; 
run;

/*7. 	Sort the new weight factor dataset and your analytic dataset, and then merge the 
		two datasets. Create the weights for your analysis by multiplying the weights of 
		participants with non-missing data by the weight factor.	*/

proc sort data=wtfactor; 
by agewt riagendr racewt;
run;

proc sort data=dbts04; 
by agewt riagendr racewt;
run;

data dbts04_2; 
merge dbts04 (in=a) wtfactor; 
by agewt riagendr racewt; 
dbts_wt_new = dbts_Wt * wt_factor_dbts; /*Adjusted Diabetes Weight*/
if dbts_wt_new = 0 then wt_dbts = 2; /*Adj-diabetes Weight Indicator*/
else wt_dbts = 1;
label 	dbts_wt_new = 'Diabetes Weight'
wt_dbts = 'Diabetes Weight Indicator';
format wt_dbts TE_3F.;
drop dbts_Wt wt_factor_dbts;
run;
/*Step 5:  Test your new weight	*/
/* 1. Make sure your adjusted weight is defined for all the cases included 
	in the analysis: 1823	*/
proc freq data = dbts04_2;
	table wt_dbts;
run;

/* 2. Make sure total weighted frequency is the total population represented by 
		NYC HANES, 5,827,719, and the gender distribution matches that of the survey, 
		46.1% male and 53.9% females */
proc freq data = dbts04_2;
	table riagendr;
	weight dbts_wt_new;
run;

/*  Sampling Design Variables	/	
/*  Create HHNest variable in 2004 dataset	*/
data dbts04_2;
set dbts04_2;
HHNest = _n_; /*This ensures the HHNest has a unique value for each case*/
run;

/*Step 6:  Generate prevalence estimates 	*/
/*SUDAAN	*/
proc descript data=dbts04_2 filetype=sas design=wr NOTSORTED;
	nest stratum psu; /*Sampling design variables*/
	weight dbts_wt_new; /*Weight*/
	var dbts dbts; /* Variable to be analyzed. State the variable as many times as the categories*/ 
	catlevel 1 2; /*State the categories of the variable*/
	tables _one_; /*Prevalence estimate for entire sample*/
	subgroup dbts ageadj; /*State all variables in your analysis*/
	levels  2 3; /* For each variable in subgroup, state the number of categories*/	
	stdvar ageadj; /*Age standardization variable*/
	stdwgt 0.396579 0.371795 0.231626; /*Standardizing proportions*/
	/*The following concern the appearance of the results*/
	rtitle "Age-Standardized Weighted Prevalence of Diabetes, NYC HANES 2004";
	setenv COLWIDTH=16 COLSPCE=2; 
	print  NSUM="SAMSIZE" WSUM="POPSIZE" total="weights for %"  
			percent lowpct uppct sepercent/  style = nchs nsumfmt=F8.0 wsumfmt=F8.0 
			totalfmt=F8.0 sepercentfmt=F3.1 percentfmt=F6.1 lowpctfmt=F6.1 uppctfmt=F6.1; 
run;

/*  Generate prevalence estimates*/
/*  SAS	*/
/*  In SAS you can use the proc surveyreg to get weighted, age-standardized prevalence 
	estimates.
	proc surveyreg is an analogous procedure to proc reg. Both will give you linear regression 
	estimates.
	Since your outcome variable (dbts) is dichotomous, you need to recode so that a negative 
	outcome is coded as 0 and a positive outcome is coded as 100.  That will make it easier to 
	read and interpret the results from the proc surveyreg.	*/
/*Recode dbts to 0 for the negative outcome and 100 for the positive one*/
data dbts04_2;
      set dbts04_2;
      if dbts = 2 then dbts_10=0;
      else if dbts= 1 then dbts_10 = 100;
run;

proc surveyreg data=dbts04_2 nomcar;
	stratum stratum ;/*Sampling design variable*/
	cluster PSU;/*Sampling design variable*/
	weight dbts_wt_new; /*weighting variable*/
	class ageadj;/*Age standardization variable*/
	model dbts_10 = ageadj/noint solution vadjust=none CLPARM; /*To get prevalence estimate*/
	estimate ageadj  0.396579 0.371795 0.231626;/*standardizing proportions*/
	/*The following lines are all about appearance of the results*/
	TITLE1 "Age-Standardized Weighted Total Prevalence of Diabetes ";
	TITLE2 "NYC-HANES 2004";
run;


/*================= STRATIFIED ANALYSIS =================*/
/*SUDAAN*/

/* 1. Gender specific prevalence estimates	*/
proc descript data=dbts04_2 filetype=sas design=wr NOTSORTED;
nest stratum psu; /*Sampling design variables*/
weight dbts_wt_new; /*Weight*/
var dbts; /*Variable to be analyzed.*/
catlevel 1; /*State the categories of the positive outcome*/
tables riagendr; /*Prevalence estimate by gender*/
subgroup dbts riagendr ageadj; /*State all variables in your analysis*/
levels  2 2 3; /* For each variable in subgroup, state the number of categories*/	
stdvar ageadj; /*Age standardization variable*/
stdwgt 0.396579 0.371795 0.231626; /*Standardizing proportions*/
/*The following concern the appearance of the results*/
rtitle "Age-Standardized Weighted Prevalence of Diabetes, NYC HANES 2004";
setenv COLWIDTH=16 COLSPCE=2; 
print  NSUM="SAMSIZE" WSUM="POPSIZE" total="weights for %"  
		percent lowpct uppct sepercent/  style = nchs  nsumfmt=F8.0 wsumfmt=F8.0 
		totalfmt=F8.0 sepercentfmt=F3.1 percentfmt=F6.1 lowpctfmt=F6.1 uppctfmt=F6.1; 
run;


/* 2. Statistical testing	*/
proc descript data=dbts04_2 filetype=sas design=wr NOTSORTED;
 nest stratum psu; /*Sampling design variables*/
 weight dbts_wt_new; /*Weight*/
 var dbts; /*Variable to be analyzed*/
 catlevel 1; /*State the category of the positive outcome*/
 class  riagendr ageadj; /*State all variables in your analysis*/
 pairwise riagendr; /*Prevalence estimate by gender*/
 stdvar ageadj; /*Age standardization variable*/
 stdwgt 0.396579 0.371795 0.231626; /*Standardizing proportions*/
/*The following concern the appearance of the results*/
rtitle "Age-Standardized Weighted Prevalence of Diabetes by Gender, NYC HANES 2004";
setenv COLWIDTH=16 COLSPCE=2; 
print  NSUM="SAMSIZE" WSUM="POPSIZE" total="weights for %"  percent 
   lowpct uppct sepercent/  style = nchs nsumfmt=F8.0 wsumfmt=F8.0 totalfmt=F8.0 
   sepercentfmt=F3.1 percentfmt=F6.1 lowpctfmt=F6.1 uppctfmt=F6.1; 
run;

/*SAS*/

/* 1. Gender specific prevalence estimates and statistical testing*/
proc surveyreg data=dbts04_2 nomcar;
stratum stratum ; /*Sampling design variable*/
cluster PSU; /*Sampling design variable*/
weight dbts_wt_new; /*Weighting variable*/
class   riagendr ageadj; /*Stratification and age standardization variables*/
model dbts_10 = riagendr ageadj riagendr*ageadj/noint solution vadjust=none CLPARM; /*To get age standardized prevalence estimate by Gender*/
estimate '1: Male'  riagendr 0 1 ageadj  0.396579 0.371795 0.231626
	riagendr*ageadj 0 0 0 0.396579 0.371795 0.231626;
estimate '2: Female' riagendr 1 0 ageadj  0.396579 0.371795 0.231626
	riagendr*ageadj 0.396579 0.371795 0.231626 0 0 0;
/*The following concern the appearance of the results*/
TITLE1 "Age-Standardized Weighted Total Prevalence of Diabetes by Gender";
TITLE2 "NYC HANES 2004";
run;



/*================== COMPARISON ACROSS YEARS==================*/

/*Step 1: Recoding variables */
/* For each survey year:
1. Keep just variables to be included in the analysis
2. Create a Year variable,
3. Rename, recode or reformat those variables, if needed.*/
/*2004 Dataset*/
data dbts_2004;
	set dbts04_2;
	rename 	psu= psunest
					riagendr= gender
					stratum= borostratum
					ageadj =  agegroup;
	key = put(sp_id, 7.);
	year = 2004;
	keep stratum psu HHNest	riagendr dbts key year dbts_wt_new wt_dbts ageadj;
	format riagendr gender. ageadj AGEGROUPFMT. dbts TE_3F. borostratum psunest key;
run;

/*2013-14 Dataset*/
data dbts_2014;
	set dbts14_2;
	year = 2014;
	keep borostratum psunest 	 HHNest gender dbts key year dbts_wt_new wt_dbts agegroup;
	format dbts TE_3F. borostratum psunest key;
run;

/*Step 2: Append the datasets	*/
/* 1. Verify names and formats*/
proc contents data=dbts_2004; run;
proc contents data=dbts_2014; run;

/* 2. Append the datasets*/
data dbts;
	set dbts_2004 dbts_2014;
	/*Restric datset to only cases to be included in the analysis*/
	if dbts =. then delete; 
run;
proc freq data=dbts;
	tables year*dbts/list missing;
run;

/* 3. Comparison Across Survey Years	*/
/*SUDAAN*/
proc descript data = dbts filetype = sas design = wr notsorted;
	nest borostratum psunest HHNest; /*Sampling design variables*/
	weight dbts_wt_new; /*Weight*/
	var dbts; /*Variable to be analyzed*/
	catlevel 1; /*State de category that refers to the positive outcome*/
	class  gender agegroup year;	/*State all variables involve in the analysis*/
	contrast year = (-1,1) / name = "2004 vs 2013-2014" ;/*Statement to compare the two survey years*/
	tables gender ;/*To get prevalence estimate by gender*/
	stdvar agegroup; ; /*Age standardization variable*/
	stdwgt 0.396579 0.371795 0.231626; ; /*standardizing proportions*/
	/*The following lines are all about appearance of the results*/
	rtitle "Total Diabetes (Age-Adjusted and weigthed): NYC-HANES 2004 vs 2013-2014";
	SETENV COLWIDTH=16 DECWIDTH=2 COLSPCE=2; 
	setenv decwidth = 1; 
	print NSUM="SAMSIZE"  percent sepercent t_pct p_pct / style = nchs
	nsumfmt = f8.0 percentfmt = f8.3 sepercentfmt = f8.2 p_pctfmt=f8.4;
run ;