| 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. |
1. WHAT THIS SOFTWARE DOES.
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:
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.
3. EXAMPLES.
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
ciseason(events)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 |
OUTPUT:
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 |