Reducing Suicide: A National Imperative (2002)

One can naturally model the unobserved population heterogeneity or extra-population heterogeneity via mixture models. The simplest and most natural occurrence of the mixture model arises when one samples from a super population g which is a mixture of a finite number, say m, of populations (g₁,…, g_m), called the components of the population. Suppose a sample from a super population g is recorded as data (Y_i, J _i) for i = 1, …, n, where Y_i = y_i is a measurement on the i^th sampled unit and J_i = j_i indicates the index number of the component to which the unit belongs. If sampling is done from the j^th component, an appropriate probability model for the sampling distribution is given by

The function f_j(y_i; θ_j) represents a density function, called the component density for the i^th observation y_i and the parameter θ_j Θ_j^d_j called the component parameter, that describes the specific attributes of the j^th component population, g_j.

We treat the J_i as missing data and define the latent random variable Φ= (Ф₁,…, Ф_n) to be the values of the parameter θ {θ₁, …,θ_m} corresponding to the sampled components J₁, …, J_n; i.e., if the i^th observation came from the j^th component, then define Ф_i = θ_j. Then the Ф_i are a random sample from the discrete probability measure Q that assigns a positive probability π_j to the support point θ_j for j = 1, …, m. That is, the latent

random variable Φ has a distribution Q, which is called the mixing distribution, defined by P(Φ =θ_j) = Q({θ_j }) = π_j. Thus the mixing distribution

puts a mass (or weight) of π_j on each parameter θ_j. We thus arrive at the mixture distribution , or equivalently

One can find the posterior probability that the i^th observation comes from the j^th component (i.e, y_iG_j) using

In turn, one can assign each y_i to the population to which it has the highest estimated posterior probability of belonging; i.e., to g_j if τ_ij > τ_it for t = 1, …, m and t ≠ j.

One has to estimate the number of subpopulations or number of mixture components m as well as the parameters in the mixture distribution g…). This can be achieved via parametric or non-parametric maximum likelihood estimation (Laird, 1978; Böhning et al. 1992, Pilla and Lindsay, 2001). Laird (1978) considered nonparametric likelihood estimation of a mixing distribution and Lindsay (1983) linked its general theory, developed from a study of the geometry of the mixture likelihood, to the problem of estimating a discrete mixing distribution. One may also use an alternative approach to maximum likelihood for which menu driven software is provided by Böhning et al. (1998).

For a given m, one can estimate the parameters in the mixture distribution via the EM algorithm (Dempster et al. 1977) and several improvements of it which are discussed in McLachlan and Krishnan (1997). Let = ( , ) be the vector of estimated parameters. Once has been obtained, estimates of the posterior probability of population membership _ij can be formed for each y_i to give a probabilistic clustering.

When m, the number of mixtures, is large, the EM algorithm can be excruciatingly slow. In that case, one can use alternative EM methods proposed by Pilla (1997) and Pilla and Lindsay (2001) in the nonparamet-

ric mixture setting. Their approach assumes that the θ’s are known and the goal is to estimate the π’s only. However, by setting up the problem in a nonparametric framework, one may find the maximum likelihood estimator (MLE) quite easily in a high-dimensional problem via Pilla and Lindsay’s method. These methods are especially useful because θ₁ < θ₂ < …, i.e., the support points are ordered and hence the corresponding nearby densities can be paired to take advantage of the correlated densities.

Suppose data are available on the spatial occurrence of suicide in N counties in a state or for the entire country over a five year period. Assume that the counties are a mixture of two groups g₁ and g₂, corresponding to normal and high risk with respect to suicide, in proportions of π₁ and π₂ respectively. The number of events y_i in the i^th county is assumed to have a Poisson distribution with a rate of μ_j in g_j for j = 1, 2; i= 1, …, n. Let = (μ₁, μ₂, π₁, π₂)^T and θ_i = (μ₁N_i, μ₂N_i)^T, where N_i= M_it_i is the mass of population at risk in the i^th county assuming that there are M_i individuals at risk over time t_i. In this scenario, the log likelihood becomes

where f_j(y_i; θ_i) = exp(–μ_jN_i) (μ_jN_i)^y_i/y_i!. One can use the EM algorithm to find the MLE of the parameters. The complete data log likelihood function becomes

