# 9.2 - SAS - Poisson Regression Model for Count Data

Printer-friendly version

### Example - Horseshoe Crabs and Satellites

Here we work through this example in SAS.

This problem refers to data from a study of nesting horseshoe crabs (J. Brockmann, Ethology 1996); see also Agresti (2007) Sec. 3 and Agresti (2013) Sec. 4.3. Each female horseshoe crab in the study had a male crab attached to her in her nest. The study investigated factors that affect whether the female crab had any other males, called satellites, residing near her.  Explanatory variables that are thought to affect this included the female crab’s color (C), spine condition (S), weight (Wt), and carapace width (W). The response outcome for each female crab is her number of satellites (Sa). There are 173 females in this study. Datafile: crab.txt.

Let’s first see if the width of female's back can explain the number of satellites attached. We will start by fitting a Poisson regression model with only one predictor, width (W) via PROC GENMOD as shown in the first part of the crab.sas SAS Program as shown below:

Model Sa=w specifies the response (Sa) and predictor width (W).  Also, note that specification of Poisson distribution are dist=pois and link=log. The obstats option as before will give us a table of observed and predicted values and residuals. You can use any additional options in GENMOD, e.g., TYPE3, etc.

Here is the output that you should get from running just this first part of the crab.sas program (found in crab.lst):

Model Information

What do you learn from the "Model Information" section? How is this different from when we fitted logistic regression models? This section gives information on the GLM that's fitted. More specifically, we see that the response is distributed via Poisson, the link function is log and the dependent variable is Sa. The number of observations in the data set used is 173 and that's all of them were used in the analysis, that is there were no missing values neither for the response variable nor the predictor width.

Does the model fit well? What does the Value/DF tell you? Given the value of deviance statistic of 567.879 with 171 df, the p-value is zero and the Value/DF is much bigger than 1, so the model does not fit well. The lack of fit maybe due to missing data, covariates or overdispersion.

Is width the significant predictor? From the "Analysis of Parameter Estimates" table, with Chi-Square stats of 67.51 (1df), the p-value is 0.0001 and this rejects the null hypothesis that $\beta_W=0$, we can conclude that the width of a female's crab back is a significant predictor.

The estimated model is:  log (μi) = -3.3048 + 0.164Wi

The ASE of estimated β = 0.164 is 0.020 which is small, and the slope is statistically significant.

Interpretation: Since estimate of β > 0, the wider the female crab the greater expected number of male satellites on the multiplicative order as exp(0.1640) = 1.18. More specifically, for one unit increase in the width, the number of Sa will increase and it will be multiplied by 1.18.

If we look at the scatter plot of W vs. Sa (see further below) we may suspect some outliers, e.g., observations #48, #101 and #165. For example, #165 has W = 33.5, and Sa = 7. But by studying the residuals, we see that this is not an influential observation, e.g., standardized deviance residual StResdev = -0.739. You can consider other types of residuals, influence measures (like we saw in linear regression), as well as residual plots. Notice that there are some other points that have large outliers, e.g., #101. Here is the first two lines of the output and then it skips down to highlight the output for observations #48, #101 and #165.  For the entire listing of the output run the SAS code:

From the above output we can also see that the predicted counts ("Pred") and the values of the linear predictor ("Xbeta") that is the log of the expected counts. For example, for the first observation, Pred = 3.810, Xbeta = 1.3377, log(Pred) = Xbeta, that is log(3.810) = 1.3377, or exp(Xbeta) = Pred, that is exp(1.3377) = 3.810.

From the above analysis, we see that the predictor is significant but the model does not fit well. We can further assess the lack of fit by plotting residuals or influential points, but let us assume for now that we do not have any other covariates and try to adjust for overdispersion to see if we can improve the model fit.

#### Change the Model: Adjusting for Overdispersion

In the above model we detect a potential problem with overdispersion since the scale factor, e.g., Value/DF, is greater than 1.

What do you think overdispersion means for Poisson Regression? What does it tell you about the relationship between the mean and the variance of the Poisson distribution for the number of satellites? Recall that one of the reasons for overdispersion is heterogeneity where subjects within each covariate combination still differ greatly (i.e., even crabs with similar width will have different number of satellites). If that's the case, which assumption of the Poisson model that is Poisson regression model is violated?

