# Calculating summary statistics on a distribution of "day of year" data [R]

Say I have a dataset that consists of 1,000 integers between 1 and 365 that represent days of the year a certain event happened. I am trying to figure out how I could calculate summary statistics of the distribution while respecting the fact that the data are cyclic (i.e. 1 and 365 are neighbors). I specifically want to calculate the median value and the values of at specific quantiles (.025 and .975).

Here is an example of a distribution I'd like to calculate these values for:

I can't quite figure out how to do this. I think I may need to do something like map the values to a circle, but this is not something I am familiar with (and R libraries like circular do not have much documentation that pertains to this sort of situation).

I am using R, so I would very much appreciate possible R code solutions!

EDIT
I am going to add a bit more information to make it clear why I want to do this.

For simplicity, say I am trying to model which day of the year kids are happiest. I have happiness data, and I generate a model that estimates happiest doy. Now I want to know the uncertainty of this estimate.

These 1000 values are draws from a bootstrap that estimates the happiest doy for each bootstrap. I want to summarize this info using the median value and the 95% CI for this estimate.

I don't want to simply show the distribution, because say I am repeating this for 100 different countries and want to just show the median+CI in a table.

• The fundamental question is: what do you hope a median or quantile would represent? How do you define them for circular data?
– whuber
Commented Aug 2, 2023 at 21:18
• I want to describe this distribution. In actuality, these values are akin to a bootstrap distribution of estimates of a day of the year at which a certain event happens. I want to be able to say what the average value across 1000 bootstraps is and what the CI is for that average value. Does that make sense? Commented Aug 3, 2023 at 2:42
• Your graph already is a great description. Generally, any statistical description focuses on a small number of salient or important properties of the dataset or its inferred data-generating process. You have now changed the properties you refer to (from percentiles to an average), so it's worth asking for confirmation: what kind of "average value" do you seek? One way to answer that would be to explain how you plan to interpret it (or how you intend your audience to interpret it).
– whuber
Commented Aug 3, 2023 at 12:16
• I have a number of these distributions, and I want to be able to show how they differ at a glance (say in a table). For this purpose, I was thinking of showing median value and CI for each one. Commented Aug 3, 2023 at 14:25
• Isn't a graphic ideal for "at a glance" assessment?
– whuber
Commented Aug 3, 2023 at 15:08

You specify more than summarizing a distribution: this is a distribution relative to the day the event happened. You wish, then, to summarize the differences. Although these are still circular data, they have a natural origin (0) and are signed.

This enables us to make sense of your original request. When you sort the signed differences (lying between $$-183$$ and $$+183$$), the median is the middle value. A nonparametric confidence interval can be found in the standard way using order statistics as described at https://stats.stackexchange.com/a/284970/919 and https://stats.stackexchange.com/a/160148/919. BTW, these are not the same as the $$1-\alpha/2$$ and $$\alpha/2$$ quantiles of the distribution.