where z_ij = 1 if y_i g_j and 0 otherwise. The MLE of the parameters are found iteratively via

and

where which is the estimated posterior probability that y_ig_j.

Testing for the number of components in a mixture is a very hard problem. Consider the problem of testing the hypothesis H₀: number of components = m against H₀: number of components = m + 1. The likelihood ratio test (LRT) statistic defined as 2 log L_n = 2 [l(_m) – l(_m+1)], where l(_m) is the log likelihood log evaluated at _m, the nonparametric MLE under H₀ and similarly l(_m+1) is the nonparametric MLE under H₁. The LRT statistic generally has an asymptotic χ² distribution with d degrees of freedom, where d equals the difference between the number of parameters under the alternative and null hypothesis. However, this theory is known to fail for the mixture problem (Titterington et al., 1985:154). The reason for this is that the null hypothesis does not lie in the interior of the parameter space. The problem is not yet solved except for some special cases. However, the simulations or bootstrap can shed some light on the distributional properties of the LRT. See Lindsay (1995) and Böhning (1999) (including the references therein) for recent developments in theory and methods to test for the number of components in a mixture. There are other techniques for identifying the number of components of a mixing distribution including a variety of graphical techniques (Lindsay and Roeder, 1992).

The continuous EM algorithm that finds the MLE of parameter estimators, for a given m, requires the specification of initial values for the number and location of the support points. Pilla and Lindsay (2001) show that if the number of support points is fixed in advance, then one could have multiple modes whereas the nonparametric ML approach corresponds to a unimodal likelihood. Böhning (1999) illustrates through an example where the choice of initial values is critical; for example, in testing for the number of components in the mixture. Through simulations, he demonstrates the effectiveness of several approaches that combat this problem. His recommendation is to choose the support points as values of certain order statistics of the data (see Böhning, 1999:68–70). In his simulation studies, the gold standard (the grid approach given below) seems to outperform the other methods. The optimization problem solved by the grid-based EM (where the support points are assumed known) proposed by Pilla and Lindsay (2001) has a unique optimum which in turn generates an estimated number of components. However, the continuous EM can converge to the suboptimal mode.

In light of the above characteristics, Pilla and Lindsay (2001) suggest the following: to obtain a continuous solution with a flexible number of support points, first start with one of the fast EM-based algorithms proposed by Pilla and Lindsay on a grid and next use the resulting solution as starting values for the continuous EM. This will protect against using an incorrect number of points or finding suboptimal solutions.

Finding a natural stopping rule for iterative algorithms in the mixture problem when the final log likelihood is unknown is an important problem. In a potentially slow algorithm like EM, a far weaker convergence criterion based on the changes in the log likelihood value or on changes in the parameter estimates between iterations can be very misleading (Titterington et al. 1985:90). See Lindsay (1995:62–63) for further details. The following gradient-based stopping criterion proposed by Lindsay (1995:131–132) can be used to detect the convergence of the EM algorithm: we stop when sup D_Q(θ) ≤ ε, with ε = 0.005, where the gradient function D_Q(θ) is defined as

then we automatically satisfy the ‘ideal stopping criterion’: |log L_obs(; y) – log L_obs(^p; y)| ≤ ε, where ^p is the estimate at the p^th iteration. The gradient function therefore creates a natural stopping rule for iterative algorithms such as the EM when the final log likelihood is unknown, although it is more stringent than the ideal stopping rule.

To illustrate use of the Poisson mixture model we analyzed countylevel suicide data for the US for the time period of 1996–1998. The three year period was used to obtain more stable estimates of the parameters. Application of the Poisson mixture model identified a mixture of four component distributions with proportions of .1125, .4523, .3119, and .1234, and means of 7.5, 10.9, 14.9, and 20.6 suicides per 100,000. Results of classifying counties into component distributions are displayed in Figure 1. Inspection of the map in Figure 1 reveals that those counties classified in the highest component distribution (20.6 suicides per 100,000), are predominantly in the western states with lowest population density, although there are exceptions throughout the country (see Figure A-1). Counties classified in the lowest component distribution (7.5 suicides per 100,000)

FIGURE A-1 Annual Suicide Rates per 100,000 (1996–1998). The Poisson Mixture Model applied to county-level data.

