Community Health Status Analysis

Step 01

Set up the CHSI library

The project uses the 2009 Community Health Status Indicators report compiled by the CHSI Project Working Group as part of the CDC. A macro sets the local data directory, LIBNAME creates the permanent CHSI library, and PROC CONTENTS inspects the demographics and risk-factor tables. The original program was written on May 28, 2017.

/*Using macro variables to upload data and create the permanent library,CHSI*/
%let path=C:\Users\ptl09\Downloads;
%put &path;
libname CHSI "&path";

/*checking to make sure everything's there*/
proc contents data=chsi.demographics;
run;

proc contents data=chsi.riskfactors;
run;

Step 02

Prepare the variables and response

The demographics table is reduced to population size, population density, poverty, and five racial or ethnic population percentages. The risk-factor table supplies the HPSA indicator and a community health center indicator. The program merges on State_FIPS_Code and defines shortage = 1 when HPSA_Ind = 2; all other values are assigned 0. Negative poverty and population-density values are removed. In SAS, numeric missing values also compare below zero, so these filters remove missing values in those two fields—not every field.

Reproducibility note: the original merge uses the state identifier alone and includes no sort step. Before rerunning county-level data, verify sorting and whether a county identifier is needed to align records correctly. Missing HPSA values also fall into the “otherwise 0” branch.

/*Creating smaller data sets with only the needed variables*/
data demographics1;
	set chsi.demographics;
	keep 'State_FIPS_Code'n 'Population_Size'n 'Population_Density'n Poverty White 
		Black 'Native_American'n Asian Hispanic;
run;

data riskfactors1;
	set chsi.riskfactors;
	keep 'State_FIPS_Code'n 'Community_Health_Center_Ind'n 'HPSA_Ind'n;
run;

/*Merging the two data sets by State_FIPS_Code*/
data health;
	merge demographics1 riskfactors1;
	by 'State_FIPS_Code'n;
run;

/*Setting up an indicator variable for my binary response and deleting all missing observations*/
Data health1;
	set health;

	if 'HPSA_Ind'n=2 then
		shortage=1;
	Else
		shortage=0;

	If poverty < 0 then
		delete;

	If 'Population_Density'n < 0 then
		delete;
Run;

Step 03

Explore relationships with LOESS

Two panels of scatterplots use LOESS curves to explore how shortage varies with community characteristics. My source comments describe increasing shortage with population size and density, a decreasing trend with poverty, and mixed or nonlinear patterns across demographic percentages. State FIPS is explored here but is not included in the later regression models. These are exploratory associations, not causal conclusions.

/*Using a panel display of scatterplots with a loess curve to see which
variables might be a good fit for logistic regression*/
proc sgscatter data=health1;
	compare y=shortage x=('State_FIPS_Code'n 'Population_Size'n 
		'Population_Density'n poverty)/loess;
run;

proc sgscatter data=health1;
	compare y=shortage x=(White Black 'Native_American'n Asian Hispanic)/loess;
run;

Step 04

Fit a poverty-only logistic model

PROC LOGISTIC models the probability of shortage with poverty as the sole predictor. A second fit requests exponentiated coefficients and pseudo-R². The original comments record an intercept of 2.76 and a poverty coefficient of −0.1154. These describe log-odds, not a linear probability: log[p / (1 − p)] = 2.76 − 0.1154 × Poverty. The script does not include the resulting pseudo-R² value.

/*Running a logistical regression on health data
with health care profession shortage as the response and poverty as the predictor*/
proc logistic data=health1;
	model shortage(event='1')=poverty;
run;

/*Intercept=2.76  Poverty=-0.1154 or Shortage=2.76-0.1154Poverty*/
/* Produce the psuedo r^2 to see how much of the variation in health care shortages
is accounted for by poverty level */
proc logistic data=health1;
	model shortage(event='1')=poverty/expb rsquare;
run;

Step 05

Check linearity with an empirical logit

PROC RANK divides poverty into four groups. PROC MEANS calculates the number of shortage events and mean poverty in each bin. The empirical logit uses a 0.5 continuity correction before PROC SGPLOT fits a line. My original interpretation describes a negatively linear pattern, supporting a linear poverty term on the log-odds scale.

