Time series for count data, with counts < 20

I recently started working for a tuberculosis clinic. We meet periodically to discuss the number of TB cases we're currently treating, the number of tests administered, etc. I'd like to start modeling these counts so that we're not just guessing whether something is unusual or not. Unfortunately, I've had very little training in time series, and most of my exposure has been to models for very continuous data (stock prices) or very large numbers of counts (influenza). But we deal with 0-18 cases per month (mean 6.68, median 7, var 12.3), which are distributed like this:

[image lost to the mists of time]

[image eaten by a grue]

I've found a few articles that address models like this, but I'd greatly appreciate hearing suggestions from you - both for approaches and for R packages that I could use to implement those approaches.

EDIT: mbq's answer has forced me to think more carefully about what I'm asking here; I got too hung-up on the monthly counts and lost the actual focus of the question. What I'd like to know is: does the (fairly visible) decline from, say, 2008 onward reflect a downward trend in the overall number of cases? It looks to me like the number of cases monthly from 2001-2007 reflects a stable process; maybe some seasonality, but overall stable. From 2008 through the present, it looks like that process is changing: the overall number of cases is declining, even though the monthly counts might wobble up and down due to randomness and seasonality. How can I test if there's a real change in the process? And if I can identify a decline, how could I use that trend and whatever seasonality there might be to estimate the number of cases we might see in the upcoming months?

• The links have died. If you're able to get hold of the images (or regenerate them), please use the new system which stores them at imgur. – Glen_b Mar 2 '15 at 1:56
• Unfortunately, these plots were from a couple of jobs ago. Apologies! – Matt Parker Mar 9 '15 at 18:28
• Given this post (#173) is from way back when the site was in beta, that's not a surprise - few people could redo a plot from posts back that far at this point. Thanks anyway. – Glen_b Mar 10 '15 at 0:55
• Now your images have turned into ad banners. – Ilmari Karonen Aug 25 '15 at 20:03

To assess the historical trend, I'd use a gam with trend and seasonal components. For example

require(mgcv)
require(forecast)
x <- ts(rpois(100,1+sin(seq(0,3*pi,l=100))),f=12)
tt <- 1:100
season <- seasonaldummy(x)
fit <- gam(x ~ s(tt,k=5) + season, family="poisson")
plot(fit)


Then summary(fit) will give you a test of significance of the change in trend and the plot will give you some confidence intervals. The assumptions here are that the observations are independent and the conditional distribution is Poisson. Because the mean is allowed to change smoothly over time, these are not particularly strong assumptions.

To forecast is more difficult as you need to project the trend into the future. If you are willing to accept a linear extrapolation of the trend at the end of the data (which is certainly dodgy but probably ok for a few months), then use

fcast <- predict(fit,se.fit=TRUE,
newdata=list(tt=101:112,season=seasonaldummyf(x,h=12)))


To see the forecasts on the same graph:

plot(x,xlim=c(0,10.5))
lines(ts(exp(fcast$fit),f=12,s=112/12),col=2) lines(ts(exp(fcast$fit-2*fcast$se),f=12,s=112/12),col=2,lty=2) lines(ts(exp(fcast$fit+2*fcast\$se),f=12,s=112/12),col=2,lty=2)


You can spot the unusual months by looking for outliers in the (deviance) residuals of the fit.

You might want to have a look at strucchange:

Testing, monitoring and dating structural changes in (linear) regression models. strucchange features tests/methods from the generalized fluctuation test framework as well as from the F test (Chow test) framework. This includes methods to fit, plot and test fluctuation processes (e.g., CUSUM, MOSUM, recursive/moving estimates) and F statistics, respectively. It is possible to monitor incoming data online using fluctuation processes. Finally, the breakpoints in regression models with structural changes can be estimated together with confidence intervals. Emphasis is always given to methods for visualizing the data."

PS. Nice graphics ;)

• I'll have to read it more thoroughly later on, but yes, this package is definitely addressing the kinds of problems I'm facing here. Thanks! And also, thanks for the kind words about the plots ;p – Matt Parker Jul 21 '10 at 15:16

Does it really need some advanced model? Based on what I know about TB, in case there is no epidemy the infections are stochastic acts and so the count form month N shouldn't be correlated with count from month N-1. (You can check this assumption with autocorrelation). If so, analyzing just the distribution of monthly counts may be sufficient to decide if some count is significantly higher than normal.
On the other hand you can look for correlations with some other variables, like season, travel traffic, or anything that you can imagine that may be correlated. If you would found something like this, it could be then used for some data normalization.

• Your point about month N's count not necessarily being correlated with N-1 is well-taken. With a slow-growing disease like TB, that's something I'd have to look at carefully, but I'm pretty sure I could identify about how much lag there is between the time we report a source case and the time we report any secondary cases. – Matt Parker Jul 21 '10 at 0:11
• However, it's your point about analyzing the distribution of monthly counts that's at the heart of my question. There is a definite decline in TB, both nationally in the US and in my district. For example, when I compare 2009 to the previous years, there are decidedly fewer cases. 2010 is on track to have fewer still. What I'm trying to identify (which I did a poor job of explaining in the question) is whether or not these declines are part of an ongoing downward trend, or just a downward wobble. Thanks - you've gotten me to think much more carefully about the problem. – Matt Parker Jul 21 '10 at 0:15

Often, disease data like this is performed with a generalized linear model, as its not necessarily a great application of time series analysis - the months often aren't all that correlated with each other.

If I were given this data, here's what I would do (and indeed, have done with data similar to it):

Create a "time" variable that's more accurately described as "Months since 1/1/2000" if I'm eyeballing your data correctly. Then I'd run a general linear model in R using the Poisson distribution (or Negative Binomial) and a log link with roughly the following form:

log(Counts) = b0 + b1*t + b2*(t^2) + b3*cos(2pi*w*t) + b4*sin(2pi*w*t)

Where t is the time described above, and w is 1/365 for a yearly disease like flu. Generally its 1/n, where n is the length of your disease's cycle. I don't know offhand what it is for TB.

The two time trends will show you - outside normal seasonal variation - if you have meaningful variation over time.

You might consider applying a Tukey Control chart to the data.

• Heh - this was actually one of the plots I made that didn't make it into the post. The problem I had is deciding how to calculate the bounds - my initial attempt was with Poisson bounds, with lambda set to the mean of my data, but variance is too high to be a proper Poisson (high enough to matter in practice? I don't know). – Matt Parker Jul 21 '10 at 15:01
• A further problem is that the center of the distribution can change over time - for example, it wouldn't make sense to set those bounds using data from the early 1900s, when Colorado was a haven for TB patients. So what's an appropriate way to keep the lines up-to-date with long-term changes in the process, while still being able to identify deviations? – Matt Parker Jul 21 '10 at 15:02

You may try to model your data using a Dynamic Generalized Linear Model (DGLM). In R, you can fit this kind of models using packages sspir and KFAS. In a sense, this is similar to the gam approach suggested by Rob, except that instead of assuming that the log mean of the Poisson observations be a smooth function of time, it assumes that it follows a stochastic dynamics.

• Cool - I didn't know about the KFAS package. There's also the dlm and dse for state space-ish approaches, and an general overview for R users here: cran.r-project.org/web/views/TimeSeries.html – conjugateprior Nov 2 '10 at 13:19
• I'd HIGHLY recommend the dlm package. DLMs are not as easy to use as other solutions, but dlm makes it as easy as possible and has a nice vignette walking you through the process. – Wayne Sep 22 '11 at 16:59

I'm going to leave the main question alone, because I think I will get it wrong (although I too analyse data for a healthcare provider, and to be honest, if I had these data, I would just analyse them using standard techniques and hope for the best, they look pretty okay to me).

As for R packages, I have found the TSA library and it's accompanying book very useful indeed. The armasubsets command, particularly, I think is a great time saver.

Escape from traditional enumerative statistics as Deming would suggest and venture into traditional analytical statistics - in this case, control charts. See any books by Donald Wheeler PhD, particularly his "Advanced Topics in SPC" for more info.

• Someone just asked a question about SPC (aka QC) charts in R in stats.stackexchange.com/questions/15809/… where I provides some package hints. I'm not sure of SPC/QC myself: I wonder if it was more useful in the pre-computer era and for workers on the shop floor, but it's worth considering. – Wayne Sep 22 '11 at 16:57
• Actually, is this answer redundant with @babelproofreader's? – Wayne Sep 22 '11 at 17:00

In response to your direct question "How can I test if there's a real change in the process? And if I can identify a decline, how could I use that trend and whatever seasonality there might be to estimate the number of cases we might see in the upcoming months?" Develop a Transfer Function Model ( ARMAX ) that readily explains period-to-period dependency including and seasonal ARIMA structure. Incorporate any Identifiable Level Shifts , Seasonal Pulses, Local Time Trends and PUlses that may have been suggested by empirical/analystical methods like Intervention Detection. IF THIS ROBUST MODEL INCLUDES A FACTOR/SERIES matching up with "declines" Then your prayers have been answerered. In the alternative simply add an hypothesized structure e.g. to test a time trend change at point T1 construct two dummies X1 = 1,1,2,3,,,,,,T and X2 = 0,0,0,0,0,0,0,1,2,3,4,5,.... WHERE THE ZEROES END AT PERIOD T1-1 . The test of the hypothesis of a significant trend change at time period T1 will be assessed using the "t value" for X2 .

Edited 9/22/11

Often, disease data like this has monthly effects since weather/temperature is often an unspecified causal . In the omission of the true caudsal series ARIMA models use memory or seasonal dummies as a surrogate. Additionally series like this can have level shifts and/or local time trends reflecting structural change over time. Exploiting the autoregressive structure in the data rather than imposing various artifacts like time and time square and time cubic etc. have been found to be quite useful and less presumptive and ad hoc. Care should also be taken to identify "unusual values" as they can often be useful in suggestng additional cause variables and at a minimum lead to robust estimates of the other model parameters. Finally we have found that the variability/paramaters may vary over times thus these model refinements may be in order.

• How does one adjust the p-value to account for identifying T1 post hoc (i.e., noticing it after reviewing the data)? – whuber May 8 '11 at 22:29
• all refinements i.e. diagnostic checking leading to augmentation are treated the same as structure suggested while the data was in the drawer. Adding a lag to a predictor variable based upon diagnostic checking simply adds another null hypothesis to the list. This doesn't differ from deleting a non-significant coefficent. – IrishStat May 9 '11 at 2:20
• @Irish I don't follow that. T1 doesn't seem to be a lag: it stipulates a changepoint. Isn't it easy to find "significant" changes if you're allowed first to inspect the data to choose where the changepoint is? – whuber May 9 '11 at 3:43
• @whuber Yes you are right. But think of a Stepwise Forward Process in Regression where you examine the alternative (omitted) variables for "potential incorporation" . This is no different in that sense that you are deciding which "omitted trend" should be added to the model in order to render the resultant error process Gaussian. – IrishStat May 9 '11 at 10:56
• @Irish That's an illuminating analogy. If I understand, in effect you contemplate having one potential variable for every period (designating a possible change at that period) and are invoking a systematic process to determine which of these should be included in the model. This suggests that some standard p-level adjustment procedures, such as the Bonferroni, might reasonably be applied. Would that be valid? – whuber May 9 '11 at 14:02