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?