appear to be clustered in the central United States. Finally, counties classified into the two intermediate groups are largely in the eastern portion of the United States. Unlike the earlier attempt at fitting a mixture distribution to temporal suicide data collected in Cook county (Gibbons et al., 1990) that was designed to identify “suicide epidemics” and failed to do so, the results of this analysis has focused on the spatial distribution of suicide and has identified qualitatively distinct geographic groupings of suicide rates across the United States. Note that this mixture distribution does not adjust for demographic disparity across the counties (e.g., age, sex, race, population density, etc.) that may have produced the mixture distribution in the first place.

In the following sections, we present statistical details of both fixed-effects and mixed-effects Poisson regression models.

In Poisson regression modeling, the data are modeled as Poisson counts whose means are expressed as a function of covariates. Let y_i be a response variable and x_i = (1, x_i1,…, x_ip)^T be a [(p+1) × 1] vector of covariates for the i^th individual. Then the Poisson regression model y_i is

(2)

where λ _i = exp(β^Tx_i) and β = (β₀, β₁, …, β_p)^T is a (p+1)-dimensional vector of unknown parameters corresponding to x_i. An important property of the Poisson distribution is the equality of the mean and variance:

From equation (2) one can write the likelihood and the log likelihood functions for the n independent observations as

and

respectively. The first and second partial derivatives with respect to the unknown parameters β are given by

and

respectively. Since the above equations are not linear, iterative procedures such as the Newton-Raphson algorithm or Fisher’s scoring algorithm are used to obtain the MLE of β, denoted by . A consistent estimator of the variance-covariance matrix of , V() is the inverse of the information matrix, I() — which is the negative of inverse of the matrix of second partial derivatives evaluated at the MLE. Inferences on the regression parameter β can be made using V().

When we have a mixture of fixed and random effects, the more general mixed-effects Poisson regression model is used. Suppose that there are n = Σ_in_i nonnegative observations y_ij for i = 1, …, n_c clusters and j = 1, …, n_i observations and a (p + 1) dimensional unknown parameter vector, β, associated with a covariate vector x_ij = (1, x_ij1, …, x_ijp)^T. Further, let the random effect υ_i be normally-distributed with mean 0 and variance σ² and independent of the covariate vector x_ij. Then the conditional density function, f(y_i; θ_i), of the n_i suicide rates in cluster i is written as:

(4)

with

where θ_i = υ_i/σ such that θ_i ~N (0,1). Thus the log-likelihood function corresponding to equation 4 becomes

The first and second derivatives of log L with respect to β and σ are

and

The marginal density of y_i can be written as

where g(θ) represents the multivariate standard normal density. One can approximate the above integral via the Gaussian quadrature as

(5)

where λ_ijq = exp(β^Tx_ij + σB_q), B_q is the quadrature node and A(B_q) is the corresponding quadrature weight. From equation (5), one can obtain the first derivatives of the approximate log likelihood as

and

Solution of the above likelihood equations can be obtained iteratively using the Newton-Raphson algorithm. On the i^th iteration, estimates for a parameter vector Θ are given by

(6)

where is the matrix of second derivatives of the log-likelihood, evaluated at Θ_i. Alternatively, one can use the Fisher scoring method. This method replaces the matrix of second derivatives in (6) by their expected values:

which is equal to the negative of the information matrix. The scoring solution is often used in cases where the matrix of second derivatives is difficult to obtain. Note that when the expected value and the actual value of the Hessian matrix coincide, the Fisher scoring method and the Newton-Raphson method reduce to the same algorithm.

The convergence of the iterative procedure depends on the initial values of the parameters. Good starting values reduce the number of iterations needed to reach the final estimates. It is sometimes possible that a poor choice of starting values may reach some local maximum instead of the global maximum. For this reason, it is often desirable to re-estimate the parameters under a variety of starting values, or to use the EM algorithm to obtain starting values that are reasonably close to the final solution (Hedeker and Gibbons, 1994).

Often, it is of interest to estimate values of the random effects θ_i within a sample. For this purpose the expected “a posteriori” (EAP) or empirical Bayes estimator _i has been suggested by Bock and Aitkin (1981). The estimator _i given y_i for cluster i can be obtained as:

