I probably wouldn't do this with a hypothesis test, but we can use some probabilistic reasoning here.
Poisson regression is a way to sort of "fit" a distribution on various covariates (like day of week for example). So here is what I propose.
- Use a couple months worth of data to fit a poisson regression
- Estimate prediction intervals via the bootstrap, and
- When an observation falls outside of these intervals, declare it as anomalous.
I'll illustrate with a simulated example. Suppose you're working with data that looks like the plot below. The data runs between Oct 1 2023 and Dec 31 2023, and let's say you're looking to make predictions on December.. Clearly, there is an anomaly sometime during December. I've purposefully made this very large for purposes of illustration.
Our first step is to fit some sort of model. There is clearly some sort of periodic behaviour here, so a good first attempt would be to fit a regression on day of week.
# Fit a model
training <- filter(d, dts<'2023-12-01')
fit <- glm(y ~ dow, data=training, family = poisson())
Now, we can plot our predictions to see how the model does. I'll focus just on December for now. The data (black) and predictions (red) look like
Predictions look good, but how can we tell if the errors in the predictions are normal versus anomalous. We need some sense of what is to be expected. Here is where we use the bootstrap to integrate out the uncertainty.
The sampling distribution of the coefficients is something like
$$\hat{\beta} \sim \mbox{Normal}(\beta, \Sigma)$$
R provides us with an estimate of $\Sigma$ via the vcov method, we can act as if $\beta = \hat{\beta}$ for now. So we can draw a bunch of samples from this distribution, compute our predictions, and sample from the likelihood to get an idea of what kind of observations we should expect to see.
The outcome is something like this
The red ribbon is the range where we expect to see 95% of observations from the process, assuming the model is correct. So one way to operationalize an outlier would be to say "if it is outside this band, then it has low probability of happening". As you can see, out outlier is above this point, so it would be deemed as an outlier.
There are various ways this could be improved (e.g. by bootstrapping the model fit within the computation of the prediction intervals), but I think this is fine for now.
Code
library(tidyverse)
# Create a dater range
dts <- seq(
ymd('2023-10-01'),
ymd('2023-12-31'),
by='day'
)
d <- tibble(dts) %>%
mutate(
dow = wday(dts, label=T)
)
# Simulate the data
X <- model.matrix(~dow-1, data = d)
beta <- log(c(200, 300, 400, 400, 400, 400, 200))
y <- rpois(nrow(X), exp(X%*%beta))
y[80] <- round(y[80] * 1.5)
d$y <- y
plot <- d %>%
ggplot(aes(dts, y)) +
geom_line() +
theme(aspect.ratio = 1/1.61)
# Fit a model
training <- filter(d, dts<'2023-12-01')
fit <- glm(y ~ dow, data=training, family = poisson())
d$pred <- predict(fit, newdata = d, type='response')
plot +
geom_line(data=d, aes(dts, pred), color='red') +
scale_x_date(limits = c(ymd('2023-12-01'), NA))
get_prediction_intervals <- function(fit){
samples <- map_dfr(1:1000, ~{
b <- coef(fit)
Sigma <- vcov(fit)
bsamp <- MASS::mvrnorm(n=1, mu=b, Sigma = Sigma)
# Make predictions for each day of week
Xd <- fit$data %>%
distinct(dow)
lambda <- predict(fit, newdata = Xd, type = 'response')
Xd$yobs <- rpois(length(lambda), lambda)
Xd
}, .id = 'sample')
samples %>%
group_by(dow) %>%
reframe(
prob = c(0.025, 0.975),
bnd = c('pred.lo', 'pred.hi'),
pct = quantile(yobs, probs = c(0.025, 0.975))
) %>%
pivot_wider(id_cols = c('dow'), names_from = 'bnd', values_from = 'pct')
}
ints <- get_prediction_intervals(fit=fit)
plot +
geom_line(data=d, aes(dts, pred), color='red') +
geom_ribbon(data=left_join(d, ints), aes(dts, pred, ymin=pred.lo, ymax=pred.hi), alpha=0.5, fill='red', inherit.aes = F) +
scale_x_date(limits = c(ymd('2023-12-01'), NA))
```
