Correspondence to Gilbert Berdine MD
Email: Gilbert.Berdine@ttuhsc.edu

SWRCCC 2015;3(11):55-59
doi:10.12746/swrccc2015.0310.149

Previously, we have introduced Poisson and Negative binomial regressions for modeling count data. Here we will use a real example to demonstrate how to use SAS software performing such analyses.

A new oral antibiotic drug Gorilacillin was developed and has had excellent effects in several clinical trials. Gorilacillin has two side effects, including rash and elevated liver function tests (ELFT). To evaluate whether different patient groups have different risks of having such side effects, a Gorilacillin side effect study was conducted. A total of 5,275 participants were recruited from five participating countries. All the patients were followed for up to one month, and the number of patients who developed rash or ELFT was recorded. The goal of the study was to investigate whether there was a significant difference in developing side effects for patients in different age groups.

Data from the study were collected and stored in an Excel table (see below for a partial view of the table).

It is typical to use a Poisson or negative binomial regression for analyzing such data since the outcome is a side effect count, and the probability for developing side effects is low – rare events. Meanwhile, because rash and ELFT are the two main side effects of Gorilacillin, an event is declared if a patient develops either rash or ELFT. Note that the numbers of patients in different age groups/countries are different; thus it makes sense to model the rate of side effect per patient, as a function of age group and country, to adjust for the differences in number of patients.

A scatter plot can usually help visualize potential relationships. As we can see from Figure 1, there does not seem to have a strong relationship between side effect rate and age group. Also the side effect rates are quite similar among countries.

To apply a Poisson regression, we have,

To adjust for the numbers of patients in different age groups/countries, we use the rate of side effect (dividing the expected number of events by the number of patients) as the outcome variable. Equivalently, the above equation can also be written as,

,

where the additional term on the right-hand side, log(n), is called an offset. Corresponding to the above two models, there are two equivalent SAS statements:

proc genmod data=data;

class country age;

model incident/n = age country / dist= poisson link=log;

lsmeans age / ilink diff cl;

run;

Or equivalently,

proc genmod data=data;

class country age;

model incident = age country / offset=logn dist=poisson link=log;

lsmeans age / ilink diff cl;

run;

Note that, in the second equation, log⁡n=log⁡(n), where n is the number of patients in a specific group. The lsmeans statement can be used to obtain the side effect rate estimates for the 10 age groups, averaged over countries. The ilink option specifies the inverse link function to be used for calculating the rate estimates, and the cl option produces the confidence intervals. In addition, the diff option provides all pairwise comparisons of side effect rates among age groups.

The above SAS output table shows that age group was not significantly associated with Gorilacillin side effect rate and there was no significant difference among the 5 countries.

From the lsmeans estimates, we see that the estimated side effect rate for patents 0-4 years old was 1.1% (the Mean column; table above) with a confidence interval (0.6%, 2.2%; the Lower Mean and Upper Mean columns), and for the 5-9 years old it was 1.9% with a confidence interval (1.1%, 3.3%), etc.

The diff option does provide all pairwise comparisons should such comparisons be of interest (table below shows part of the comparisons).

Now, recall that we previously explained that a negative binomial regression model might be more appropriate should data overdispersion exist. To test overdispersion, an easy way is to apply a negative binomial regression with scale=0 and noscale options in the model statement. These options test whether overdispersion of the form μ+kμ^2 exists by testing whether the dispersions parameter equals to 0.

proc genmod data=data;

class country age;

model incident/n = age country / dist= nb link=log scale=0 noscale;

run;

From the above test of overdispersion result, we can see that the p value is less than 0.0001, and thus it is appropriate to use a negative binomial regression.

The result from the negative binomial regression (table above) is similar to that from the Poisson regression. We did not detect any difference in side effect rate between the reference and other age groups. Looking at the raw data in the scatter plot (Figure 1), one might think that Americans of age 25-29 were at risk for adverse effects of the drug, but the statistical analysis shows the result to be within the 95% confidence limit for purely random effect compared to the reference group. The American 25-29 data point appears, at first glance, to be an outlier with some non-random effect, but, in fact, it is a purely random walk from the other data points.

The statistical analysis is consistent with the reality of the situation. Gorilacillin does not exist and the data was simulated by sampling rare occurrence events from an online game in which the game developers assure us that the events are, indeed, random. The game has generated all sorts of “theories” about how to elicit these rare events more often, but the statistical analysis shows the “theories” to be no more substantial than Americans age 25-29. This example illustrates how rare events can seem to generate “outliers” that are merely results of small samples and rare occurrence rates.

Many times, rare events are hard to observe, and it might take quite some time before one event is observed. If feasible, one alternative strategy of studying an association between a rare event and potential risk factors is to collect data retrospectively. For example, identify the list of patients who had the event, match them with those who did not have the event, then collect all the necessary data and perform data analysis.

...................................................................................................................................................................................................................................................................................................................................

Return to top