Similarly, the variance of the posterior distribution of _i, which may be used to make inferences regarding the EAP estimator (i.e., the posterior variance is an estimate of precision of the estimate of _i), is given by

An alternative approach to the analysis of clustered suicide rate data is based on Generalized estimating equations (GEEs) models, which were introduced by Liang and Zeger (1986) and Zeger and Liang (1986).

Let y_ij be the j^th response for the i^th unit (e.g., county ) for i = 1, …, n_c and j =1, …, n_i observations and x_ij = (1, x_ij1, …, x_ijp ) be a vector of covariates for the i^th unit. Let y_i be the vector of the responses for the i^th unit with corresponding mean vector μ_i and let V_i be an estimate of the covariance matrix of y_i. Then the GEE marginal model for estimating β is given by

(7)

and

where g is the link function. Common choice for the link function might be the identity link for continuous data, log link for count data, and logit link for binary data. For example, the link functions for the Poisson and logistic regression models are g(a) = log(a) and g(a) = log(a/(1 – a)), respectively.

In addition to the marginal model, the covariance structure of the correlated observations for a given unit of y_i is modeled as

(9)

where A_i is a diagonal matrix of variance functions and R(a) is the working correlation matrix of y_i specified by the vector of parameters a. Various types of working correlation structures such as exchangeable or autoregressive can be used in the model.

The maximum likelihood estimator can be obtained by solving the above estimating equations iteratively:

(10)

Returning to the national suicide data from the previous section, we now illustrate how Poisson regression models can be used to estimate the effects of age, race, and sex on clustered (i.e., within counties) suicide rate data. For the purpose of illustration, we considered the effects of age divided into five categories (5-14, 15-24, 25-44, 45-64, and 65 and older), sex, and race (African American versus Other) in the prediction of suicide rates across the U.S. for the period of 1996-1998. These categories were used so that there would be sufficient sample sizes available to compare observed and expected annual suicide rates for both GEE and mixed-effects (maximum marginal likelihood - MMLE) Poisson regression models. To this end, we fit a Poisson regression model with all main effects and two-way interactions using both GEE and a full likelihood mixed-effects model. Given the large sample sizes almost all terms in the model were statistically significant although of widely varying effect sizes. Table 2 displays a comparison of parameter estimates and standard errors for the GEE and mixed-effects models. In general, the GEE and MMLE parameter estimates were remarkably similar. The only nonsignificant terms were two terms in the race by age interaction. The comparison of rates for ages 5-14 versus 15-24 and 5-14 versus 25-44, did not depend on race.

Based on the parameter estimates in Table A-1, we estimated marginal suicide rates for both GEE and mixed-effects models. Table A-2

displays observed and expected annual suicide rates for both methods of estimation, broken down by age, sex, and race. Inspection of Table A-2 reveals several interesting results. In general, suicide increases with age, is higher in males, and is lower in African Americans. Black females have the lowest suicide rates across the age range. In white males, the suicide rate is increasing with age whereas in all other groups, the suicide rate either is constant or decreases after age 65. Comparison of the expected frequencies for the GEE and mixed-effects models reveal that they are quite similar and the GEE does a slightly better job of recovering the observed marginal rates.

A special feature of the mixed-effects model is the ability of estimating county-specific rates using empirical Bayes estimates of the random effects as described in the previous section. Using these estimates, we can accomplish two goals. First, we can estimate county-specific expected suicide rates, which directly incorporate the effects race, sex, and age of that county. Table A-3 provides a comparison of observed and expected numbers of suicides (1996-1998) for 100 randomly selected counties. In-

spection of Table A-3 reveals remarkably close agreement between observed and expected numbers of suicides.

Second, we can use the Bayes estimates directly to obtain county-level suicide rates adjusted for the effects of race, sex, and age. In the case of a Poisson model, the Bayes estimate for a given county is a multiple of the national suicide rate adjusted for the case mix in that county (i.e., race, sex and age). For example, a Bayes estimate of 1.0 represents an adjusted rate that is equal to the national rate. By contrast, a Bayes estimate of 2.0 represents a doubling of the national rate, and a Bayes estimate of 0.5 represents one-half of the national rate. Figure A-2 displays the Bayes estimates by county across the U.S. Inspection of Figure A-2 reveals that even after accounting for these important demographic variables, considerable spatial variability remains. Again, the highest adjusted rates are typically found in the less densely populated areas of the western U.S.

