Petridou et al. 1 recently evaluated the seasonality of suicides, showing evidence of association with months of higher sunshine. A reviewer of that paper had pointed out an interesting issue: how can we construct a confidence interval (CI) for the relative risk (RR) of suicide (maximum over minimum frequency) when the null value of RR (unity) is also a boundary of the possible values of RR? Although approaches exist for testing some boundary null hypotheses, 2 we are unaware of an application in this area that computed CIs when the null value is on the boundary. Therefore, in the original article, we first tested the null, finding strong evidence for seasonality for all countries, and then constructed the CIs using the maximum asymptotic margin of error across countries. This minimized potential false positive findings. Since then, we have developed and applied a more rigorous method to address the problem, which we briefly outline below. The results from the new method strengthen the substantive findings of the original paper.
Data and Model
The data are the numbers of suicides in each of the 12 months for twenty countries, for a period covering 4 to 24 of the most recent years. To demonstrate the new model, we focus on one such country, and let ni denote the number of suicides for i = 1, . . . , 12 months. We assume that the counts (n1, . . . , n12) are from a multinomial distribution with probabilities of suicide (p1, . . . , p12), where pi are proportional to the periodic factor ek cos (2π(i − β)/12);e2k is the ratio of maximum over minimum frequency of suicide and β is the month of maximum frequency of suicides, where k is nonnegative and β is an integer between 1 and 12. This model is a discrete version of the circular normal distribution for continuous data. 3
The model above has several notable features. It allows for the null hypothesis of no seasonality of suicides, when k = 0, whereby all the probabilities, pi, are equal to 1/12. Furthermore, the above model will be more powerful in detecting true seasonality compared to a less structured model, such as the saturated model that fits 11 different probabilities. The model also possesses a useful property that allows us to construct CIs for the risk parameter, k, when standard methods are not generally appropriate, as discussed next.
Confidence Intervals when the Null Risk Is on the Boundary
Suppose (k0, β0) are the true values of the parameters (k, β) that generate the observed suicide counts (n1,…n12). We can write the likelihood function for any parameter values (k, β) as
We wish to construct CIs for k0, based on the maximum likelihood estimates (MLE), ˆ and ˆ;&bgr, which are the values of k and β that maximize the likelihood function in Eq 1. However, because the null value, k = 0, is a boundary value, the usual theory for CIs of the form ˆ ± 1.96 se (ˆ), obtained using a normal approximation of ˆ, is not applicable here. Moreover, although approaches exist for testing some boundary null hypotheses, 2 we are unaware of a previous application in this area that has offered CIs when the null value is on the boundary. Here we propose an approach for valid CIs with our model based on first principles.
For this approach, it is important to consider, for every pair (k, β), the quantile function q (k, β) defined so that Pr (‖ˆ−k ‖≤q (k, β)‖k, β) is 95%. In words, the quantile q (k, β) is the number q such that, when counts of total size n are generated repeatedly from model 1 with true parameter values (k, β), the fraction of times that ‖ˆ −k ‖ is less that q is 95%. For a fixed pair (k, β), we can calculate q (k, β) as follows: (i) simulate a large number of sets of (n1, . . . , n12) from model 1, each set totalling the fixed size n; (ii) for each set, calculate its MLE ˆ and its distance from k, ‖ˆ −k ‖ and (iii) set q (k, β) equal to the smallest number below which fall 95% of the calculated distances ‖ˆ −k ‖ across the simulations.
The quantile function gives an exact reference for how far the estimated risk parameter is expected to be from the true value. In particular, if for any month β we define the set S (ˆ, β) of values of the risk parameter k to be MATH
then, by result of inversion of tests, we have that the range of values k defined by the set S (ˆ, β0) is an exact 95% CI for the true value k0, whether or not k0 lies on the boundary.
Generally, finding the CI for k0 using S (ˆ, β0) is difficult because the true peak month β0 is unknown. However, model 1 has the property that, for fixed risk parameter k, the quantile function q (k, β) is not a function of β. The intuition for this is that, given the risk parameter k, the likely values of its estimator (ˆ) are invariant of the month at which the peak occurs (a proof is given in the Appendix). It follows that S (ˆ, β0) =S (ˆ, ˆ;&bgr) so that, instead of the unknown month β0, we can simply use the MLE, ˆ;&bgr, to calculate the set S (ˆ, ˆ;&bgr), which will then be the exact 95% CI for k0.
We have developed an algorithm to compute S (ˆ, ˆ;&bgr), the exact 95% CI for k0, that involves three distinct computational tasks. The first is to compute the MLEs, ˆ and ˆ;&bgr. The second is to obtain the two roots of the equation, f (k) = 0, where f (k) = ‖k − ˆ‖ minus;q (k, ˆ;&bgr); the larger root is the upper endpoint of the 95% CI; the lower endpoint is the smaller root, provided that the root is positive, and is zero otherwise. In this task, the roots of the equation are obtained iteratively using the secant method. 4 For the larger root, ku, of the equation f (ku) = 0, the secant method starts with two points, ku(1) and ku(2), and proceeds iteratively to obtain a sequence ku(i), i = 3, 4,. . ., that converges to ku. Here we choose ku(1) and ku(2) to be around the normal approximations ˆ + 1.96 se (ˆ) and ˆ + 1.96 se (ˆ) + Δ, respectively, where Δ is a small number. The smaller root, kl, is obtained similarly where we choose the two starting points as kl(1) = ˆ − 1.96 se (ˆ) and kl(2) =kl(1) + Δ. A third task of our algorithm is required inside the secant method: because the quantile function, q, is not available in closed form, we compute q by simulation from the model at the values required by the iterations of the secant method. Although the parameter space for β is the integers, for ease of computation here we use the continuous scale.
In some simulations, we found that the endpoints of the exact 95% CI were relatively close to their starting normal approximations. However, as also discussed ear- lier, validity of the normal approximation cannot be guaranteed for arbitrary data with this model, and so the use of the exact 95% CIs here is generally preferable.
We applied our method to the suicide counts used in Petridou et al. 1 across the 12 months, for each country. The sample counts for all the months within a country ranged between 3,770 (Greece) and 222,227 (Japan), as summarized over recent years. (For further details, see Petridou et al. 1)
The CIs for the relative risk computed with the new method are presented in the table. As expected, the point estimates of RR remain essentially the same as in Petridou et al. 1 However, the exact CIs for RR computed with the new method are much shorter than the conservative ones originally reported and are theoretically valid for the boundary constraint, a fact that further strengthens the original findings. Because we find strong evidence for seasonality (k > 0) with this method, as well, we also report the normal approximation to the CIs for the date of peak incidence of suicide, β0 (in degrees).
We present a method that defines the relative risk in seasonality through a discrete version of the circular normal distribution, providing CIs that are valid no matter the underline value of that risk. Moreover, because the circular normal model is correct under the null hypothesis of all seasons having the same risk, our CIs are guaranteed to have the right coverage under the null hypothesis, regardless of model or sample size assumptions.
The method we discussed, illustrated here for 12 seasons (here, months) can be used for any arbitrary number of seasons when the researcher expects a single peak in the circular frequency of the studied event. Extensions of the method for CIs for multiple peaks, and for accounting for other covariates simultaneously, may be possible by introducing appropriate terms in the likelihood.
The computer program implementing the new method for obtaining CIs for a seasonal risk with general number of seasons is available from the Web site: http://biosun01. biostat.jhsph.edu/∼cfrangak/papers/suicide/
The authors thank the Editor and three reviewers for constructive comments, and also Drs. D. Trichopoulos, E. Petridou, and K. Rothman for fruitful discussions.
To show that the quantile function q (k, β) is free of β, it suffices to show that the distribution of the MLE (ˆ) is the same when the data are generated by (k0, β0) or by (k0, β0 + 1), ie, that for any fixed k*, MATH
The key to showing (2) is to note that, for any 12-count n = (n1,. . ., n12) and (k, β), the model 1 is circularly invariant, that is MATH
where ni(+1) (i = 1,...,12) are the counts of the translated configuration n(+1) defined as (n12, n1,..., n11). This means that, given a fixed true value of the risk parameter k, an observation of a specific ordered pattern for the counts n when the true peak date is β is equally likely to an observation of that pattern rotated by an amount of time when the true peak date is also rotated by the same amount of time. In particular, denote by Nk* to be all configurations n that give ˆ(n) = k*. Then, the RHS and LHS of (2), respectively, become pr (n ∈Nk* ‖k0, β0 + 1) and pr (n ∈Nk* ‖k0, β0). But, by invariance (3), the latter equals pr (n ∈Nk*(+1) ‖k0, β0 + 1), where Nk*(+1) has the n(+1) configurations of Nk*. Therefore, to show (2), it suffices to show that Nk* =Nk*(+1), or, that if n belongs to Nk* then it also belongs to Nk*(+1) and vice versa. We show the first.
Take n in Nk*, ie, ˆ(n) =k*, and say the MLE of the date, ˆ;&bgr(n) = β*. Now, letting n(−1) = (n2,..., n12, n1), we have that:MATH
(by invariance of model in ) MATH
(by invariance property of MLEs) MATH
(by invariance of model in )
which shows that k* is then also the MLE of k with data n(−1), so n(−1) belongs in Nk*, and so n(−1 + 1) =n must also belong in Nk*(+1), as well as in Nk*. With a similar argument we can show that if n is in Nk*(+1) then n must also be in Nk*, which proves (2).
1. Petridou E, Papadopoulos FC, Frangakis CE, Skalkidou A, Trichopoulos D. A role of sunshine in triggering of suicide. Epidemiology 2002; 13: 106–109.
2. Self SG, Liang KY. Asymptotic properties of maximum likelihood estimators and likelihood ratio tests under nonstandard conditions. J Am Statist Assoc 1987; 82: 605–610.
3. Gumbel EJ, Greenwood JA, Durand D. The circular normal distribution: theory and tables. J Am Statist Assoc 1953; 48: 131–152.
4. Burden RL, Faires JD. Numerical Analysis
. 7th ed. Pacific Grove, CA: Brooks/Cole Publishing, 2000.