PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
 
Biometrics. Author manuscript; available in PMC 2010 July 2.
Published in final edited form as:
PMCID: PMC2889169
NIHMSID: NIHMS117728

Covariate Adjustment and Ranking Methods to Identify Regions with High and Low Mortality Rates

SUMMARY

Identifying regions with the highest and lowest mortality rates and producing the corresponding color-coded maps help epidemiologists identify promising areas for analytic etiologic studies. Based on a two-stage Poisson-Gamma model with covariates, we use information on known risk factors, such as smoking prevalence, to adjust mortality rates and reveal residual variation in relative risks that may reflect previously masked etiologic associations. In addition to covariate adjustment, we study rankings based on standardized mortality ratios, empirical Bayes estimates, and a posterior percentile ranking method and indicate circumstances that warrant the more complex procedures in order to obtain a high probability of correctly classifying the regions with the upper 100γ% and lower 100γ% of relative risks for γ= 0.05, 0.1 and 0.2. We also give analytic approximations to the probabilities of correctly classifying regions in the upper 100γ% of relative risks for these three ranking methods. Using data on mortality from heart disease, we found that adjustment for smoking prevalence has an important impact on which regions are classified as high and low risk. With such a common disease, all three ranking methods performed comparably. However, for diseases with smaller event counts, such as cancers, and wide variation in event counts among regions, empirical Bayes and posterior percentile ranking methods outperform ranking based on standardized mortality ratios.

Keywords: Incidence rate maps, Mortality rate maps, Poisson-Gamma model, Ranking procedures, Standardized incidence ratio, Standardized mortality ratio

1. Introduction

Epidemiologists study maps of disease incidence rates and mortality rates to find clues to disease etiology. Regions with the highest and lowest rates are of particular interest, because analytic epidemiologic studies can be conducted in such regions to identify etiologic factors. Standardized incidence ratios (SIRs) or standardized mortality ratios (SMRs), namely the ratio of the number of deaths, Oi, in region i divided by the expected deaths, Ei, are often ranked, and the corresponding regions are color coded according to decile or other quantile of the rank. One often focuses on the regions with the highest risk (Devesa et al., 1999).

Two important issues arise in identifying regions with the highest or lowest mortality rates. (Hereafter we refer to mortality rates, but the same ideas apply to incidence rates). One issue is whether the variation in SMRs is partly explained by known factors such as smoking prevalence. One would like to adjust for such factors to identify regions with the highest or lowest residual mortality rates where one might find new etiologic factors apart from e.g. smoking. A second issue is what procedures to use for ranking. The SMR is unstable in areas with small values of Ei, and a two-stage Poisson-Gamma model (Clayton and Kaldor, 1987; Pennello, Devesa, and Gail, 1999) can be adapted to produce smoother Bayes or empirical Bayes (EB) estimates for ranking the regions. Although EB estimates produce good estimates of the ranks overall (Laird and Louis, 1989), they are not optimal for identifying the regions with the upper 100γ% of rates (Lin et al., 2006). For that purpose, Lin et al. (2006) recommended ranking the posterior probability that a region’s rank will be in the top 100γ% of ranks, which we call fully Bayesian posterior percentile ranking (B-PPR). In this paper, we propose an empirically approximated posterior percentile ranking procedure (PPR) that implicitly assumes vague priors on certain parameters, and we compare PPR with B-PPR. We use the abbreviation PPR for ranking either in the upper 100γ%of ranks or lower 100γ% of ranks and the abbreviations PPHR and PPLR for the upper 100γ% and lower 100γ% of ranks respectively.

We describe methods for covariate adjustment for the SMR and for the EB and PPR approaches under the Poisson-Gamma model (Section 2) and compare these procedures in an analysis of heart disease mortality data (Section 3). In this example, we find that covariate adjustment has a greater impact than the choice of ranking method, which motivates simulations and analytic studies to determine under what circumstance the choice of ranking method has a material impact (Section 4). The choice of ranking method is only important when there is wide variation in Ei. In this circumstance, PPR outperforms EB slightly, and both can be noticeably better than SMR if the expected variance of the SMR is large. We present accurate analytic approximations to guide the choice of ranking procedure, both for the Poisson-Gamma and Gaussian-Gaussian models (Section 4). We conclude with a discussion in Section 5.