The map in Figure A-2 also provides a useful tool for qualitative research into the etiology of suicide. A natural approach is to examine the

FIGURE A-2 Bayes Estimates of County-Level Deviations from the National Suicide Rates per 100,000 (1996–1998), Adjusted for Age, Sex, and Race.

spatial distribution of Bayes estimates in Figure A-2 for outliers. For example, in the western U.S. and Alaska where suicide rates are typically high there are a few counties that have Bayes estimates that are consistent with the national average. Similarly, in the central U.S. where there is a high concentration of counties with the lowest suicide rates, there are a few counties that exhibit the highest suicide rates. As an example, Table A-4 displays the Bayes estimates, observed and expected numbers of suicides and suicide rates for all counties in Alaska and New Mexico and the 8 counties with estimated adjusted suicide rates less than or equal to half of the national average (national average is 12.33 suicides per 100,000 during the period of 1996–1998). Table A-4 reveals that although there are generally elevated suicide rates in Alaska and New Mexico relative to the national average, there are several counties with low rates, similar to the national average. What are the protective factors that have produced these spatial anomalies? Are these spatial anomalies simply due to reporting bias or some other unmeasured characteristic that has produced the outlying Bayes estimate. Examining these spatial anomalies in greater detail is certainly a fruitful area for further research.

In motivating probit and logistic regression models, it is often assumed that there is an unobservable latent variable (y) which is related to the actual response through the “threshold concept.” For a binary response, one threshold value is assumed, while for an ordinal response, a series of threshold values γ₁, γ ₂, …, γ _K–1, where K equals the number of ordered categories, γ₀ = –∞, and γ _K = ∞. Here, a response occurs in category k (Y = k) if the latent response process y exceeds the threshold value γ _k–1, but does not exceed the threshold value γ _k.

The traditional ordinal regression model is parameterized such that regression coefficients α do not depend on k, i.e., the model asumes that the relationship between the explanatory variables and the cumulative logits does not depend on k. McCullagh (1980) calls this assumption of identical odds ratios across the K – 1 cut-offs, the proportional odds assump-

tion. This assumption implies that the effect of a regressor variable is the same across the cumulative logits of the model, or proportional across the cumulative odds. In the present context, this implies that the effect of an intervention is the same on the shift from no ideation to suicidal ideation, as it is from suicidal ideation to a suicidal attempt, an assumption that is questionable at best. As noted by Peterson and Harrell (1990), hovever, examples of non-proportional odds are not difficult to find. Recently Hedeker and Mermelstein (1998) have described an extension of the random effects proportional odds model to allow for non-proportional odds for a set of explanatory variables. For example, it sems reasonable to assume that the effects of the model covariates on suicidal ideation and suicidal attempts may not be the same. The resulting mixed-effects partial proportional odds model follows Peterson and Harrell’s (1990) extension of the fixed-effects proportional odds model.

Assuming that there are i = 1, …, N level-2 units and j = 1, …, n_i level-1 units nested within the i^th level-2 unit. The cumulative probabilities for the k ordered categories (k = 1, … K) are defined for the ordinal outcomes Y as:

where the mixed-effects logistic regression model for these cumulative probabilities is given as

with K – 1 strictly increasing model intercepts γ_k (i.e., γ₁ > … > γ_K–1). As before, w_ij is the p × 1 covariate vector and x_ij is the design vector for the r random effects, both vectors being for the j^th level-1 unit nested within level-2 unit i. Also, α is the p × 1 vector of unknown fixed regression parameters, and β_i is the r × 1 vector of unknown random effects for the level-2 unit i. Since the regression coefficients ,α do not depend on k, the model assumes that the relationship between the explanatory variables and the cumulative logits does not depend on k.

As in most mixed-effects models, it is convenient to orthogonally transform the response model. Letting β = Λt, where ΛΛ^T=Σ_β is the Cholesky decomposition of Σ_β, the reparameterized model is then written as log

where t_i are distributed according to a multivariate standard normal distribution.

To allow for a partial proportional odds model the intercepts γ_k are denoted instead as γ_0k, and the following terms are added to the model:

or absorbing γ_0k and into γ_k′

