/*Clinic Only Weight Adjusted*/

/*Program highlights:
1. First uses the Clinic Only weight because the ACASI component was only conducted at the 
clinic.
2. An Adjustment Weight for the Clinic Only Weight was then created due to missing data in the
subsetted sample.
2. Calculates Crude prevalence estimates of ever had sex and corresponding Relative 
Standard Errors (RSEs).
3. Calculates Age-Adjusted prevalence estimates of ever had sex and corresponding RSEs. 
*/


proc format;

/*gender*/
value riagendr
	1='Male'
	2='Female';

/*Race_eth*/
value race_eth
	1='Non-Hisp White'
	2='Non-Hisp Black'
	3='Non-Hisp Asian' /*(Non-Hispanic Asian and Hawaiian/PI)*/
	4='Hispanic'
	5='Non-Hisp Other' /*(American Indian, Other, Multiracial who could not select one race, and
	SPs missing race but Non-Hispanic)*/;


libname hanes 'C:\Data';

proc sort data=hanes.spfile out=spfile; by sp_id;
proc sort data=hanes.capi out=capi; by sp_id;
proc sort data=hanes.acasi out=acasi; by sp_id;

data acasi_ex; merge spfile (in=a) capi acasi (in=b); by sp_id;
if a and b; /*Keeps SPs who are in both ids.spcomplete AND ACASI*/

/*Ever had sex*/
if sxq020=1 then sex=1;
else if sxq020=2 then sex=2;
else delete; /*Limits analysis to SPs who were not missing/Ref/DK for this 
	initial question to Sexual Behavior section*/
run;

/*Merge in Adjustment Weights for 'CLINIC ONLY' weight - Refer to Analytic Guidelines for
	further details.*/

/*Sum Clinic Only weights of the respondents (in subsetted dataset) within each category of age,
gender, race/ethnicity - 'SUM B'*/
proc summary data=acasi_ex nway;
var wtsf1c; /*Clinic Only Weight*/
class agewt riagendr racewt;
output out=samplewts sum=samplewt;
run;

/*The file IDS.clin_con_tots contains the sum of the Clinic Only Control Totals within each
category of age, gender, race/ethnicity (from Westat) - 'SUM A'*/
/*Merge SUM B to SUM A to create Weight Factor (A/B)*/
proc sort data=samplewts; by agewt riagendr racewt;
proc sort data=hanes.clin_con_tots out=clinicwts; by agewt riagendr racewt;
data wtfactor; merge samplewts clinicwts; by agewt riagendr racewt; 
wtfactor=clinic_control_total/samplewt;
run;

/*Merge dataset with new weight factor to subsetted dataset in order to adjust the weight by the
weight factor A/B for each SP in the cell*/
proc sort data=wtfactor; by agewt riagendr racewt;
proc sort data=acasi_ex; by agewt riagendr racewt; run;
data ac_adwt; merge acasi_ex (in=a) wtfactor; by agewt riagendr racewt; if a;
acwt = wtsf1c * wtfactor; /*New weight = Old weight * (A/B)*/
run;

/*Check adjusted weight*/
proc summary data=ac_adwt nway;
var acwt;
output out=popchk sum=pop;
run;
proc print data=popchk; run; /*5,827,719*/


/*SUDAAN CODE*/

/*Ever had sex - CRUDE (using PROC CROSSTAB)*/
proc sort data=ac_adwt; by stratum psu; run;
PROC CROSSTAB data=ac_adwt;
NEST STRATUM PSU;
WEIGHT acwt; /*Adjusted Clinic Only Weight*/
subgroup sex riagendr agegroup race_eth;
LEVELS    2     2     	 5		  5;
TABLES  RIAGENDR*sex;
RTITLE "Ever had Sex Among NYC Population, by Gender - Crude";
PRINT NSUM="SAMSIZE" WSUM="POPSIZE" ROWPER SEROW lowrow uprow/
TESTS=ALL WSUMFMT=F9.0 SEROWFMT= F6.3 style=nchs;
output/filename=sexcrude filetype=sas TABLECELL=DEFAULT REPLACE;
SETENV COLWIDTH=14 DECWIDTH=2 COLSPCE=2; run;

/*Calculate Relative Standard Error of the estimates using the SUDAAN output dataset.  Typically,
estimates with RSE >=30% are considered unreliable and notated that it should be interpreted with
caution.*/
data rse; set sexcrude;
RSE=round(serow/rowper*100,0.1);
run;


/*Ever had sex - Age-Adjusted (using PROC DESCRIPT)*/
proc sort data=ac_adwt; by stratum psu; run;
PROC descript data=ac_adwt;
NEST STRATUM PSU;
WEIGHT acwt; /*Adjusted Clinic Only Weight*/
var sex sex;
catlevel 1 2;
subgroup sex riagendr ageadj race_eth;
LEVELS    2     2        3		  5;
TABLES  RIAGENDR;
stdvar ageadj;
stdwgt 0.396579 0.371795 0.231626;
RTITLE "Ever had sex Among NYC Population, by Gender - Age-Adjusted";
PRINT NSUM="SAMSIZE" WSUM="POPSIZE" total percent sepercent lowpct uppct/ 
WSUMFMT=F8.0 totalfmt=f8.0 percentfmt=f8.1 lowpctfmt=f8.1 uppctfmt=f8.1 
SEpercentFMT=F3.1 style=nchs;
output/filename=sexagead filetype=sas TABLECELL=DEFAULT percentfmt=f8.3 lowpctfmt=f8.3
uppctfmt=f8.3 REPLACE;
SETENV COLWIDTH=14 DECWIDTH=2 COLSPCE=2; run;

/*Calculate Relative Standard Error of the estimates using the SUDAAN output dataset.  Typically,
estimates with RSE >=30% are considered unreliable and notated that it should be interpreted with
caution.*/
data rse; set sexagead;
RSE=round(sepercent/percent*100,0.1);
run;