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?
 A: 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 ;)
A: 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.
A: 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.
A: 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.
A: You might consider applying a Tukey Control chart to the data.
A: 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.
A: 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.
A: 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.
A: 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.