2. Notation and Methods for Covariate Adjustment and Ranking

Consider I regions, i = 1, 2,(...), I with J age-groups j = 1,(...),J in each. Let yij be the number of person-years exposure and Oij be the observed deaths in region i and age-group j. Note that j might also index a cross-classification on age, race and gender. Suppose we know the smoking prevalence Pij in region i and age-group j. Then the expected numbers of deaths are modeled as

Eij=yijξj{Pijr+(1Pij)},
(1)

where {ξj}, the age-specific mortality rates for non-smokers, and r, the relative risk associated with smoking, are to be estimated.

We would like to know the unmeasured random effect θi which is assumed to come from a gamma distribution Γ(α,1/α) with mean α(1/α)=1 and variance 1/α, where α is unknown. The quantity θi represents the intrinsic relative risk (or frailty) of region i apart from the measured factors encompassed in Eij. Ideally we would be able to rank the {θi} and identify those regions with ranks in the upper 100γ% of ranks and those regions with ranks in the bottom 100γ% of ranks. Conditional on θi, Oij have the Poisson distribution with mean Eijθi,

Oij|θi~Poi(Eijθi).
(2)

We assume θi is independent of Eij. Thus, unconditionally, E(Oij)=EijE(θi)=Eij. This model for covariates emphasizes that we wish to identify regions with high and low intrinsic θi apart from measured factors such as age and smoking prevalence. Clayton and Kaldor (1987) and Christiansen and Morris (1997) incorporate covariates into the mean of θi. If they had used the mixture model (1), their methods would yield the same marginal distribution of Oij as we obtained, namely the negative binomial,

Oij~NB{success=α,prob=α/(α+Eij)},
(3)

with mean = success(1-peob)/prob and variance = success(1-prob)/prob2. The marginal likelihood is the product over regions and age groups of (3), which can be maximized to yield maximum likelihood estimates [psi]={[alpha], [r with circumflex] and [Xi w/ hat]j=exp([delta with circumflex]j)}. The parameterization in terms of δj avoids numerical instabilities. Assuming ψ={α, r and ξj} are known, the posterior distribution of θi given Oi=∑jOij is

θi|Oi~iidΓ{Oi+α,1/(Ei+α)}.
(4)

Thus the EB estimate of θi is

θ^iEB=(Oi+α^)/(E^i+α^),
(5)

where Êi = ∑jÊij, and Êij is obtained from (1) with [r with circumflex] and [Xi w/ hat]j in place of r and ξj.

The fully Bayesian posterior percentile ranking (B-PPR) method was first proposed in Lin et al. (2006). They assumed a Poisson-Gamma model with vague and independent priors on the unknown parameters. To implement their procedures, we use α ~ LN(3.69, 100), r ~ LN(0.69, 100), and ξ~iidLN(6,100). Here LN(η, [var phi]2) is the lognormal density of a random variable whose logarithm has mean η and variance [var phi]2. The top 100γ% highest risk regions are identified by ranking the following posterior probabilities,

P(rank(θi)(1γ)I|O).
(6)

These posterior probabilities can be estimated by sampling from the joint posterior distribution of {θi} using a Gibbs or Metropolis-Hastings algorithm. In this study, the B-PPR analysis was conducted in WinBugs1.4.3 (Lunn et al.; 2000). Based on graphical assessment of convergence, the following procedure was used for each simulation. We discarded the initial 2500 samples and used every fifth sample thereafter to obtain 1000 samples of parameters.

In this paper, we propose an approximate posterior percentile ranking (PPR) method and compare it with B-PPR. For PPR, we approximate the joint posterior distribution of ψ={α, r and ξj} as a multi-normal distribution with means [psi] and with covariance given by the inverse of the observed Fisher information matrix from the marginal likelihood based on (3), and we obtain an approximate posterior sample of ψ, namely ψ* from this multi-normal distribution. Given ψ*, the approximate posterior {θi*} are conditionally independent and can be independently generated from (4). In order to approximate (6), this sample of {θi*} is ranked and the regions with ranks in the top 100γ% of ranks (and bottom 100γ% of ranks) are recorded. We repeated this approximated posterior sampling procedure N = 1000 times and selected as “high” the γI regions that most often had among the top 100γ% of ranks. Likewise we determined which regions were “low” as those that most often had among the lowest 100γ% of ranks. To be precise, let {Si} be a sample from the posterior distribution of {θi*} and for a particular region i consider P{rank(Si) > (1−γ)I}. We estimate this probability as above and rank the regions according to PPHRi = rank{P(rank Si > (1−γ)I)}. The “high” regions are those with PPHRi > (1 − γ)I, and the “low” regions are those with PPLRi > (1 − γ)I, where now PPLRi = rank{P(rank Si ≤ γ I)}.

