How to add a confidence interval to a Circular Histogram (von Mises distribution)

user_5810852 0
Viewed 409

I have time data, and I want to plot the frequency per hour on a 24hr clock.

The data are transformed to circular, and the estimates for 'periodic mean' mu and 'concentration' kappa are calculated with mle.vonmises().

The graph is generated with ggplot2, using geom_hist() and coord_polar(). The periodic mean is drawn on the plot with a simple call to geom_vline().

Question

I want to draw a confidence interval of 95% around the mean. Then, I would like to visually check whether a given timestamp (e.g. "22:00:00") lies within the CI or not. How do I do this with a von mises distribution and ggplot2?

The code below shows how far I got.

The data

timestamps <- c("08:43:48", "09:17:52", "12:56:22", "12:27:32", "10:59:23", 
                "07:22:45", "11:13:59", "10:13:26", "10:07:01", "06:09:56", 
                "12:43:17", "07:07:35", "09:36:44", "10:45:00", "08:27:36", 
                "07:55:35", "11:32:56", "13:18:35", "11:09:51", "09:46:33", 
                "06:59:12", "10:19:36", "09:39:47", "09:39:46", "18:23:54")

The code

library(lubridate)
library(circular)
library(ggplot2)

## Convert from char to hours
timestamps_hrs <- as.numeric(hms(timestamps)) / 3600

## Convert to class circular
timestamps_hrs_circ <- circular(timestamps_hrs, units = "hours", template = "clock24")

## Estimate the periodic mean and the concentration 
## from the von Mises distribution
estimates <- mle.vonmises(timestamps_hrs_circ)
periodic_mean <- estimates$mu %% 24
concentration <- estimates$kappa

## Clock plot // Circular Histogram
clock01 <- ggplot(data.frame(timestamps_hrs_circ), aes(x = timestamps_hrs_circ)) +
  geom_histogram(breaks = seq(0, 24), colour = "blue", fill = "lightblue") +
  coord_polar() + 
  scale_x_continuous("", limits = c(0, 24), breaks = seq(0, 24), minor_breaks = NULL) +
  theme_light()

clock01

## Add the periodic_mean
clock01 + 
  geom_vline(xintercept = as.numeric(periodic_mean), color = "red", linetype = 3, size = 1.25)

This yields the following graph:

enter image description here

confidence-interval ggplot2 r
1 Answers

I think I found an approximation of the solution. As we know the parameters mu and kappa (resp. the periodic mean and the concentration), we know the distribution. This, in turn, means we know the densities of given timestamps, and we can calculate the cutoff for the 95% confidence level.

Once we have that, we can generate the timestamps for every minute of a day. We transform the timestamps as needed, calculate the densities, and compare against the cutoff value.

This way we know on the 1 minute level, whether we are in the Confidence Interval or not.

The code

(it is assumed that the code in the question has been run)

quantile <- qvonmises((1 - 0.95)/2, mu = periodic_mean, kappa = concentration)
cutoff <- dvonmises(quantile, mu = periodic_mean, kappa = concentration)

## generate a timestamp for every minute in a day
## then the transformations needed
ts_1min <- format(seq.POSIXt(as.POSIXct(Sys.Date()), 
                             as.POSIXct(Sys.Date()+1), 
                             by = "1 min"), 
                  "%H:%M:%S", tz = "GMT")
ts_1min_hrs <- as.numeric(hms(ts_1min)) / 3600
ts_1min_hrs_circ <- circular(ts_1min_hrs, units = "hours", template = "clock24")
## generate densities to compare with the cutoff
dens_1min <- dvonmises(ts_1min_hrs_circ, mu = periodic_mean, kappa = concentration)
 
## compare: vector of FALSE/TRUE
feat_1min <- dens_1min >= cutoff
df_1min_feat <- data.frame(ts = ts_1min_hrs_circ, 
                             feature = feat_1min)

## get the min and max time of the CI
CI <- df_1min_feat %>% 
  filter(feature == TRUE) %>%
  summarise(min = min(ts), max= max(ts))

CI
#   min      max
# 5.283333 14.91667

With the information above, and using geom_rect(), we can get what we want:

ggplot(data.frame(timestamps_hrs_circ), aes(x = timestamps_hrs_circ)) +
  coord_polar() +
  scale_x_continuous("", limits = c(0, 24), breaks = seq(0, 24), minor_breaks = NULL) +
  geom_vline(xintercept = as.numeric(CI), color = "darkgreen", linetype = 1, size = 1.5) +
  geom_rect(xmin = CI$min, xmax = CI$max, ymin = 0, ymax = 5, alpha = .5, fill = "lightgreen") +
  ggtitle(label = "Circular Histogram", subtitle = "periodic mean in red,\n95%-CI in green" ) +
  geom_histogram(breaks = seq(0, 24), colour = "blue", fill = "lightblue") +
  geom_vline(xintercept = as.numeric(periodic_mean), color = "red", linetype = 2, size = 1.5) +
  theme_light()

Resulting in the following graph:

enter image description here

I hope someone may benefit from this as well.

user_58108520
Related