2
$\begingroup$

I have some time-series data showing the monthly counts of hospital admissions. It has both a long-term trend (increasing) and seasonality (highest in summer). I am trying to measure the magnitude of the seasonality (i.e. how much greater is the rate of admission in summer than in winter?).

Here is some sample data (this is R code):

ts <- structure(list(year = c(2010L, 2010L, 2010L, 2010L, 2010L, 2010L, 
                              2010L, 2010L, 2010L, 2010L, 2010L, 2010L, 2011L, 2011L, 2011L, 
                              2011L, 2011L, 2011L, 2011L, 2011L, 2011L, 2011L, 2011L, 2011L, 
                              2012L, 2012L, 2012L, 2012L, 2012L, 2012L, 2012L, 2012L, 2012L, 
                              2012L, 2012L, 2012L, 2013L, 2013L, 2013L, 2013L, 2013L, 2013L, 
                              2013L, 2013L, 2013L, 2013L, 2013L, 2013L, 2014L, 2014L, 2014L, 
                              2014L, 2014L, 2014L, 2014L, 2014L, 2014L, 2014L, 2014L, 2014L), month = structure(c(1L, 2L, 3L, 4L, 5L, 6L, 7L, 8L, 9L, 10L, 
                       11L, 12L, 1L, 2L, 3L, 4L, 5L, 6L, 7L, 8L, 9L, 10L, 11L, 12L, 
                       1L, 2L, 3L, 4L, 5L, 6L, 7L, 8L, 9L, 10L, 11L, 12L, 1L, 2L, 3L, 
                       4L, 5L, 6L, 7L, 8L, 9L, 10L, 11L, 12L, 1L, 2L, 3L, 4L, 5L, 6L, 
                       7L, 8L, 9L, 10L, 11L, 12L), .Label = c("Jan", "Feb", "Mar", "Apr", 
                                                              "May", "Jun", "Jul", "Aug", "Sep", "Oct", "Nov", "Dec"), class = "factor"), period = 1:60, N = c(539L, 528L, 584L, 628L, 644L, 670L, 
                     657L, 680L, 629L, 660L, 614L, 606L, 564L, 699L, 623L, 621L, 
                     692L, 683L, 644L, 715L, 639L, 646L, 625L, 609L, 547L, 594L, 
                     653L, 682L, 724L, 737L, 671L, 667L, 698L, 688L, 621L, 643L, 
                     573L, 643L, 675L, 722L, 669L, 659L, 737L, 692L, 647L, 701L, 
                     653L, 603L, 564L, 638L, 657L, 681L, 720L, 713L, 674L, 788L, 
                     697L, 717L, 680L, 648L)), row.names = c(NA, -60L), class = "data.frame")

We can fit Poisson models with and without seasonality, and plot the predicted values:

model1 <- glm(N ~ period, family = 'poisson', data = ts)
model2 <- glm(N ~ period + month, family = 'poisson', data = ts)

ts$pred1 <- predict(model1, newdata = ts, type = 'response')
ts$pred2 <- predict(model2, newdata = ts, type = 'response')

plot(1, type = 'n', xlim = c(0, 60), ylim = c(0, max(ts$N) * 1.2), xlab = 'time', ylab = NA)
with(ts, {
  points(period, N)
  lines(period, pred1, col = 'blue')
  lines(period, pred2, col = 'red')
})

Time series

Clearly there is strong evidence of seasonality:

anova(model1, model2, test = 'LRT')

Analysis of Deviance Table

Model 1: N ~ period
Model 2: N ~ period + month
  Resid. Df Resid. Dev Df Deviance  Pr(>Chi)    
1        58    202.753                          
2        47     62.334 11   140.42 < 2.2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

However, how do we measure the magnitude of this seasonality (rather than the statistical evidence for seasonality or goodness-of-fit of the model)?

My first thought was to calculate the extent to which residuals are reduced by adding a term for seasonality into the model, in comparison to the data.

# calculate residuals:
ts$resid1 <- ts$pred1 - ts$N
ts$resid2 <- ts$pred2 - ts$N

# calculate the reduction in residuals in comparison to the data
(sum(abs(ts$resid1)) - sum(abs(ts$resid2))) / sum(ts$N)

This suggests there is 2.7% seasonality in this data. Does this make sense? Is there a standard approach to this problem?

$\endgroup$
4
  • 1
    $\begingroup$ It might be more "statistically natural" to consider the reduction in squared residuals, similar to an $F$ statistic. The forecaster in me would use a holdout sample and forecast into it with a seasonal and a non-seasonal model (e.g., exponential smoothing), and report the percentage reduction in the mean squared error. ... $\endgroup$ Jun 24 at 12:26
  • 1
    $\begingroup$ ... Whether in-sample or out-of-sample, this is a case of finding out how much including a predictor improves a model - there should definitely be something on this, but I wasn't able to find it in a short search here. Related: Importance of regressors in time series data $\endgroup$ Jun 24 at 12:27
  • 1
    $\begingroup$ You might find this paper useful. They used similar methods to measure the amplitude of seasonality. Link: sciencedirect.com/science/article/pii/S1198743X14610934 $\endgroup$ Jun 25 at 7:10
  • $\begingroup$ Thanks @SpurEconomics, that was a very useful paper. I used the idea of a "peak-to-low" ratio and used a Monte-Carlo method to estimate it - see answer below. $\endgroup$
    – Dan
    Jun 26 at 18:27

1 Answer 1

0
$\begingroup$

Here's one solution I've tried for estimating the peak-to-low ratio. My main problem was that the month variable in the Poisson model had an arbitrary reference month (eg. January) and you might not know which month was the highest (only the month with the highest point estimate). Therefore the 'biggest' association from the Poisson model didn't seem the right solution. I used the Poisson model to predict the number of admissions in a 'typical' January, February, etc., and calculated the ratio of peak-to-low months. I repeated this for 1000 simulated datasets in which the count of admissions for each month was sampled from a Poisson distribution:

nd <- data.frame(mth = month.abb, p = 1)
p <- ts$period
mth <- ts$month

set.seed(10)
bs <- sapply(1:1000, function (x) {
  if (x %% 100 == 0) print(x)
  a <- rpois(nrow(ts), ts$N)
  m <- glm(a ~ p + mth, family = 'poisson')
  b <- predict(m, newdata = nd, type = 'response')
  max(b) / min(b)
})
quantile(bs, probs = c(0.5, 0.025, 0.975))

This gives a peak-to-low ratio and confidence interval:

     50%     2.5%    97.5% 
1.260151 1.209427 1.322026 
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.