*************************************************************************************
* Program:  C:\NHANES\LogisticRegression_Uni.sas                                    *
* Proposal: Logistic regression analysis with one independent variable using SUDAAN * 
*************************************************************************************;

LIBNAME NH "C:\NHANES\DATA";
*LIBNAME NH "C:\myfiles";
OPTIONS NODATE NOCENTER;
options ls=78;

proc format;                                       
                                               
  VALUE sexfmt  1 = 'Male'                          
                2 = 'Female' 
                 ;   
  VALUE agefmt
  	            1='20-39 yrs'
	            2='40-59 yrs'
	            3='60 + yrs'; 
  VALUE bpfmt
  	            1='hypertension'
	            0='no hypertension'
                ;
  VALUE chfmt
  	            1='high cholesterol'
	            0='not high cholesterol'
                ;
  VALUE bmifmt
  	            1='BMI<25'
	            2='25<=BMI<30'
	            3='BMI>=30'
                ;
 run;

data analysis_data;
	set NH.analysis_data;

*create an indicator of hypertension based on blood pressure and medication;
	if (bpxsar >= 140 or bpxdar >= 90 or bpq050a = 1) then Hyper = 1;
  	else if (bpxsar ne . and bpxdar ne .) then Hyper = 0;

*create an indicator of high cholesterol based on total cholesterol and medication;
	if (lbxtc>=240 or bpq100d = 1) then HiChol = 1;
  	else if (lbxtc ne . ) then HiChol = 0;

	*create bmi groups;
	if 0<=bmxbmi<25 then bmigrp=1;
	else if 25<=bmxbmi<30 then bmigrp=2;
	else if bmxbmi>=30 then bmigrp=3;

	*create age groups;
	if 20<=ridageyr<40 then age=1;
	else if 40<=ridageyr<60 then age=2;
	else if ridageyr>=60 then age=3;

*create log of fasting triglycerides;
	logtrig=log(lbxtr);

if hyper ne . and hichol ne . and bmigrp ne . and age ne . and logtrig ne . and wtsaf4yr ne 0 then eligible=1;

	label 	Hyper='Indicates if have hypertension'
	        HiChol='Indicates if have High Cholesterol'
			BMIGRP='BMI groups, lean, overweight, obese'
			Riagendr = 'Gender'
			Age = 'Age Group' ;
run;

proc sort data=analysis_data;
	by sdmvstra sdmvpsu;
run;

proc crosstab 
data=analysis_data design=wr;
  nest sdmvstra sdmvpsu;
  weight wtmec4yr;
  subpopn eligible=1;
  class age riagendr hichol bmigrp hyper/NoFreq;
  table hyper*(age riagendr hichol bmigrp);
  print colper /
     nohead notime                                                      
     style=nchs ;                              
rformat riagendr sexfmt.;
rformat hyper bpfmt.;
rformat age agefmt.;
rformat hichol chfmt.;
rformat bmigrp bmifmt.;
rtitle "Cross tabulations for categorical variables: NHANES 1999-2002";
run;

%macro logit (var,fmt,ref);
proc rlogist data=analysis_data;
	nest sdmvstra sdmvpsu;
	weight wtmec4yr;
	*subpopn ridageyr>=20;
	subpopn eligible=1;
	class &var/ nofreq;
	reflevel  &var=&ref;
	model hyper=&var;
	test waldf satadjf satadjchi;
	rtitle "Univariate Association of &var and Odds of Hypertension";
	rformat &var &fmt;
    rformat hyper bpfmt.;
run;
%mend;
%logit(Age,agefmt.,2);
%logit(Riagendr,sexfmt.,2);
%logit(Hichol,chfmt.,1);
%logit(Bmigrp,bmifmt.,2);

proc rlogist data=analysis_data;
	nest sdmvstra sdmvpsu;
	weight wtsaf4yr;
	subpopn eligible=1;
	model hyper=logtrig;
	test waldf satadjf satadjchi;
	rtitle "Univariate Association of Log of Triglycerides and Odds of Hypertension";
	rformat hyper bpfmt.;
run;