where, u_ij is a [(h+1) × 1] vector containing (in addition to a 1 for γ_0k) the values of observation ij on the subset of h covariates for which proportional odds are not assumed (i.e., u_ij is a subset of w_ij). In this model, γ_k is a (h+1) × 1 vector of regression coefficients associated with the h variables (plus the intercept) in u_ij. Notice that the effects of these h covariates vary across the (K – 1) cumulative logits. This extension of the model follows similar extensions of the ordinary fixed-effects ordinal logistic regression model discussed by Peterson and Harrell (1990) and Cox (1995). Terza (1985) discusses a similar extension for the ordinal probit regression model.

There is a caveat in this extension of the proportional odds model. For the explanatory variables without proportional odds, the effects on the cumulative log odds, namely result in (K – 1) non-parallel regression lines. These regression lines inevitably cross for some values of u*, leading to negative fitted values for the response probabilities. For u* variables contrasting two levels of an explanatory variable (e.g., gender coded as 0 or 1), this crossing of regression lines occurs outside the range of admissible values (i.e., < 0 or > 1). However, if the explanatory variable is continuous, this crossing can occur within the range of the data, and so, allowing for non-proportional odds is problematic. For continuous explanatory variables, other than requiring proportional odds, a solution to this dilemma is sometimes possible if the variable u has, say, m levels with a reasonable number of observations at each of these m levels. In this case (m – 1) dummy-coded variables can be created and substituted into the model in place of the continuous variable u.

With the above mixed-effects regression model, the probability of a response in category k for the i-th level-2 unit, conditional on γ_k, α, and t is given by:

where, under the logistic response function,

with Note that P_j0 and P_jK = 1.

Estimation of the p covariate coefficients α, the population variance-covariance parameters Λ (with r(r+1)/2 elements), and the (h + 1)(K –1) parameters in γ_k (k = 1, …, K – 1) is described in detail by Hedeker and Gibbons (1994).

To illustrate application of the mixed-effects ordinal logistic regression model, we reanalyzed longitudinal data from Rudd et al. (1996) on suicidal ideation and attempts in a sample of 300 suicidal young adults (personal communication, M.D. Rudd, Baylor University, October 2001). 180 subjects were assigned to an outpatient intervention group therapy condition and 120 subjects received treatment as usual. In this analysis, we used an ordinal outcome measure of 0=low suicidal ideation, 1=clinically significant suicidal ideation, and 3=suicide attempt. Suicidal ideation was defined as a score of 11 or more on the Modified Scale for Suicide Ideation (MSSI; Miller et al., 1986). Model specification included main effects of month (0, 1, and 6) and treatment (0=control, 1=intervention), and the treatment by month interaction. Although data at 12, 18 and 24 months were also available, the drop-out rates at these later months were too large for a meaningful analysis. In addition, to illustrate the flexibility of the model, depression as measured by the Beck Depression Inventory (BDI; Beck and Steer, 1987) and anxiety as measured by the Millon Clinical Multiaxial Inventory (MCMI-A; Millon, 1983) were treated as time-varying covariates in the model, to relate fluctuations in depressed mood and anxiety to shifts in suicidality.

Table A-5 displays the observed proportions in each of the three ordinal suicide categories as a function of month and treatment. Table A-5 reveals that, if anything, the observed rate of attempts is as high or higher in the intervention group than the control group. Note that this is also true at baseline, therefore, these difference may be an artifact of randomization.

Table A-6 displays the model parameter estimates, standard errors, and associated probability estimates. A model with both random intercept and slope (i.e., month) effects failed to converge because the variance component for the intercept was close to zero. In light of this, the model was specified with a single random effect (i.e., a random slope model). Table A-6 reveals that both of the time-varying covariates, depression and anxiety, were significantly associated with suicidality. By contrast, the treatment by time interaction was not significant, indicating that the intervention did not significantly affect the rates suicide ideation or at-

tempts over time. The random month effect was significant, indicating appreciable inter-individual variability in the rates of change over time.

