COUNTREG Procedure

(View the complete code for this example.)

This example demonstrates how to use the COUNTREG procedure to model count data with spatial effects. The data set Shigdata contains county-level data related to shigellosis, an infectious bacterial disease. The dependent variable Shigellosis records the number of shigellosis cases reported in each California county in 2011. The data are from the California Department of Public Health. Additional variables include PopDensity (population in thousands per square mile), Hispanic (percentage of Hispanic population), and BigHousehold (percentage of households with six or more persons); these data are from the 2010 United States Census. The following statements compute summary statistics for these variables and the frequency of each observed shigellosis count:

proc means data=Shigdata;
   var Shigellosis PopDensity Hispanic BigHousehold;
run;
proc freq data=Shigdata;
   table Shigellosis;
run;

The results are presented in Output 11.4.1.

The number of observations is 58, which is the number of California counties. The variance of Shigellosis (1653.99) is much larger than its mean (16.36), suggesting overdispersion. Moreover, the fact that 18 counties (31.03%) have counts of 0 indicates zero-inflation in the data.

A map of California counties that shows their shigellosis counts would help visualize any spatial effects. You can generate such a map by using the following statements:

ods graphics on;
goptions reset=all border;
data ca;
   set maps.counties;
   where state = 6;
run;

proc sort data=ca out=ca;
   by county;
run;

pattern value=mempty color=blue;

/* The COUNTY and COUNTIES data sets are unprojected */
/* Without the gproject procedure, you would get an inverse map */
proc gproject data=ca out=caproj;
   id state county;
run;

proc gmap map=caproj data=Shigdata all;
   id county;
   choro Shigellosis/discrete coutline=black;
run;

The map is shown in Output 11.4.2.

The following statements fit a Poisson model with spatially weighted regressors:

proc countreg data=Shigdata Wmat=W;
   model Shigellosis = PopDensity BigHousehold / dist=poisson;
   spatialeffects PopDensity BigHousehold;
run;

The SPATIALEFFECTS statement accounts for spatial effects on regressors in the MODEL statement. Because there might be local spillovers in the two explanatory variables, PopDensity and BigHousehold, both variables are included in the SPATIALEFFECTS statement. The data set W contains the spatial weights matrix. The model fit summary, convergence status, and parameter estimation results are shown in Output 11.4.3.

As shown in the "Parameter Estimates" table, all estimates are significant at the 5% level. The spatially weighted regressors are labeled with the prefix "W_". Because of the likely overdispersion in the data, you might consider using a negative binomial model. The following statements fit a negative binomial model with spatially weighted regressors:

proc countreg data=Shigdata Wmat=W;
   model Shigellosis = PopDensity BigHousehold / dist=negbin;
   spatialeffects PopDensity BigHousehold;
run;

The model fit summary, convergence status, and parameter estimation results are listed in Output 11.4.4.

The AIC and SBC values in the "Model Fit Summary" table in Output 11.4.4 are both smaller than those in Output 11.4.3, indicating that the negative binomial model provides a better fit than the Poisson model. The parameter estimate for _Alpha is an estimate of the dispersion parameter in the negative binomial distribution. The null hypothesis that _Alpha is 0 can be tested against the alternative hypothesis that _Alpha is positive, by using the likelihood ratio test. The likelihood ratio test statistic is , where and are the log likelihoods for the Poisson (null) and negative binomial (alternative) models, respectively. The likelihood ratio test is highly significant at the 5% level, providing strong evidence of overdispersion.

An alternative model to account for overdispersion is the Conway-Maxwell-Poisson (CMP) model, which uses dispersion regressors to model the parameter that controls the degree of dispersion. The following statements fit a CMP model that accounts for local spillovers in both model regressors and includes a DISPMODEL statement for the dispersion model:

proc countreg data=Shigdata Wmat=W;
   model Shigellosis =  PopDensity / dist=compoisson;
   spatialeffects PopDensity;
   dispmodel Shigellosis ~ Hispanic;
   spatialdispeffects Hispanic;
run;

The model fit summary, convergence status, and parameter estimation results are provided in Output 11.4.5.

The regression coefficients in the DISPMODEL statement are labeled with the prefix "Dsp_", and their spatially weighted counterparts in the DISPSPATIALEFFECTS statement are labeled with the prefix "Dsp_W". Both coefficients—Dsp_Hispanic and Dsp_W_Hispanic—are significant at the 5% level. The AIC and SBC values in Output 11.4.5 are both slightly higher than those from the previous negative binomial model. Therefore, the negative binomial model with spatially weighted regressors provides the best fit thus far.

Because 31% of observations have counts of 0, you might consider a zero-inflated negative binomial (ZINB) model. The following statements fit a ZINB model with local spillovers in the model regressors, with a ZEROMODEL statement to model zero inflation:

proc countreg data=Shigdata Wmat=W;
   model Shigellosis = PopDensity BigHousehold / dist=ZINB;
   zeromodel Shigellosis ~  Hispanic;
   spatialeffects PopDensity BigHousehold;
   spatialzeroeffects Hispanic;
run;

The model fit summary, convergence status, and parameter estimation results are shown in Output 11.4.6.

The "Parameter Estimates" table gives the estimate for the spatially weighted regressor from the ZEROMODEL statement, labeled with the prefix "Inf_W_". The estimates of regression coefficients for Inf_Hispanic and Inf_W_Hispanic are insignificant at the 10% level.