There's one tiny technical issue: what do you do when your bootstrap is diametrically opposite the reference value? That is, it's exactly 183 days away (let's suppose this is even possible because perhaps you are considering leap years) and you need to represent it as a deviation of either $$+183$$ or $$-183.$$ I would suggest, even though this is ad hoc, splitting the difference and introducing both values into your dataset, each with weight $$1/2.$$

Why?

1. This is superior to just dropping such values, because it retains all the data, which will yield more accurate values for the confidence interval.

2. It implicitly supposes your error distribution is approximately symmetric. Thus, the likelihood that this is a positive error equals the likelihood it is a negative error.

3. It's not really any more difficult to compute weighted quantiles and weighted order statistics: you just put the weights of your data in the same order as the data and compute cumulative sums.

If your bootstrap is producing a highly asymmetrical error distribution, then you would need to do something more sophisticated. But in that case you likely need a more complicated summary of the results in the first place.

An alternative, which is fairly robust because these are circular data (thereby with bounded, compact support), is to use the square root of a circular variance. Convert each of your $$n$$ deviations (a signed number of days) into an angle $$\theta$$ (divide by 365, 366, or 365.25 as appropriate, then multiply by $$2\pi$$), convert that to Cartesian coordinates $$(x,y)=(\cos\theta, \sin\theta),$$ compute their means $$(\bar x,\bar y)$$ separately for each coordinate, and convert that to an angle and back to a number of days. This estimates the bias.

For the circular dispersion, set $$R^2 = (\bar x)^2 + (\bar y)^2$$ to be the squared length of that mean vector and report $$\sqrt{-\log(R^2)},$$ converted back to days.

To illustrate, and provide all details, here is R code to compute the circular mean and dispersion. The function expects your data vector in theta and the length of your year (in any units you desire) in modulus. It returns a vector of two values: the circular mean and circular dispersion.

summary.circular <- function(theta, modulus = 2 * pi) {
angles <- theta / modulus * 2 * pi # Convert to angles
xy <- cbind(cos(angles), sin(angles))         # Cartesian coordinates
xy.bar <- colMeans(xy)
R2 <- sum(xy.bar^2)
c(Mean = atan2(xy.bar[2], xy.bar[1]) / (2 * pi) * modulus,
Dispersion = sqrt(-log(R2)) / (2 * pi) * modulus)
}


As an example of how this works, let's generate some realistic looking data. I have chosen a true date of January 30 (day of year is 30) and a bias of +5 days.

set.seed(17)
day.actual <- 30 # True date
bias <- 5
sigma <- 25      # Nominal dispersion (of the normal values used for data generation)
n <- 1e3         # Dataset size
U <- 183 * (runif(n) < 0.05) # Used to generate some outlying values
day.simulated <- round(day.actual + bias + rnorm(n, U, sigma)) %% 365
hist(day.simulated, breaks = seq(0, 365, length.out = 51),
main = "Simulated Data", xlab = "Day of year")


Here is a summary of a thousand such values:

On this histogram I have drawn the circular mean as a red dashed line and the true reference date as a blue line, illustrating the use of summary.circular:

stats <- summary.circular(day.simulated, 365)
abline(v = stats[1], lwd = 2, col = "Red", lty = 3)


The value in stats for these data is

      Mean Dispersion
35.33973   38.02667


Because the input data are in days, these statistics are in days, too. Their bias is $$35.33973 - 30 = 5.33973,$$ differing from the intended value of $$5$$ only because of the random generation of these data.

• This is fantastic! Thank you for the thorough explanation! I will put it to the test :) Commented Aug 3, 2023 at 20:19

You probably want statistics on when the observed events $$E$$ happened relative to a reference point $$R$$. If so, this can be calculated as:

((E-R+182.5) %% 365) - 182.5


where %% is the function known mathematically as $$\text{mod}$$. This formula produces a number of days between $$-182$$ and $$+182$$, negative for events up to half a year before $$R$$, and positive for events up to half a year after $$R$$. All the usual statistics can be calculated based on these differences.

It probably does not make sense to use the mean of your data as a reference point, precisely because means of circular statistics are not obvious.

A good possibility for a reference point could be:

• a holiday
• a solstice or equinox
• the first day of a month
• the first day of a hunting season
• the mode of your data, calculated as names(sort(-table(E)))[1].

Any of these can work for $$R$$, so long as there are few events which are half a year from $$R$$, and so long as a phrase like “five days before $$R$$” will make intuitive sense.

EDIT In the example above, you can:

• Use the modal happiest day in each country as its reference point,
• e.g.: Christmas.
• Convert calendar days to signed integers using the first formula above,
• e.g.: two days before Christmas or three days after Christmas.
• Calculate statistics on those integers as usual,
• e.g.: a mean of 4 days before Christmas and sd of 5 days.
• Convert the statistics back to calendar days,
• e.g.: a mean of December 21, and a standard deviation of 5 days.
• Oh, that is a really interesting idea! The use of a reference point could be useful in ways that go beyond the scope of this question. However, your last point would still be an issue--if there are points that are a half year from R. In actuality, the distribution I have plotted in my post are draws from a bootstrap procedure that is trying to show the uncertainty of my doy estimate. I am doing this procedure for many different types of events, and I want to be able to report them succinctly. Hence why I want to get the mean/median and a CI of these values. Some may hit the 1/2 year... Commented Aug 3, 2023 at 14:37
• If, indeed, you are researching the uncertainty (variability) in an estimate of a circular variable, then you need to change your question, because the answer is to summarize that variability in a circular metric, and not just to describe the distribution of the results without reference to the estimand.
– whuber
Commented Aug 3, 2023 at 15:37
• I tried to edit the question a bit to help. What I want to do is provide the mean/median of these bootstrap estimates of doy while also providing an estimate of the uncertainty. Commented Aug 3, 2023 at 16:37