/*Create a empirical Logit plot to determine linearity*/
proc rank data=health1 groups=4 out=out;
	/*4 groups in the
	output file out*/
	var Poverty;
	ranks bin;

	/*Create the variable bin for grouping*/
run;

/*Summarize the data by each bin*/
proc means data=out nway noprint;
	class bin;
	var shortage Poverty;
	output out=bins sum(shortage)=Shortage mean(Poverty)=Poverty;
Run;

/* Empirical Logit Plot*/
data bins1;
	set bins;
	logit=log((shortage+.5)/(_FREQ_ - shortage + .5));
run;

Title 'Empirical Logit Plot';

proc sgplot data=bins1;
	reg x=Poverty y=Logit;
run;

Title;

/* The model appears to be negatively linear*/

Step 06

Review diagnostics and predictor correlation

The poverty model requests diagnostic index plots using IPLOTS. My source notes report no significant outliers. PROC CORR then examines all eight candidate predictors; the comments flag a negative association between Black and White population percentages without identifying other multicollinearity concerns. These are the interpretations recorded in the program; diagnostic plots and correlation tables are not included in the linked file.

/*Performing a residual analysis to determine if there are any outliers*/
proc logistic data=health1;
	model shortage(event='1')=Poverty/iplots;
run;

/*no significant outliers*/
/*Use a correlation analysis to determine whether there are any
multicollinearity problems with all the predictors included*/
proc corr data=health1;
	var 'Population_Size'n 'Population_Density'n Poverty White Black 
		'Native_American'n Asian Hispanic;
run;

/*Black and White are negatively correlated as expected but nothing else indicates a multicollinearity problem*/

Step 07

Select a reduced model

Backward elimination begins with eight candidate predictors: population size, population density, poverty, White, Black, Native American, Asian, and Hispanic percentages. My source comments identify a six-predictor model retaining population size, population density, poverty, Black, Native American, and Hispanic percentages. White and Asian percentages are not retained. The selection statistics and coefficient table are not present in the source file.

/* Variable selection in Logistic Regression using backwards elimination*/
proc logistic data=health1;
	model shortage (event='1')='Population_Size'n 'Population_Density'n Poverty 
		White Black 'Native_American'n Asian Hispanic/selection=backward;

	/*all sebsets selection*/
run;

/* The best model includes population size, population density, Poverty, Black, Native_American, Hispanic*/

Step 08

Test the omitted predictors

The full model is refitted and TEST Asian, white checks whether the two coefficients are jointly zero. Although the original comment calls this a nested F test, PROC LOGISTIC TEST performs a joint Wald chi-square test, not an F test or likelihood-ratio test. My source interpretation says these variables did not contribute to the model; no test statistic or p-value is included.

/* Produce a nested F test to compare the full model with the reduced model*/
proc logistic data=health1;
	model shortage(event='1')='Population_Size'n 'Population_Density'n Poverty 
		White Black 'Native_American'n Asian Hispanic;
	test Asian, white;
run;

/*'Asian' and 'White' had no effect on the model*/

Step 09

Compare the original ROC models

The final procedure requests ROC curves for the full eight-predictor model, the six-predictor model labeled “best,” and the poverty-only model labeled “single.” ROCCONTRAST uses “best” as the reference and requests estimates. The GitHub file contains the SAS instructions, but not the rendered curves, AUC values, or comparison results. This walkthrough preserves that code without reconstructing or simulating its outputs.

/*Compare ROC curves for the ful model, 'reduced' best model, and the single 'Poverty' model*/
proc logistic data=health1 plots=ROC;
	model shortage (event='1')='Population_Size'n 'Population_Density'n Poverty 
		White Black 'Native_American'n Asian Hispanic;
	ROC 'best' 'Population_Size'n 'Population_Density'n Poverty Black 
		'Native_American'n Hispanic;
	ROC 'single' Poverty;
	ROCcontrast reference('best')/estimate e;
run;