We compare our PPR approximation with B-PPR using the heart disease example described in the following section, and we find that PPR performs as well as or better than B-PPR in ranking regions. Moreover, PPR takes much less computing time than B-PPR.

3. Example of Mortality from Heart Disease

We studied mortality from heart disease (ICD codes I00-I09, I13, I20-I51 based on ICD-10, 1992) in white U.S. males for the period 1990–2005. The data were presented separately in 805 health service areas (HSAs), which are defined as a single county or cluster of contiguous counties that are relatively self-contained with respect to hospital care. Heart disease death counts and person-years for each HSA and 13 age intervals (25–29, 30–34, …, 80–84, 85+ years) were obtained from the National Cancer Institute’s Surveillance Epidemiology and End Results (SEER) Program data at http://seer.cancer.gov.

We estimated the prevalence of current smoking {Pij} from aggregated 1995–2005 data from the Behavioral Risk Factor Surveillance System (BRFSS) (http://www.cdc.gov/brfss/technical_infodata/surveydata/1995.htm). Four HSAs in Alaska were eliminated because smoking data were too sparse. Twenty HSAs had total sample sizes less than 40; to estimate {Pij} for these HSAs, we aggregated their data with that of adjacent counties within the same state. For some HSAs, there were too few samples in particular age groups to estimate {Pij} reliably. Therefore, in all HSAs we aggregated data into three age groups (25–44, 45–64, and 65+ years) to estimate {Pij}. Thus, several values of j share the same estimate Pij. In the analysis we assume that the estimated {Pij} are known without error (see Section 5).

In this example, we use the ranking procedures presented in the last section to identify regions in the top 10% and bottom 10% of relative risks. To highlight the “high” and “low” HSAs, we plotted those in the top decile in black and those in the bottom decile in grey in Figure 1a, 1b, 1c and 1d. Other areas are shown in white, except the four excluded HSAs in Alaska, which are light grey.

Figure 1
Rankings of heart disease mortality relative risks in 805 Health Service Areas (HSAs) for white men, 1990–2005. Those HSAs in the top decile are shown in black; those in the bottom decile are shown in grey; others are in white. Four HSAs in Alaska ...

Figure 1a depicts SMRs that are not adjusted for smoking. Adjustment for smoking (Figure 1b) produces noticeable changes in the pattern of “high” risk HSAs especially in and above Texas. Likewise areas of “low” risk are changed in California and elsewhere. Of those 81 HSAs classified as “high” risk in Figure 1b, only 60 were classified as high without adjustment for smoking (Figure 1a). Of those 80 classified as “low” risk in Figure 1b, only 62 were so classified in Figure 1a. Thus, adjustment for smoking changes the identification of high and low risk HSAs appreciably.

The choice of ranking procedure had less impact than adjustment for smoking. Figure 1b, 1c and 1d depict HSA deciles of risk based on the SMR, EB and PPR methods of ranking respectively, all adjusted for smoking. Of the 81 “high” risk HSAs identified by the PPR method, 80 were also ranked “high” by the EB method, and 73 by the SMR method. Likewise, of the 80 HSAs ranked “low” by the PPR method, 79 were also ranked low by the EB method and 70 by the SMR method. Thus the choice of ranking procedure had relatively little impact.

To study factors that might influence the relative performance of the ranking procedures and compare the performance of B-PPR and PPR methods, we conducted simulations based on the heart disease data. Using the estimates from the original data, [alpha]=42.43, [r with circumflex]= 1.71, and [Xi w/ hat]=(0.281, 0.487, 0.987, 2.17, 4.65, 9.22, 15.94, 26.15, 43.01, 73.17, 114.48, 181.57, 303.08, 622.98)×10−4, and yij and Pij, we simulated 1000 data sets. In these simulations we knew that “true” ranking of the {θi} and could therefore calculate the proportion of the truly high and truly low {θi} that were correctly classified by the ranking procedures, with γ = 0.05, 0.10, and 0.20. Heart disease mortality rates are higher than rates for most diseases including cancers. We therefore divided yij by 10 and repeated the simulations.

Table 1 reports results for ranking procedures for smoking adjusted B-PPR, PPR, EB and SMR, and for SMR without smoking adjustment. Across different values of γ and yij, PPR performed slightly better than B-PPR. The B-PPR method took about 7 minutes per simulation, whereas PPR took less than one minute. We found that a slight change in the mean of the vague prior could change the ranking results for B-PPR slightly for a given data set. Because PPR was faster than B-PPR, performed at least as well as B-PPR and does not require special MCMC expertise, we focus on PPR hereafter.

Table 1
Mean percentages (with standard error) of correctly classified HSAs in simulated cardiovascular mortality dataa.

In simulations reflecting the original heart disease mortality data (columns 2 and 3 in Table 1), the proportions of HSAs correctly classified as low or high risk were lower without smoking adjustment. The PPR and EB procedures performed comparably and very slightly better than SMR with adjustment for smoking. The advantage of the PPR and EB procedures was more visible, however, in the scenario in which person-years were reduced by a factor of ten. In this setting, PPR performed slightly better than EB, and both outperformed SMR noticeably. The performance of all the procedures was less satisfactory than with the original yij, because the variance of SMRi was larger relative to the true dispersion of the {θi}.

4. Simulations and Analytical Studies to Compare Procedures for Selecting High Risk Groups

The example in Section 3 indicates that the proportion of truly high (or low) regions correctly detected decreases as the expected conditional variance of SMRi which is E(Eiθi)/Ei2=1/Ei, increases. Unreported simulations as in Section 3 indicate that the superiority of PPR over EB increases with increasing variability of 1/Ei across regions. These findings are analogous to results in Table 2 of Lin et al. (2006) for the Gaussian-Gaussian model. We now explore these relationships more systematically for the Gaussian-Gaussian and Poisson-Gamma models and provide excellent analytic approximations for the probability of correct classification of regions in the upper 10% of risk for the MLE (or SMR), EB, and PPHR ranking procedures. In what follows, we ignore the effects of smoking or other covariates.

4.1 Gaussian-Gaussian Model

The Gaussian mean is normally distributed μi ~ N(μ0, τ2). Given μi and the conditional variance σi2, we observe Xi~N(μi,σi2). The conditional variance σi2 might depend on the sample size in a survey in region i, for example. We assume that log(σi2)~N(η,[var phi]2) is independent of μi. We assume μ0, τ2, η and [var phi]2 are known. The joint density is

f(Xi,σi2,μi)=ϕ(μi;μ0,τ2)LN(σi2;η,[var phi]2)ϕ(Xi;μi,σi2),
(7)

where ϕ(μi; μ0, τ2) is a normal density with mean μ0 and variance τ2 and LN(σi2;η,[var phi]2) is the log-normal density of a random variable whose logarithm has mean η and variance [var phi]2. In a sample of I independent regions, we wish to identify those regions whose values μi have ranks [0.9I], [0.9I]+1, ……, I, where [ ] is the integer part of a real number.

It is possible to estimate the proportion of regions correctly classified as having μi among the top 10% of ranks by simulation. The MLE procedure (analogous to SMR) ranks the Xi. The empirical Bayes procedure ranks μ^iEB=E(μi|Xi)=ρiμ0+(1ρi)Xi, where ρi=σi2/(σi2+τ2). The PPHR procedure is based on 1000 simulations. In each simulation, a sample Si is drawn from the posterior distribution of μi, namely ϕ(μi; ρiμ0+(1−ρi)Xi, κi), where κi = τ2 ρi is the posterior variance. The {Si} are ranked and the regions in the top 10% of ranks recorded. This process is repeated 1000 times, and the 0.1I regions most often in the top 10% of ranks are designated as “high” regions.

The left three graphs in Figure 2 show the simulation results of the percentages of correctly identified high regions for three procedures. Without loss of generality, we have set μ0 = 0 and τ2 = 1. We can see that the percentage correctly classified increases with k=τ2/GM(σi2), which is a measure of signal to noise. GM(σi2) is the geometric mean of {σi2}, which equals exp(η). For [var phi]2 = 0.1, there is little variation among the {σi2}, and the PPHR, EB and MLE rankings perform similarly. With [var phi]2 = 1, PPHR and EB outperform MLE slightly. With [var phi]2 = 10, PPHR and EB are markedly superior to MLE, and PPHR outperforms EM slightly. Thus, when the variability in σi2 is small, all ranking methods perform comparably, whereas when the variability in σi2 is large, the PPHR method and EB are superior, and the more complex PPHR method is recommended when the variability of σi2 is very large (cf. Lin et al., 2006).

Figure 2
Gaussian-Gaussian model: left hand panels plot percent of regions correctly classified as having among the top 10% of means versus the signal to noise ratio, k; right hand panels plot the percent correctly classified estimated from simulations with 1000 ...

The following approximations are close to the results from simulations in Figure 2. For X ~ N(μ, σ2), with μ ~ N(μ0, τ2) and log(σ2) ~ N(η, [var phi]2), and assuming μ0, τ2, η and [var phi]2 are known, the joint density f (X, σ2, μ) is given by equation (6). Since the three ranking methods we consider are all monotone functions of X, we use R(X) to denote a generic statistic with distribution

GR(t)=P{R(X)t}=P{XR1(t)}=Eμ,σ2[P{XR1(t)|μ,σ2}]=σ2=0μ=x=R1(t)f(x,σ2,μ)dxdμdσ2,
(8)

and with pth percentile ζpR=GR1(p).

Suppose μ* is drawn from a N(μ0, τ2) distribution but truncated to be at or above the pth percentile δ = μ0 + zpτ, where zp = Φ−1(p) = 1.645, and Φ is the standard normal distribution. We approximate the probability that a region in the top 10% of risk will be correctly classified as high risk as the probability that an R(X*) corresponding to μ* exceeds ζpR, where p = 0.9. This probability is

P{R(X*)>ζpR}=1P{R(X*)ζpR}=1P{X*R1(ζpR)}=1Eμ*,σ2[P{X*R1(ζpR)|μ*,σ2}]=111pσ2=0μ*=δx*=R1(ξpR)f(x*,σ2,μ*)dx*dμ*dσ2.
(9)

f (X*, σ2, μ*) is given by equation (6) except that μ* is from the truncated normal distribution. The specific functions R(X) are defined as follows. The MLE procedure ranks RMLE(X) = X, the EB procedure ranks REB(X) = ρμ0 +(1−ρ)X, where ρ = σ2/(σ2 + τ2), and the PPHR procedure ranks RPPHR(X)=1Φ[{δμ0ρ(1ρ)X}/κ], with κ = ρτ2. As seen in the right three graphs in Figure 2, the approximation (8) gives results very close to the simulated values, as all pairs fall on the equiangular line by visual inspection.

4.2 Poisson–Gamma Model

We simplify the model in equations (1) and (2) to eliminate the covariate effect by setting r = 1 and by assuming ξj and hence Eij and Ei=j=1JEij are known. The underlying parameter of interest θi ~ Γ(α,1/α) as before, and conditional on θi, Oi is Poisson with mean θiEi. We further assume log(Ei) are normally distributed N(η, [var phi]2), so that Ei are lognormally distributed and independent of θi. The joint density is thus

g(Oi,Ei,θi)=Γ(θi;α,1/α)LN(Ei;η,[var phi]2)Poi(Oi;Eiθi).
(10)

We assume α, η, [var phi]2, and Ei are known. Selection of high risk regions is based on the SMRi = Oi/Ei, on the empirical Bayes estimate E(θi|Oi) = (α+Oi)/(α+Ei), and on the PPHR procedure. The posterior distribution of θi given Oi is Γ(α+Oi,1/(α+Ei)), and Oi are conditionally independent given {Ei} and α. Hence we can draw samples of Si from these posterior distributions and rank them. From 1000 sets of samples of this type, we determine which regions are most often in the top 10% of ranks and call them “high risk” regions.

In the simulation, we fixed α = 42.43, which is the MLE estimate from the heart disease mortality data, and we used [var phi]02=1.55, which was obtained by matching the first two moments of Ei from the heart disease mortality data to those of a lognormal distribution. The unconditional variance of SMRi = Oi/Ei is E(θiEi)/Ei2 = 1/Ei which has geometric mean exp(−η). The variance of θi is 1/α, yielding a signal to noise ratio k = (1/α)/exp(−η) = exp(η)/α.

In Figure 3, the left three graphs illustrate the simulation results corresponding to three degrees of variability in 1/Ei, namely [var phi]2=1.55/5=0.31, 1.55, and 1.55×=7.75. The percentage of high risk regions correctly classified increases with k for all procedures across different values of [var phi]2. When there is little variability in 1/Ei (i.e. [var phi]2=1.55/5), PPHR, EB and SMR perform comparably, but for [var phi]2=1.55, PPHR and EB outperform SMR slightly, especially for small k. For [var phi]2=7.75, the PPHR and EB are much better than SMR, even for large k. In Figure 3, we also indicated three points that correspond to 1/100, 1/10, and 1 times the death counts expected for heart disease. These points correspond respectively to an uncommon disease (e.g. esophagus cancer), a moderately common disease (e.g. breast cancer) and a common disease (e.g. heart disease). These results indicate that when there is variability in 1/Ei and the Poisson variation in the SMR estimate is large, PPHR or EB rankings outperform the SMR ranking, and , if the variability of 1/Ei is large, PPHR outperforms EB.

Figure 3
Poisson-Gamma model: left hand panels plot percent of regions correctly classified as having among the top 10% of relative risks versus the signal to noise ratio, k; right hand panels plot the percent correctly classified estimated from simulations with ...

For O ~ Poi(Eθ), with θ ~ Γ(α,1/α) and log(E) ~ N(η, [var phi]2), and assuming α, η and [var phi]2 are known, we approximate the performance of the three ranking methods, which are all monotone functions of O, denoted R(O) with

GR(t)=P(R(O)t)=P(OR1(t))=Eθ,E{P(OR1(t)|θ,E)}=E=0θ=0O=0{R1(t)}g(O,θ,E)dθdE,
(11)

and with pth percentile ζpR=GR1(p).

Suppose θ* is drawn from a Γ(α,1/α) distribution but truncated to be at or above the pth percentile δ = Γ−1(p; α,1/α). We approximate the probability that a region in the top 10% of risk will be correctly classified as high risk as the probability that an R(O*) corresponding to θ* exceeds ζpR, where p = 0.9. This probability is

P{R(O*)>ζpR}=1P{R(O*)ζpR}=1P{O*R1(ζpR)}=1Eθ*,E{P{O*R1(ζpR)|θ*,E}}=111pE=0θ*=δO*=0{R1(ζpR)}g(O*,θ*,E)dO*dθ*dE.
(12)

Specifically, the SMR procedure ranks RSMR (O) = O/E, the EB procedure ranks REB(O) = (O+α)/(E+α), and the PPHR procedure ranks RPPHR(O)=10δΓ(s;O+α,E+α)ds. As seen in the right three graphs in Figure 3, the approximation (12) gives results visually indistinguishable from the simulated data.

5. Discussion

We investigated procedures for identifying the regions with the upper 100γ% and lower 100γ% of standardized mortality rates. In an example of heart disease mortality (Section 3), we found that adjusting for smoking prevalence produced notable changes in the regions identified as high and low risk. Once smoking prevalence was accounted for, however, ranking based on the SMR, on empirical Bayes estimates (EB), or on posterior percentile ranking (PPR) yielded very similar classifications (Figure 1b, 1c and 1d). The EB and PPR methods outperformed SMR, however, when the expected counts were reduced by a factor of 10, (Table 1). In scenarios representing less common diseases, PPR and EB outperformed SMR when expected counts varied over regions (Section 4), and when there was wide variation in expected counts, the PPR method was better than EB. We developed very accurate approximations for the probability of correct classification in the upper 100γ% of risk for Gaussian-Gaussian and Poisson-Gamma models (Section 4).

The Poisson-Gamma empirical Bayesian framework was used by Clayton and Kaldor (1987). The choice of Gamma distribution for the random effect θi is mathematically convenient, because it is the conjugate prior for the Poisson distribution, and the posterior of θi is also gamma distributed. Clayton and Kaldor (1987) and Christiansen and Morris (1997) modeled the log-mean of θi as a linear combination of covariates. We suggested incorporating the effects of a covariate like smoking prevalence as a mixture that represented smokers and non-smokers. The Poisson-Gamma model helps to smooth θi by borrowing strength from the other regions. One could achieve further smoothing by allowing for spatial correlations. However, the Gamma prior cannot be easily generalized to for this purpose. In our analysis of mortality from heart disease, we aggregated data from nearby counties, reducing the number of units from 3,141 counties to 805 health service areas, and achieving some local smoothing.

Lin et al. (2006) presented data from a Gaussian-Gaussian model (their Table 2) that indicated better performance for picking high and low regions with B-PPR than by classification based on their posterior expected ranks. Conlon and Louis (1999) and Shen and Louis (1998, 2000) discussed procedures that simultaneously gave good estimates of region-specific rates, good estimates of the ranks of the rates, and good estimates of the cumulative distribution of the rates (“triple goal estimates”). Stern and Cressie (1999) described a loss function to identify regions with extreme risks.

Although the original B-PPR method was based on a fully Bayesian analysis, simulated comparisons using heart disease data in Section 3 suggests that our approximate PPR performs as well as or better than B-PPR with reduced computation time and easier implementation. Perhaps the fully Bayesian approach would perform better with different choices of parameters for the vague prior distributions and Gibbs sampling regime.

We assumed that the prevalence of smoking was known without error. Further work is needed to evaluate the impact of error and the use of hierarchical models to improve prevalence estimates. More experience is needed to assess the usefulness of the covariate adjustment and ranking procedures in this paper. The present work illustrates, however, that it is computationally feasible to adjust for covariates and use a preferred ranking procedure, PPR, which may be particularly advantageous for less common diseases with wide variation in expected counts.

ACKNOWLEDGEMENTS

The authors thank Professor Thomas A. Louis for pointing out some key references and providing reprints and the reviewers for helpful suggestions. This work was supported by the Intramural Research Program of the Division of Cancer Epidemiology and Genetics of the National Cancer Institute.

REFERENCES

  • Christiansen CL, Morris CN. Hierarchical Poisson Regression Modeling. Journal of the American Statistical Association. 1997;92:618–632.
  • Clayton DG, Kaldor J. Empirical Bayes estimates of age-standardized relative risks for use in disease mapping. Biometrics. 1987;43:671–682. [PubMed]
  • Conlon EM, Louis TA. Addressing multiple goals in evaluating region-specific risk using Bayesian methods. New York: Wiley; 1999.
  • Devesa SS, G DG, Blot WJ, Pennello G, Hoover RN, Fraumeni JF., Jr Atlas of Cancer Mortality in the United States, 1950–94. Washington, DC: 1999. US Govt Print Off [NIH Publ No. (NIH) 99-4564]
  • Laird NM, Louis TA. Empirical Bayes ranking methods. Journal of Educational Statistics. 1989;14:29–46.
  • Lin R, Louis TA, Paddock SM, Ridgeway G. Loss function based ranking in two-stage, hierarchical models. Bayesian Analysis. 2006;1:915–946. [PMC free article] [PubMed]
  • Lunn DJ, Thomas A, Best N, Spiegelhalter D. WinBUGS -- a Bayesian modelling framework: concepts, structure, and extensibility. Statistics and Computing. 2000;10:325–337.
  • Pennello GA, Devesa SS, Gail MH. Using a mixed effects model to estimate geographic variation in cancer rates. Biometrics. 1999;55:774–781. [PubMed]
  • Shen W, Louis TA. Triple-goal estimates for disease mapping. Statistics in Medicine. 2000;19:2295–2308. [PubMed]
  • Shen W, Louis TA. Triple-goal estimates in two-stage hierarchical models. Journal of the Royal Statistical Society Series B. 1998;60:445–471.
  • Stern HS, Cressie N. In: Inference for Extremes in Disease Mapping and Risk Assessment for Public Health. Lawson A, Biggeri A, Bohning D, Lesaffre E, Viel J-F, Bertollini R, editors. Chichester: John Wiley and Sons; 1999. pp. 63–84.