As we saw in logistic regression, if we want to test and adjust for overdispersion we need to add the scale parameter by changing scale=none to scale=pearson; see the third part of the SAS program crab.sas labeled 'Adjust for overdispersion by "scale=pearson" '.  (As stated earlier we can also fit a negative binomial regression instead (also see the crab.sas code)

Below is the output when you run the second part of the program labeled 'Adjust for overdispersion by "scale=pearson" '.  How does this compare to the output above from the earlier stage of the code? Do we have a better fit now? Consider the "Scaled Deviance" and "Scaled Pearson chi-square" statistics? The overall model seems to fit better when we account for possible overdispersion with p-value about 0.33 for $G^2=178.454, df=171$.

With this model the random component does not have a Poisson distribution any more where the response has the same mean and variance. From the estimate given (e.g., Pearson X2 / 171 = 3.1822), the variance of random component (response, the number of satellites for each Width) is roughly three times the size of the mean.

The new standard errors (in comparison to the model where scale = 1), are larger, e.g., 0.0356 = 1.7839 × 0.02 which comes from the $scaled SE=\sqrt{Value/df}*original SE$. Thus the Wald X2 statistics will be smaller, e.g., 21.22 = 67.21 / 3.1822. Note that sqrt(3.1822) = 1.7839.

What could be another reason for poor fit besides overdispersion? How about missing other explanatory variables? Can we improve the fit by adding other variables?

Since we did not use the \$sign in the input statement to specify that the variable "c" was categorical, we can now do it by using class c as seen below. The following change is reflected in next section of the crab.sas program labeled 'Add one more variable as a predictor, "color" '. Again compare the parts of this output with the initial crab.lst output. Note the "Class level information" on 'color' indicating that this variable has 4 levels and thus are we are introducing "dummy variables" into the model. From the "Analysis of Parameter Estimates" output below we see that the reference level is level #4 (see the shaded output below). We are also adjusting for overdispersion but by using deviance instead of X2, although scale=pearson is preferred; we are doing this to demonstrate options in SAS and since the values are close, it doesn't matter which option we are using! The estimated model is:$\log{\hat{\mu_i}}$= -3.0974 + 0.1493W + 0.4474(C="1") + 0.2477(C="2") + 0.0110(C="3"). There does not seem to be a difference in the number of satellites between any color class and the reference level 4 according to the chi-squared statistics for each row in the table above. Furthermore, by Type 3 Analysis output below we see that the color is not statistically significant predictor after we take the width into consideration.  LR Statistics For Type 3 Analysis Chi- Source Num DF Den DF F Value Pr > F Square Pr > ChiSq W 1 168 14.96 0.0002 14.96 0.0001 C 3 168 0.85 0.4661 2.56 0.4640 Does this model fit the data better, with and without the adjusting for overdispersion? #### Change the Model: Include ‘color’ as a Numerical Predictor We will run another part of the crab.sas program labeled "Include color as a Numerical Predictor" by running the following change, Compare these partial parts of the output with the output above where we used color as a categorical predictor. We are doing this just to keep in mind that different coding of the same variable will give you different fits and estimates. What is the estimated model now?$\log{\hat{\mu_i}}$= -2.520 + 0.1496W - 0.1694C. Since adding a covariate doesn't help, the overdispersion seems to be due to heterogeneity. Is there something else we can do with this data? We can either (1) consider different methods, e.g., small area estimation, etc.. , (2) collapse over levels of explanatory variables, or (3) transform the variables. Let's consider grouping the data by the widths and then fitting a Poisson regression model for the rates. ### Change the Model: Model the Rate Data Here is the sorted data by W. The columns are in the following order: Widths, # Satellites, and Cumulative # of Satellites: Sorted  W Sa Cum Sa W Sa Cum Sa W Sa Cum Sa 1 21.0 0 0 25.3 2 103 27.3 1 270 2 22.0 0 0 25.4 6 109 27.4 5 275 3 22.5 0 0 25.4 4 113 27.4 6 281 4 22.5 1 1 25.4 0 113 27.4 3 284 5 22.5 4 5 25.5 0 113 27.5 6 290 6 22.9 4 9 25.5 0 113 27.5 9 299 7 22.9 0 9 25.5 0 113 27.5 1 300 8 22.9 0 9 25.6 0 113 27.5 6 306 9 23.0 1 10 25.6 7 120 27.5 0 306 10 23.0 0 10 25.7 8 128 27.5 3 309 11 23.1 0 10 25.7 5 133 27.6 4 313 12 23.1 0 10 25.7 0 133 27.7 6 319 13 23.1 0 10 25.7 0 133 27.7 5 324 14 23.2 4 14 25.7 0 133 27.8 0 324 15 23.4 0 14 25.7 0 133 27.8 3 327 16 23.5 0 14 25.8 10 143 27.9 7 334 17 23.7 0 14 25.8 0 143 27.9 6 340 18 23.7 0 14 25.8 0 143 28.0 0 340 19 23.7 0 14 25.8 0 143 28.0 1 341 20 23.8 0 14 25.8 3 146 28.0 4 345 21 23.8 0 14 25.8 0 146 28.2 6 351 22 23.8 6 20 25.8 0 146 28.2 8 359 23 23.9 2 22 25.9 4 150 28.2 11 370 24 24.0 0 22 26.0 4 154 28.2 1 371 25 24.0 10 32 26.0 3 157 28.3 8 379 26 24.1 0 32 26.0 14 171 28.3 15 394 27 24.2 0 32 26.0 9 180 28.3 0 394 28 24.2 2 34 26.0 5 185 28.4 3 397 29 24.3 0 34 26.0 3 188 28.4 5 402 30 24.3 0 34 26.1 5 193 28.5 0 402 31 24.5 5 39 26.1 3 196 28.5 1 403 32 24.5 1 40 26.2 0 196 28.5 9 412 33 24.5 1 41 26.2 3 199 28.5 3 415 34 24.5 6 47 26.2 3 202 28.7 0 415 35 24.5 0 47 26.2 0 202 28.7 3 418 36 24.5 0 47 26.2 0 202 28.9 4 422 37 24.5 0 47 26.2 0 202 29.0 1 423 38 24.7 0 47 26.2 2 204 29.0 4 427 39 24.7 5 52 26.2 2 206 29.0 10 437 40 24.7 0 52 26.3 1 207 29.0 3 440 41 24.7 4 56 26.5 1 208 29.0 1 441 42 24.7 4 60 26.5 4 212 29.0 1 442 43 24.8 0 60 26.5 0 212 29.3 4 446 44 24.9 0 60 26.5 0 212 29.3 12 458 45 24.9 6 66 26.5 4 216 29.5 4 462 46 24.9 0 66 26.5 7 223 29.7 5 467 47 25.0 3 69 26.7 5 228 29.8 4 471 48 25.0 2 71 26.7 2 230 30.0 8 479 49 25.0 8 79 26.7 0 230 30.0 9 488 50 25.0 5 84 26.8 5 235 30.0 5 493 51 25.0 6 90 26.8 0 235 30.2 2 495 52 25.0 4 94 26.8 0 235 30.3 3 498 53 25.1 5 99 27.0 3 238 30.5 3 501 54 25.1 0 99 27.0 3 241 31.7 4 505 55 25.2 1 100 27.0 6 247 31.9 2 507 56 25.2 1 101 27.0 6 253 33.5 7 514 57 27.0 0 253 58 27.1 8 261 59 27.1 0 261 60 27.2 5 266 61 27.2 3 269   Plot of Number of Satellites by Width of Crab—All Data Plot of Average Number of Satellites by Width of Crab—Distinct Widths Plot of Average Number of Satellites by Width— Widths Grouped #### To account for different widths, in this section we will group the Widths into intervals and re-analyze by using an OFFSET option in Model statement in SAS. The data have been grouped into 8 intervals, as shown in the (grouped) data below:  Interval NumCases AverWt AverSa SDSa VarSa =23.25 14 22.693 1.0000 1.6641 2.7692 23.25-24.25 14 23.843 1.4286 2.9798 8.8792 24.25-25.25 28 24.775 2.3929 2.5581 6.5439 25.25-26.25 39 25.838 2.6923 3.3729 11.3765 26.26-27.25 22 26.791 2.8636 2.6240 6.8854 27.25-28.25 24 27.738 3.8750 2.9681 8.8096 28.25-29.25 18 28.667 3.9444 4.1084 16.8790 >29.25 14 30.407 5.1429 2.8785 8.2858 Let’s run crabrate.sas (this is a new program not found in crab.sas). In this program we entered the grouped data above. In this case, each observation within a category is treated as if it has the same width. We also create a variable LCASES=log(CASES) which takes the log of the number of cases within each grouping. This is our adjustment value 't' in the model that represents the fixed space, in this case the group (crabs with similar width). We thus form a rate of satellites for each group by dividing by each group size, and are fitting a log-linear model to rate of satellites incidence given the crab's width. Here is the output. Note "Offset variable" under the "Model Information". Does the model now fit better or worse than before? It clearly fits better. For example the Value/DF for the deviance statistic now is 1.0861. The estimated model is:$log (\hat{\mu}_i/t)\$ = -3.535 + 0.1727widthi

As the width increases, the rate of satellites cases changes by exp(0.1727).

The residuals analysis indicate the good fit as well. From the table above we also see that the predicted number of counts ('Pred') cells correspond a bit better to the observed counts in the "SaTotal" cells. Let's compare the observed and fitted values in the plot below:

This last part is the output from crabrate.sas where just for demonstration we fit the Poisson regression model with the identity link for the rate data. Notice that this model does NOT fit well for the grouped data as the Value/DF for deviance statistic is about 11.649.