It might be argued, that by including the effects of depression and anxiety in the model, the effect of treatment and the treatment by time interaction may be masked. For example, the intervention may work by modifying depression and anxiety which are in turn related to suicide ideation and attempts. In this case, including depression and anxiety in the model may mask the effect of the intervention. Reanalysis excluding depression and anxiety, however, also yielded nonsignificant treatment related effects, suggesting that this is not the case. Finally, a nonproportional odds model, in which the treatment by time interaction was allowed to be category-specific, also failed to identify any positive treatment related effects on suicide ideation or attempts.

A (1 –α)100% confidence interval for the true suicide rate λ can be approximated as

(11)

where = y/n, y is the total number of suicides over n time intervals (e.g., months) or geographic areas (e.g. counties), and z is an upper percentage point of the normal distribution. In general, the approximation will work extremely well for y ≥ 20. When y < 20 tabled values are available in Hahn and Meeker (1991).

As an example, consider the case in which after 61 months of observation for a particular county, 123 suicides are recorded; therefore, y = 123 and n = 61. The 95% confidence limit for the population suicide rate λ is therefore

(12)

In light of this, we can have 95% confidence that the true suicide rate is between 1.659 and 2.373 suicides per month.

Cox and Hinkley (1974) consider the case in which y has a Poisson distribution with mean µ. Having observed y their goal is to predict y*, which has a Poisson distribution with mean cµ, where c is a known constant. In the present context, y is the number of suicides observed in n previous time-periods, and y* is the number of suicides observed in a single future time-period, therefore, c = 1/n. Alternatively, for a given time period, n may represent a number of counties, and our inference is to the predicted suicide rate in a single future county for the same length of time.

Following Cox and Hinkley (1974), Gibbons (1987) derived a prediction limit using as an approximation the fact that the left side of:

(13)

is approximately a standard normal deviate, therefore, the 100 (1 –a)% prediction limit for y* is formed from percentage points of the normal distribution. Upon solving for y* the upper limit value is found as the positive root of the quadratic equation:

(14)

Returning to the previous example in which after 61 months of observation for a particular county, 123 suicides are recorded, we have y = 123 and c = 1/61. The 95% confidence Poisson prediction limit is given by

The prediction limit reveals that we can have 95% confidence that no more than four suicides will occur in the next month.

The uniformly most accurate upper tolerance limit for the Poisson distribution is given by Zacks (1970). In terms of a future time interval, we can construct an upper tolerance limit for the Poisson distribution that will cover P(100)% of the population of future measurements with (1 – a)100% confidence. The derivation begins by obtaining the cumulative probability that a or more suicides will be observed in the next time period:

(15)

which can be approximated by:

where χ²[f] designates a chi-square random variable with f degrees of freedom. This relationship between the Poisson and chi-square distribution was first described by Hartley and Pearson (1950). Given n independent and identically distributed Poisson random variables the sum

(17)

also has a Poisson distribution. Substituting T_n for μ, we can find the value for which the cumulative probability is 1 – α, that is,

(18)

The P(100)% upper tolerance limit is therefore Pr^–1[P; K_1–α(T_n)] = smallest j (≥ 0) such that

(19)

The required probability points of the chi-square distribution can be most easily obtained using the Peizer and Pratt approximation described by Maindonald (1984:294).

Returning to the previous example, where after 61 months of observation for a particular county, 123 suicides are recorded, we have, T_n = 123 and n = 61. Furthermore, we require the resulting tolerance limit to have 99% coverage and 95% confidence. The cumulative 95% probability point is

The 99% upper tolerance limit is obtained by finding the smallest nonnegative integer j such that:

Inspection of standard chi-square tables (extracted below for j = 5 to 7) reveals that the value of j that satisfies this equation is 7.

Therefore, the P^–1 [.99 ; K_.95(123)] upper tolerance limit is 7 suicides per month in the example posed.

Beck AT, Steer RA. 1987. Beck Depression Inventory: Manual. San Antonio, TX: Psychological Corporation.

Bock RD, Aitkin M. 1981. Marginal maximum likelihood estimation of item parameters: an application of the EM algorithm. Psychometrika, 46: 443–459.

Böhning D, Schlattmann P, Lindsay BG. 1992. Computer-assisted analysis of mixtures (C.A.MAN): Statistical algorithms. Biometrics, 48: 283–303.

Böhning D, Dietz E, Schlattmann P. 1998. Recent developments in computer-assisted analysis of mixtures. Biometrics, 54: 525–536.

