Software for confidence intervals for seasonal realitive risk
Frangakis, CE, and Varadhan,  R (2002).  On Confidence Intervals for Seasonal Risk of Suicides, with Null Values on the Boundary (Epidemiology,  forthcoming [ Full text - PDF]). The authors have no responsibility for any  consequences of use of this software.
SOFTWARE  FILES:   routinesu.s, routines.o   routinesw.s, routines.dll

With count data (n1,...,n{seasons} )  for ordered seasons 1,...,maxseason, we are often interested in estimating the Relative Risk  (RR)=max/min frequency of the counts. In this case, the null value RR=1 is a boundary value, and although methods for testing exist, standard methods for confidence intervals (CI) are generally not appropriate. This software uses the method based on a circular model described in the above paper and obtains CIs that are valid:

  1. for all values of the underlying seasonal risk if the model is correct;
  2. for the null boundary value of the seasonal risk, regardless of model assumptions and sample size.
The software estimates: the Relative Risk  (RR)=max/min frequency of the counts; the exact 95% intervals for RR; the season, say beta, of maximum peak in degrees  (/360) on the circle; and asymptotic confidence intervals for beta.

2. HOW TO USE THE SOFTWARE. The current version is for unix with S-plus,
or windows with R.

If you have unix with S-plus:                            If you have windows with R:

1) save the files routinesu.s,                                   1) save the files routinesw.s,
and routines.o in an empty                                     and routines.dll in an empty
directory in your computer                                    directory in your computer
(in either case, make sure the saving of the files does not add any other extensions such as .txt)

2) from S-plus type:                                                2)  from R type
source("routinesu.s")                                               source("routinesw.s")
dyn.load("routines.o")                                            dyn.load("routines.dll")
The last command is needed every time you enter S-plus or R.

3) In S-plus or R, put in memory the counts you want to analyze, and then type the main command ``ciseason", to get the results.

example (i).  Suppose we want to assess the seasonality for the suicide counts: 2  ,7 ,14 ,14 , 9 ,10 , 8 ,10 , 6  ,6,  5 , 9 obtained say in the months 1-12. By plotting first the counts against the months, we notice an apparent peak around months 3 and 4. The method of the paper can then  assess if the apparent seasonal pattern could be due to chance. Storing the counts and calling the function ciseason is done by

events <- c(2,7,14,14,9,10,8,10,6,6,5,9)
ciseason(events,maxiter=10,nsim=40000,.001,exact=T), or simply
The last command produces, after the computer stops the calculation, the table
RR lower(RR) upper(RR) P(RR) date lower(date) upper(date)
2.17 1.18 3.83 0.028 128 86.7 170


RR:  the estimate for the relative risk (maximum over minimum probability of having the event across seasons, which equals exp(2k)).

lower(RR), upper(RR):     the lower and upper endpoints of the 95% CI for RR using the method of Frangakis and Varadhan (2002).

P(RR): p-value (note: rarely, it may happen that P(RR) <.05 with lower(RR) =1.00 to 2 d.p. due to rounding).

date: the estimate for beta (date, or angle in degrees) of maximum probability of having the event. Although, for ease of computation, the estimate of beta is given in the scale 0 to 360, the actual parameter space for beta is multiple integers of 30 (0, 30, 60, etc.), so the estimate should be transformed to that space.

lower(date), upper(date): the lower and upper endpoints of the 95% CI for beta using the normal approximation. These have meaning only when there is evidence for seasonality.

(Note: because the result is based on a simulation, it is recommended that the user repeats the calculation a second
time: if the results between the two computations are not close enough to the desired level of the user, then the user
can set  the control parameter ``nsim'' to a higher number.)

example (ii).  When we apply the same method to the counts:  9  6 10 13 13  5  9  8  7  8  9  3  we get the following results.
RR lower(RR) upper(RR) P(RR) date lower(date) upper(date)
1.61 1 2.82 0.252 125 57.5 192

Return to Dr. Frangakis's  Home Page