Using a Negative Binomial model with smooth splines to predict infectious disease case counts
Nick Reich, adapted from work by Steve Lauer et al
Introduction
This analysis is based on research that is publicly available on GitHub. The code and description have been adapted for this write-up. Here is are some expository excerpts from the paper:
Dengue, a mosquito-borne virus prevalent throughout the tropics and sub-tropics, infects an estimated 390 million people every year. While the majority of infections are mild or asymptomatic, the more severe forms of dengue infection – dengue shock syndrome (DSS) and dengue hemorrhagic fever (DHF) – can result in organ failure or death.
… In Thailand, dengue infection is endemic with substantial annual and geographic variation in incidence across its 76 provinces and 13 health regions. Over the past 15 years, an average of 43,137 (range 14,952-106,320) DHF cases have been reported to the Thailand Ministry of Public Health (MOPH) each year. Within a typical year, incidence rates in different provinces can vary by an order of magnitude, with some provinces experiencing less than 10 DHF cases per 100,000 population and others over 100 per 100,000 population.
… Whether an infectious disease spreads within a population depends on the transmission rate of the disease and the number of susceptible individuals, thus long-term forecasting models for DHF incidence may need to account for climatic factors that could affect transmission as well as population susceptibility. Climatic factors, such as temperature, rainfall, and humidity, may impact both the prevalence and distribution of the dengue vector, the mosquito, as well as the transmission efficiency of dengue virus. Such climatic factors vary significantly across Thailand.
… We obtained data on DHF cases (from the MOPH), population (National Statistical Office of Thailand), and weather (NOAA). These data were summarized across time frames ranging from one month to one year to create covariates for consideration by our model selection algorithm. … We made forecasts using the data available in April of each year, the month when the MOPH has historically finalized the incidence reports obtained from all provinces for the prior calendar year. Hence, all ``annual’’ forecasts are for DHF incidence between April and December of the year they are made.
… For the NOAA and NCDC weather station data, we found the most consistently reported weather station for each province and extracted the daily maximum and minimum temperature, maximum humidity, and rainfall. We aggregated these measures into monthly covariates for maximum, minimum, and mean temperature, maximum and mean humidity, and maximum and total rainfall across the low-season, January, February, and March.
The data
The dataset shown here and provided in this repository is the same one used in the analysis. Some columns have been removed for simplicity. Here, we split the dataset into a training dataset from 2000 through 2012 and a testing dataset from 2013 and 2014.
dat <- read.csv("dengue-data.csv")
train_dat <- filter(dat, year<2013) ## all but 2013 and 2014
test_dat <- filter(dat, year>=2013)
str(train_dat)'data.frame': 988 obs. of 27 variables:
$ pid : Factor w/ 76 levels "TH01","TH02",..: 1 1 1 1 1 1 1 1 1 1 ...
$ year : int 2000 2001 2002 2003 2004 2005 2006 2007 2008 2009 ...
$ obs : int 1 63 28 35 49 18 13 61 92 173 ...
$ population : int 210537 210399 210260 210121 209983 209844 209706 209568 209430 209291 ...
$ population_density: num 16.6 16.6 16.6 16.6 16.6 ...
$ pre_season_rate : num 0 0.475 0 0.476 0.476 ...
$ last_post_rate : num 0 0 0.951 1.427 2.38 ...
$ last_high_rate : num 22.8 0 29 11.9 14.3 ...
$ sus_pct : num 17.1 36.6 29.4 16.1 28.9 ...
$ mean_temp_grid : num 22.5 23.9 23.1 23.2 23.2 ...
$ total_rain_grid : num 213.4 95.7 69.6 149 60.1 ...
$ Jan_temp_grid : num 22.1 22.1 21.9 21.1 21.3 ...
$ Feb_temp_grid : num 23.8 25.1 24 24.4 25.2 ...
$ Mar_temp_grid : num 26.8 28.6 27.2 28.6 28.6 ...
$ Jan_rain_grid : num 34.15 2.02 14.1 9.32 5.38 ...
$ Feb_rain_grid : num 42.179 71.425 7.046 26.814 0.958 ...
$ Mar_rain_grid : num 125.8 12.5 30 59.1 48.7 ...
$ max_humid_ncdc : num 91.7 94.6 91.6 94.4 92.8 ...
$ min_humid_ncdc : num 42.1 47 41.9 52.8 38.9 ...
$ Jan_humid_ncdc : num 93.6 96.4 94.6 96.4 98.9 ...
$ Feb_humid_ncdc : num 89.8 94 91.2 94.4 91.4 ...
$ Mar_humid_ncdc : num 84.4 88.1 82.3 89 80.8 ...
$ max_temp_ncdc : num 30.8 31.7 31.3 30.7 32 ...
$ min_temp_ncdc : num 17.1 17.8 17.4 18.6 17 ...
$ Jan_temp_ncdc : num 30.5 30.8 28.8 27.8 30.4 ...
$ Feb_temp_ncdc : num 32 33.6 33.2 31.6 32.8 ...
$ Mar_temp_ncdc : num 34.4 34.3 36.3 35.6 37.1 ...ggplot(dat, aes(x=year, y=reorder(pid, population), fill=obs/population)) + geom_tile() +
scale_fill_gradient(low="grey100", high="firebrick") +
ylab("Province (by population, large at top)") + xlab(NULL) 
Before we begin modeling, let’s see if there is evidence that we will need to compensate for overdispersion in a Poisson model.
disp_dat <- train_dat %>%
group_by(pid) %>%
summarize(
avg_inc_rate = mean(obs),
var_inc_rate = var(obs),
var_over_avg = var_inc_rate/avg_inc_rate
)
quantile(disp_dat$var_over_avg) 0% 25% 50% 75% 100%
5.489567 108.919350 192.394499 345.826018 1724.248223 There sure is. The per-province variance ranges from 7 to over 1500 times the accompanying mean.
Modeling
Negative binomial model
The model that we used to forecast annual DHF incidence for this study is a generalized additive model . Specifically, we use a generalized additive model with a negative binomial family, separate penalized smoothing splines for each covariate, and province-level random effects:
\[ Y_{i,t} \sim \text{NB}(n_{i,t}\lambda_{i,t}, r) \] \[ \log \Big [ \mathbf{E}(Y_{i,t}) \Big ] = \beta_0 + \log(n_{i,t}) + \alpha_i + \sum_{j=1}^J g_j(x_{j,i,t}|\boldsymbol{\theta}) \] \[ \alpha_i \sim \text{Normal}(\mu, \sigma^2) \]
We model the incidence (\(Y_{i,t}\)) for province \(i\) in year \(t\) as following a negative binomial distribution with the mean equal to the province population (\(n_{i,t}\)) times the incidence rate (\(\lambda_{i,t}\)) and a dispersion parameter \(r\). After a log transformation, we model the mean of this distribution using an intercept (\(\beta_0\)), a random effect for each province (\(\alpha_i\)) and a cubic spline for each of \(J\) covariates (\(g_j(x_{j,i,t}|\boldsymbol{\theta})\)).
Let’s fit a model and look at the results.
library(mgcv)
fmla <- formula("obs ~ offset(log(population)) +
s(pid, bs = 're') +
s(pre_season_rate, k = 3, m=2, bs = 'cr') +
s(Jan_temp_grid, k = 3, m=2, bs = 'cr') +
s(last_high_rate, k = 3, m=2, bs = 'cr')")
mod_nb <- gam(fmla, data = train_dat, family = nb())
summary(mod_nb)
Family: Negative Binomial(2.74)
Link function: log
Formula:
obs ~ offset(log(population)) + s(pid, bs = "re") + s(pre_season_rate,
k = 3, m = 2, bs = "cr") + s(Jan_temp_grid, k = 3, m = 2,
bs = "cr") + s(last_high_rate, k = 3, m = 2, bs = "cr")
Parametric coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -7.64116 0.03424 -223.1 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Approximate significance of smooth terms:
edf Ref.df Chi.sq p-value
s(pid) 50.403 75.000 144.64 2.07e-15 ***
s(pre_season_rate) 1.991 1.999 539.05 < 2e-16 ***
s(Jan_temp_grid) 1.953 1.990 22.20 3.04e-05 ***
s(last_high_rate) 1.774 1.944 23.23 1.75e-05 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
R-sq.(adj) = 0.776 Deviance explained = 73.5%
-REML = 6523.7 Scale est. = 1 n = 988The random effects term is significant but notice that the degrees of freedom is fairly high relative to the number of provinces (76), indicating quite a bit of variation between the provinces and not a lot of shrinkage or pooling of information. It also looks like all of our chosen covariates are “highly significant”, but because these relationships are modeled as smooth functions of the covariates, there are no coefficients to interpret. What do these relationships look like?
par(mfrow=c(2,2))
plot(mod_nb)
We can also fit some reduced models to compare deviance:
mod_ref <- gam(obs ~ offset(log(population))+1, data = train_dat, family = nb)
mod_randint <- gam(obs ~ offset(log(population)) + s(pid, bs = 're'), data = train_dat, family = nb)
anova(mod_ref, mod_randint, mod_nb)Analysis of Deviance Table
Model 1: obs ~ offset(log(population)) + 1
Model 2: obs ~ offset(log(population)) + s(pid, bs = "re")
Model 3: obs ~ offset(log(population)) + s(pid, bs = "re") + s(pre_season_rate,
k = 3, m = 2, bs = "cr") + s(Jan_temp_grid, k = 3, m = 2,
bs = "cr") + s(last_high_rate, k = 3, m = 2, bs = "cr")
Resid. Df Resid. Dev Df Deviance
1 987.00 1097.24
2 914.41 1018.17 72.5871 79.063
3 912.85 997.45 1.5616 20.719And here is a Poisson model for comparison. It’s interesting to see how the Poisson model has quite a bit less uncertainty in the estimated relationships.
mod_pois <- gam(fmla, data = train_dat, family = poisson)
par(mfrow=c(2,2))
plot(mod_pois)
Obtaining predictive distribution
Originally, this model was fit with the primary goal of predicting incidence rather than describing relationships between variables. This of course serves the purpose of doing “predictive model checking”. In particular, we wanted to obtain predictive distributions of incidence. This turned out to be not a straight-forward process.
To obtain predictive distribution samples, we use a two-stage procedure to incorporate the uncertainty from our model parameter estimates and from the negative binomial distribution. We first draw 100 sample parameter sets from a multivariate normal distribution with mean equal to the point estimates of the parameters (\(\boldsymbol{\theta}, \mu, \sigma^2\)) from Equations ()-() and covariance equal to the matrix of standard errors. Each of these sampled parameter sets yields a corresponding \(\widehat{\lambda}_{i,t}\). We then draw 100 samples from the negative binomial distribution given in Equation () for each \(\widehat{\lambda}_{i,t}\) with the fixed estimate of \(r\) to obtain a sample of size 10,000 from the predictive distribution for \(Y_{i,t}\). We calculate the point estimate for each province-year, \(\widehat Y_{i,t}\), as the median of these samples from the predictive distribution. The lower and upper limits of the 80% prediction intervals were defined by taking the 10\(^{th}\) and 90\(^{th}\) percentiles of these samples from the predictive distribution.
Above, the predicted incidence is in black (with 80% prediction interval) and the red dots are the observed data. In the paper itself, there are similar graphics with comparisons to a simple model that just predicts the average incidence in a given year.
The 80% coverage for the two test seasons can be calculated as:
test_dat$coverage <- with(test_dat, (obs >= ci_80_lb) & (obs <= ci_80_ub))
sum(test_dat$coverage)/nrow(test_dat)[1] 0.7236842