Böhning D. 1999. Computer-Assisted Analysis of Mixtures and Applications. Meta-Analysis, Disease Mapping and others. Chapman and Hall CRC: Boca Raton.

Cox C. 1995. Location-scale cumulative odds models for ordinal data: A generalized nonlinear model approach. Statistics in Medicine, 14: 1191–1203.

Cox DR, Hinkley DV. 1974. Theoretical Statistics. Chapman and Hall: London.

Dempster AP, Laird NM, Rubin DB. 1977. Maximum likelihood from incomplete data via the EM algorithm (with discussion). Journal of the Royal Statistical Society, Series B, 39: 1–38.

Gibbons RD. 1987. Statistical models for the analysis of volatile organic compounds in waste disposal sites. Ground Water, 25: 572–580.

Gibbons RD, Clark DC, Fawcet J. 1990. A statistical method for evaluating suicide clusters and implementing cluster surveillance. American Journal of Epidemiology, 132 (Suppl 1): S183–S191.

Hahn GJ, Meeker WQ. 1991. Statistical Intervals: A Guide for Practitioners. Wiley: New York.

Hartley HO, Pearson ES. 1950. Tables of the chi-squared integral and of the cumulative Poisson distribution. Biometrika, 37: 313–325.

Hedeker D, Gibbons RD. 1994. A random effects ordinal regression model for multilevel analysis. Biometrics, 50: 933–944.

Hedeker D, Mermelstein RJ. 1998. A multilevel thresholds of change model for analysis of stages of change data. Multivariate Behavioral Research, 33: 427–455.

Laird N. 1978. Nonparametric maximum likelihood estimation of a mixing distribution. Journal of the American Statistical Association, 73: 805–811.

Liang KY, Zeger SL. 1986. Longitudinal data analysis using generalized linear models. Biometrika, 73: 13–22.

Lindsay BG. 1983. The geometry of mixture likelihoods, part I: A general approach. Annals of Statistics, 11: 783–792.

Lindsay BG. 1995. Mixture Models: Theory, Geometry and Applications. NSF-CBMS Regional Conference Series in Probability and Statistics, Vol. 5. Hayward, CA: Institute of Mathematical Statistics.

Lindsay BG, Roeder K. 1992. Residual diagnostics for mixture models. Journal of the American Statistical Association, 87: 785–794.

Maindonald JH. 1984. Statistical Computation. Wiley: New York.

McCullagh P. 1980. Regression models for ordinal data (with discussion). Journal of the Royal Statistical Society, Series B, 42: 109–142.

McLachlan GJ, Krishnan T. 1997. The EM Algorithm and Extensions. New York: Wiley.

Miller IW, Norman WH, Bishop SB, Dow MG. 1986. The Modified Scale for Suicidal Ideation: Reliablity and validity. Journal of Consulting and Clinical Psychology, 54(5):724– 725.

Millon T. 1983. Millon Clinical Multiaxial Inventory Manual. 2nd ed. Minneapolis, MN: National Computer Systems.

Peterson B, Harrell FE. 1990. Partial proportional odds models for ordinal response variables. Applied Statistics, 39: 205–217.

Pilla RS. 1997. Improving the Rate of Convergence of EM in High-Dimensional Finite Mixtures. Ph.D. dissertation, The Pennsylvania State University.

Pilla RS, Lindsay BG. 2001. Alternative EM methods for nonparametric finite mixture models. Biometrika, 88: 535–550.

Rudd MD, Rajab MH, Orman DT, Joiner T, Stulman DA, Dixon W. 1996. Effectiveness of an outpatient problem-solving intervention targeting suicidal young adults: Preliminary results. Journal of Consulting and Clinical Psychology, 64(1):179–180.

Terza JV. 1985. Ordinal probit: A generalization. Communications in Statistical Theory and Methods, 14: 1–11.

Titterington DM, Smith AFM, Makov UE. 1985. Statistical Analysis of Finite Mixture Distributions. New York: Wiley.

Zacks S. 1970. Uniformly most accurate upper tolerance limits for monotone likelihood ratio families of discrete distributions. Journal of the American Statistical Association, 65: 307–316.

Zeger SL, Liang KY. 1986. Longintudinal data analysis for discrete and continuous outcomes. Biometrics, 42: 121–130.