Ottar N. Bjørnstad, Bärbel F. Finkenstädt, and Bryan T. Grenfell. 2002. Dynamics of measles epidemics: estimating scaling of transmission rates using a time series SIR model. Ecological Monographs 72:169-184.

The epidemic process in Eqs. 1 to 4 is defined by the observed variables I_t  (subject to measurement error) and B_t, the partially observed variable S_t, and an unobserved variable . A complete treatment will thus involve estimating the parameters in a hierarchical nonlinear model that involves three levels of stochasticity: the migration process nested within the birth-death process nested within the measurement process. We are currently working on a direct likelihood-based framework (using MCMC) to estimate all parameters simultaneously (Finkenstädt et al., unpublished manuscript). The reported result, however, arise from a more simplistic and more approximate procedure were we estimate  marginally on all the other parameters (using profile likelihood); The and the seasonal parameters -- a total of 27 parameters -- are estimated conditionally on  and marginally on m using standard linear likelihoods. Finally, we estimate m using a somewhat ad hoc nonlinear minimization. We sketch out each step below:

Step 1. On a log-scale Eq. 1 is:

.

(B.1)

The first term is linear after transformation (and represents an intercept). The last term represents an offset that depends on one unknown parameter, . If the model were to have been linear in all other parameters, it would be easy to estimate this parameter using profile likelihood marginal on the other parameters; i.e., finding the MLE for each of the linear parameters for each candidate of , then trace out the likelihood across the grid of candidates. The maximum of this marginal likelihood profile is the optimal estimate of . The word marginal refers to the fact that all the linear parameters are allowed to vary freely to maximize the likelihood. (Since the marginal likelihood follows the laws of the likelihood [Murphy and Van der Vaart 2000], the minimum -2 log Likelihood is asymptotically -distributed, so it is easy to erect a confidence envelope [e.g., Aitkin et al. 1989].) However, the second term complicates this, since it is nonlinear in the random coefficient, . To deal with this we write , where m is the mean immigration intensity as defined in the main text, and  is the stochastic fluctuation around this mean. Note that . We can thus moment-expand the second term of Eq. A.1. Whenever the migration rate is relatively modest, a first order expansion will provide a fair approximation according to:

,

(B.2)

where the stochastic component of  has zero mean, but the variance clearly depends on the variance of  and the number of previous cases. Substituting B.2 into Eq. 1 leads to Eq. 5. The error around this expectation has zero expectation but a compound variance partially determined by the epidemic birth-death process and partly due to the variability caused by the stochastic migration. Since Eq. 5 is linear in the covariates log I_t, and 1/I_t, and the error term has zero-expectation, we are permitted to use this estimating equation to provide the marginal profile likelihood for . (We assume the variance on [on the log-scale] to be constant. This is not exact. However, we have played with various heterogeneous error models and the parameter estimates appear robust).

Step 2. Given , the parameters , , and apparently m, are easily identifiable from Approximation 5, since this provides a conditionally (on ) linear model. This conjecture holds -- also in the presence of measurement error on I_t -- except for the parameter m; Let us assume that we observe I_t with (log-normal) error, so that we observe , where U_t is a zero-mean normal variate. If we estimate the parameters of Approximation 5 but substituting  for I_t throughout, we find that:

.

(B.3)

If we take the expectation it is easy to see that all parameters remain unbiased in the presence of measurement error on I_t, except for the fudge-factor , which will be inflated by a factor of . Since this expectation is always greater than unity (but depends on the measurement error variance), the estimated migration rate will be biased upwards. We thus need some alternative method to for estimating m.

Step 3. The previous steps provide consistent estimates of all parameters except for m. We have attempted a variety of approaches to obtaining unbiased estimates of this last parameter. Our current understanding is that this is not possible within the regression setting without incorporating all three stochastic levels (including the observational process) simultaneously. We are currently working on developing such an explicit framework. In the meanwhile, we estimate m through stochastic simulation and by noting that the distribution of epidemic sizes is critically dependent on the migration rate. Frequent migration will necessarily result in small but regular epidemics every year. Very weak coupling, in contrast, results in infrequent but massive outbreaks. We use this property to obtain an estimate of m by running a large number of stochastic simulations of the full model (1-4). For each of these we calculate the distributional distance between the (marginal distribution of the) simulated and the observed time series. We estimate these distributions nonparametrically using kernel density estimation with bandwidth optimized using cross-validation (Silverman 1986) . We select the value for m that minimizes the Kullback-Leibler (K-L) distance between the observed and simulated (marginal) distribution. (We chose this distance measure because of its general information theoretical relationship to the likelihood, but the choice of distance measure is arbitrary). We have carried out a sequence of numerical experiments to confirm that the correct value can be retrieved using this procedure (Finkenstädt et al., unpublished manuscript). However, we have not been able to develop any tight links between the K-L profile and any type of profile likelihood. This step is therefore rather tentative. We erect approximate confidence envelopes for the estimates by spline smoothing the scatterplot of K-L distances against candidate m-values (across a grid from 0.01 to 20 per biweek). The CI’s were erected as the inverse problem involving the spline smooth and the mean-squared error around the smooth (much akin to the profile likelihood, but using the empirical variance rather than some asymptotic -distribution). Work on this is still in progress.

Literature cited:

Aitkin, M., D. Anderson, B. Francis, and J. Hinde. 1989. Statiscal modelling in GLIM. Clarendon Press, Oxford, UK.

Murphy, S. A., and A. W. Van der Vaart. 2000. On profile likelihood. Journal of the American Statistical Association 95:449-465.

Silverman, B.W. 1986. Density estimation for statistics and data analysis. Chapman and Hall, London, UK.