Confidence Intervals for Seasonal Relative Risk with Null Boundary Values

Abstract

In evaluating the relative risk of seasonality, the null value is a boundary value. In this case, tests for the null hypothesis exist, but standard methods for confidence intervals are not appropriate. We provide a method for constructing confidence intervals under the circular normal model. The proposed confidence intervals are valid for all values of the underlying seasonal risk if the model is correct and for the null boundary value of the seasonal risk, regardless of model assumptions and sample size. We apply our method to seasonal suicide data from a recent report.

Petridou et al. ¹ 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, ² 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.

Methods

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 n_i denote the number of suicides for i = 1, . . . , 12 months. We assume that the counts (n₁, . . . , n₁₂) are from a multinomial distribution with probabilities of suicide (p₁, . . . , p₁₂), where p_i are proportional to the periodic factor e^{k cos (2π(i − β)/12)};e^2k 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. ³

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, p_i, 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 (k₀, β₀) are the true values of the parameters (k, β) that generate the observed suicide counts (n₁,…n₁₂). We can write the likelihood function for any parameter values (k, β) as

We wish to construct CIs for k₀, 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, ² 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 (n₁, . . . , n₁₂) 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 (ˆ, β₀) is an exact 95% CI for the true value k₀, whether or not k₀ lies on the boundary.

Generally, finding the CI for k₀ using S (ˆ, β₀) is difficult because the true peak month β₀ 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 (ˆ, β₀) =S (ˆ, ˆ;&bgr) so that, instead of the unknown month β₀, we can simply use the MLE, ˆ;&bgr, to calculate the set S (ˆ, ˆ;&bgr), which will then be the exact 95% CI for k₀.

Algorithm

We have developed an algorithm to compute S (ˆ, ˆ;&bgr), the exact 95% CI for k₀, 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. ⁴ For the larger root, k_u, of the equation f (k_u) = 0, the secant method starts with two points, k_u⁽¹⁾ and k_u⁽²⁾, and proceeds iteratively to obtain a sequence k_u⁽ⁱ⁾, i = 3, 4,. . ., that converges to k_u. Here we choose k_u⁽¹⁾ and k_u⁽²⁾ to be around the normal approximations ˆ + 1.96 se (ˆ) and ˆ + 1.96 se (ˆ) + Δ, respectively, where Δ is a small number. The smaller root, k_l, is obtained similarly where we choose the two starting points as k_l⁽¹⁾ = ˆ − 1.96 se (ˆ) and k_l⁽²⁾ =k_l⁽¹⁾ + Δ. 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.

Results

We applied our method to the suicide counts used in Petridou et al. ¹ 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. ¹)

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. ¹ 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, β₀ (in degrees).

Table: Table. Results for Relative Risk of Suicide (Maximum over Minimum Frequency) Within Countries, Using the New Method Applied to the Data of Petridou et al. ¹

Discussion

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/

Acknowledgments

The authors thank the Editor and three reviewers for constructive comments, and also Drs. D. Trichopoulos, E. Petridou, and K. Rothman for fruitful discussions.