As general approaches, all of the methods mentioned earlier face severe limitations. Some must modify the model, discretizing time to the scale of observations (Finkenstädt & Grenfell 2000; Grenfell et al. 2001, 2002; Bjørnstad et al. 2002; Finkenstädt et al. 2002; Clark & Bjørnstad 2004; Koelle & Pascual 2004). They rely upon fortuitous congruence between some time scale of the process under study (e.g. infection generation time) and that of the sampling interval. These methods are computationally cheap, but can lead to severely biased conclusions (Glass et al. 2003). Others rely on approximations that are valid only in restricted circumstances. For example, Cauchemez & Ferguson (2008) adopt a diffusion approximation to enable successful implementation of a Markov chain Monte Carlo approach; this approximation breaks down in small populations. On the other hand, straightforward Bayesian approaches have been demonstrated for small populations with purely demographic stochasticity (Gibson & Renshaw 1998, 2001); however, these methods break down for large populations (Roberts & Stramer 2001; Beskos et al. 2006) or when environmental stochasticity is non-negligible (Bretó et al. 2009). Finally, exact methods based on numerical solution of the so-called master equations (Keeling & Ross 2008) are available, but limited to low-dimensional models and small populations. In summary, each of these approaches makes demands upon the form of the model and/or the nature of the data and thus necessitates scientifically irrelevant or even inappropriate assumptions that can interfere with the hypotheses of interest.

In this paper, we use a relatively new approach with none of the aforementioned limitations (Ionides et al. 2006, 2007, 2009; King et al. 2008; Bretó et al. 2009). In particular, the method works on continuous time models, in large or small populations, with demographic or environmental stochasticity or both. It has the important advantage of operation based only on simulation from the dynamic model, without a need for explicit expressions for transition probabilities. The latter property has been called plug-and-play (Bretó et al. 2009). Plug-and-play inference methodology provides a powerful tool to enable carrying out data analysis via flexible, scientifically motivated classes of models. In particular, plug-and-play is extremely useful when one wishes to entertain multiple working hypotheses translated into multiple mechanistic models. Two properties analogous to plug-and-play have been considered in the contexts of optimization and complex system theory. These analogies suggest that plug-and-play is a subject whose importance extends beyond our present topic of inference for partially observed stochastic dynamic systems. The conceptual commonality is that all these methodologies can be applied to models specified as ‘black-box’ computer codes. This enables flexible model development and therefore avoids confounding methodological issues with scientific hypotheses. It also facilitates both the comparison of alternative explanations and the transfer of methodological expertise between different systems. The price of plug-and-play is primarily in terms of computational effort: when properties of a model are available in closed form, computationally more efficient methods typically exist. As computational capabilities continue to grow, the ease with which new scientific questions can be asked and answered via plug-and-play methodology will make such techniques increasingly attractive.

We used the implementation of iterated filtering contained in the software package Pomp (King et al. 2009) written for the R statistical computing environment (R Development Core Team 2006). This implementation follows the algorithm described by King et al. (2008, supplementary information therein). Heuristics and diagnostics, together with a discussion of computational issues, can be found in Ionides et al. (2006); a theoretical treatment is given by Ionides et al. (2009).

For complex dynamical systems, not all questions that one might wish to ask will be clearly answered by the available data. Typically, a question posed of the dynamic system is formalized as a hypothesis about a combination of parameters in a model for the system. If the likelihood function does not discriminate well between different values of some parameter combination, this parameter combination is said to be poorly identified. In such a case, one might proceed by making additional assumptions: adding constraints typically improves parameter identifiability at the cost of potential bias in the estimates of other parameters. In §3, we focus on questions the data are able to answer without further assumptions. A basic tool for this is the profile likelihood function, in which one parameter is varied over a range of values, while the likelihood is maximized over all remaining parameters (Barndorff-Nielsen & Cox 1994; Hilborn & Mangel 1997). In practice, we work with a smoothed estimate of the likelihood surface calculated via Monte Carlo (Ionides 2005).

The data we analyse consist of weekly measles case reports on 20 representative communities with population sizes ranging from 2000 to 3 000 000. The towns and cities were selected from 954 urban locations for which case reports have been compiled (Grenfell & Bolker 1998); specifically, we chose the largest 10 cities by population and sampled 10 more towns at random. Figure 1 shows the case reports and annual birth rates for two of these: London and Hastings. These data and their limitations have been discussed previously (Fine & Clarkson 1982; Bjørnstad et al. 2002). In particular, although measles became a notifiable disease in England and Wales in 1940, we follow Fine & Clarkson (1982) by starting our analysis in 1950 to avoid the disruptions caused by the second world war and the transition period following the introduction of the National Health Service in 1948. From 1950 onwards, it is believed that reporting was carried out fairly consistently, with between half and two-thirds of all cases being counted.

A basic susceptible–exposed–infectious–recovered (SEIR) model for measles divides the population into those individuals susceptible to infection, exposed (i.e. infected but not yet infectious), infectious and removed (i.e. quarantined or recovered and subsequently immune). The model is specified by describing the rates at which individuals move between compartments (Jacquez 1996; Matis & Kiffe 2000); figure 2 shows the transitions for the SEIR model of measles. A diagram such as figure 2 can be unambiguously interpreted as a deterministic system, a system with demographic stochasticity or a system with both demographic and extra-demographic stochasticity (Bretó et al. 2009). We give the equations for each of these interpretations in appendix A. Here, we will employ the stochastic interpretation, with both demographic and extra-demographic stochasticity, described in box 1 (in appendix A). The transition rates in a model such as figure 2 may, in general, depend on the state of the system and/or covariates. We specify the force of infection as

Here β(t) is the transmission rate; ι the mean number of infectives visiting the population at any given time; α a mixing parameter, with α = 1 corresponding to homogeneous mixing (Liu et al. 1986; Bjørnstad et al. 2002), and N(t) the population size, treated as known via interpolation from census data. Since transmission rates are closely linked to contact rates among children, which are higher during school terms (Schenzle 1984; Keeling & Grenfell 2002; Bauch & Earn 2003; Conlan & Grenfell 2007), we assume that β(t) reflects the pattern of British school terms and holidays. Specifically, we take where p = 0.759 is the proportion of the year taken up by the school term, the mean transmission rate and a the relative effect of school holidays on transmission. For ease of interpretation, it is sometimes convenient to reparameterize in terms of the annual average basic reproductive ratio, R₀, defined as the expected number of secondary infections engendered by an infective introduced into a fully susceptible population. When the duration of infection is much shorter than life expectancy, R₀ ≈ /μ_IR; here we employ a modification of this formula, given in appendix A.

A novel feature of our model is that we add white noise with intensity σ_SE to the transmission process. Some empirical support for the choice of white (i.e. temporally uncorrelated) noise is given by Lande et al. (2003). This variability in the rates allows the possibility of over-dispersion (McCullagh & Nelder 1989; Bretó et al. 2009). It is discussed at more length, together with a full specification of all the model equations, in appendix A.

The other transition rates are specified as follows: μ_EI is the rate at which exposed individuals become infectious, thus μ_EI⁻¹ is the mean latent period; μ_IR is the recovery rate; μ_SD = μ_ED = μ_ID = μ_RD = m denotes a constant per-capita death rate, ignoring the negligible effect of disease-induced mortality; μ_BS(t) is the per-capita rate of recruitment of susceptibles, depending on the birth rate b(t), which is assumed to be known via interpolation from birth records. Births enter the susceptible class with a delay τ corresponding to the age at which children enter the high-risk school-age population. We also consider a cohort-entry effect that reflects the fact that a large cohort of first-year students—the majority of them serologically naive—enters the schools each autumn (Schenzle 1984; He 2006). Specifically, a fraction c of recruits into the susceptible class enter on the school admission day, and the remaining fraction (1–c) enter the susceptible class continuously. The cohort effect provides a parsimonious and mechanistically plausible alternative to previous suggestions that the transmission rate β(t) may increase following the start of the academic year (Fine & Clarkson 1982; Grenfell et al. 2002). The reporting process is taken to be on over-dispersed binomial, with reporting rate ρ and overdispersion parameter ψ (appendix A, equation (A 6)). Quantities appearing in the model are summarized in table 1.

In table 2, the homogeneity parameter α is consistently found to be close to 1. At first glance, this may suggest that there is little evidence for inhomogeneity of transmission in the data. Profile likelihood plots over α for London and Hastings are presented in figure 3. We see that the statistical uncertainty in the estimates for these two places is comparable to the geographic variability in the point estimates in table 2 (mean ± s.d. of 0.99 ± 0.03). This finding confirms the suggestion by Glass et al. (2003) that such exponents may be artefacts of the time discretization employed in earlier approaches (Bjørnstad et al. 2002; Grenfell et al. 2002). We can conclude that the phenomenon of inhomogeneous mixing—inasmuch as it plays a dynamic role—is best captured in some other way.

The imported infection parameter ι is well identified in these data. There is a strong log-linear relationship between ι and N, consistent with the earlier findings of Bjørnstad et al. (2002) and with the results obtained by Xia et al. (2004), who analysed a much more extensive dataset using a gravity model. The evidence for such a relationship is stronger here than in a previous empirical investigation (Finkenstädt et al. 2002, fig. 6).

The extra-demographic stochasticity parameter σ_SE is consistently found to be around 0.06 (0.063 ± 0.015). The profile likelihood plots in figure 4 also indicate that this parameter is fairly well estimable, and the data strongly argue against the choice σ_SE = 0. Figure 4 suggests that a role for extra-demographic stochasticity is not only strongly indicated for large cities such as London but also clearly evidenced in smaller populations for which one might expect demographic stochasticity to dominate. One interpretation of this is that models with purely demographic stochasticity (i.e. all discrete-population continuous-time dynamic models that have been previously fitted to disease data) fail to capture an extra-demographic source of variability that has an important dynamic role. A priori exclusion of the possibility of such additional variability can lead to spurious confidence in the conclusions of the analysis (McCullagh & Nelder 1989) and to biased parameter estimates. Figure 4 shows that, in the absence of extra-demographic stochasticity, estimates of the latent and infectious periods are substantially lower when one allows for such stochasticity. This appears to arise because short transition times lead to fewer individuals in compartments E and I on average and therefore to a larger dynamic role for demographic stochasticity. In other words, the model compensates for the lack of extra-demographic stochasticity by attempting to increase the intensity of demographic stochasticity.

Each parameter of the model has an interpretation in terms of the basic biology of transmission and infection at the level of individual hosts. To the extent that the model is a faithful description of the population-level processes, one expects that parameter estimates should agree closely with estimates based on individual-level observations. When discrepancies occur, they are evidence for model misspecification or oversimplification. Strikingly, table 2 shows a strong log-linear relationship between population size and infectious period and, to a lesser extent, latent period. Since there is no evidence that the natural history of measles varies with population size, this relationship is indicative of such oversimplification. Figure 1 suggests that this simple model fails to account adequately for heterogeneity in transmission. Specifically, the model calls for short infectious periods and high transmissibility β to reproduce the explosive epidemics and interepidemic fadeouts typical of small communities. The data from large populations, in contrast, represent the aggregate of many local epidemics, in each of the many schools and local communities contained in a large conurbation. In these aggregated data, fadeouts are rare and the epidemic curve is generally shallower. The oversimple model can reproduce these features only by extending the infectious period and reducing transmissibility. Moreover, while the infectious period estimated in the small populations (2–4 days) is consistent with evidence from household studies (Hope Simpson 1948, 1952; Bailey 1956), the longer infectious periods estimated in the large populations are not. Previous studies have noted that deterministic versions of the mass-action model analysed here do a remarkably good job of qualitatively reproducing the dynamics of measles in large communities (Earn et al. 2000). Our analysis shows, however, that when both continuous time dynamics and stochasticity are taken into account, data from small populations afford a less-distorted view of the disease dynamics, at least when viewed through the lens of relatively simple models.

Clinical and household studies of infectious diseases yield information regarding the transmission and progression of infection on the scale of individuals or families (Hope Simpson 1948, 1952; Bailey 1956). Such studies are complemented by data on disease prevalence or incidence on larger spatial scales. The recent development of likelihood-based approaches capable of dealing with the full spectrum of unavoidable complexities (Cauchemez & Ferguson 2008; Cauchemez et al. 2008; King et al. 2008)—unobserved variables, measurement error, process noise, nonlinearity, non-stationarity and covariates—means that we can now view a disease from the population point of view with something like the same clarity that we have for many years been able to view individual infections. Specifically, although we have long understood how to describe the population dynamics of infectious diseases using mathematical models, we have only recently gained the ability to fit these models to data using statistically sound methods. When biological quantities estimated in this way agree with those estimated from clinical or household studies, it can be interpreted as confirmation of the assumptions embodied in the model. Disagreement between the small- and large-scale points of view, however, raises interesting scientific questions. In the present study, we find broad agreement with previous studies concerning many of the model's parameters. On the other hand, we find significant departures with respect to three important quantities—R₀, the infectious period and the intensity of extra-demographic stochasticity.

The quantity R₀ is central in epidemiological theory because it has interpretations in terms of so many quantities of interest, including mean age of first infection, mean susceptible fraction, exponential-phase epidemic growth rate and vaccination coverage required for eradication. It is important to realize that the conjunction of these interpretations occurs only in the context of very simple models. Simple models necessarily lack flexibility. In reality, these interpretations diverge owing to heterogeneities in age, spatial location, host genetics, etc. It is therefore unrealistic to expect that estimates of R₀ (or any other single quantity) derived from fitting a simple model to one sort of data should agree with estimates derived from other data sources. Rather, the key question should be: which biological interpretations are relevant in the context of the data used to inform the model? From a statistical point of view, this corresponds to the question: to what features of the data are the parameter estimates sensitive? Our estimate of the basic reproductive ratio, R₀, like those of earlier studies focusing on aggregated case-count time series, is somewhat greater than estimates based on serological surveys and age-stratified incidence data. It is possible that this reflects the sensitivity of the time-series-derived estimates of R₀ to peak incidence. Specifically, peak incidence strongly influences the estimates of both R₀ and reporting rate ρ, but since ρ is well identified by other features of the data (namely, the long-run cumulative incidence), we suspect that the model requires a high R₀ to match the observed peak epidemic case counts. Alternatively, the high R₀ estimates may reflect the sensitivity of estimated R₀ to early phase epidemic growth. If this is in fact the case, our estimates of R₀ may more accurately reflect contact rates among the core group of school-age children than they do those of the population at large. In this case, the model is effectively extrapolating these rates to the adult population, about which these data have little direct information since so few adult cases occur. In contrast, from this point of view, the interpretation of R₀ in terms of its definition as mean number of secondary infections in a fully susceptible population is hopelessly extrapolated. Likewise, the interpretation of R₀ in terms of mean age of first infection is unjustified, since it necessarily depends on the age structure of transmission, which is not part of the model. Perhaps the most direct information pertaining to mean age at first infection in these data comes from the lag between changes in birth rate and their subsequent effects on incidence. This lag is explicitly captured in our model by the delay, τ, between birth and recruitment into the susceptible pool.

Since the comparatively high estimate of R₀ does not appear to be a mere artefact of time discretization, the question remains as to why the population-level data suggest a more communicable disease than do the individual-level data. Existing estimates of R₀ are sensitive to assumptions about the age structure of transmission that have not yet been fully resolved (Wallinga et al. 2001). Perhaps transmission in schools is relatively more effectual than within households. Alternatively, it may be that this discrepancy points to a need for a better description of infectious- and/or latent-period distribution, of the disease's age structure or both. In our results, the departure of the estimated latent and infectious periods from plausible values obtained from household studies grows with city size, and so the most likely explanation for this correlation is the spatial heterogeneity of transmission within towns and cities. The latter, we have shown, cannot usefully be captured via the simple device of an exponent in the transmission term. More detailed explorations of disaggregated data and/or models with explicit spatial structure have the potential to shed light on this question. Finally, the strong evidence in favour of extra-demographic stochasticity raises the question of precisely why such stochasticity aids in the explanation of the data. To what extent does this finding indicate the presence of genuine environmental stochasticity? To what extent does it indicate model misspecification? To address these questions, again, analysis based on more detailed models is called for. In the case of measles, there are sufficient data to entertain models featuring spatial inhomogeneity (Grenfell et al. 1995; Xia et al. 2004; Bjørnstad & Grenfell 2008), age structure inhomogeneity (Schenzle 1984; Bolker & Grenfell 1996; Keeling & Grenfell 1997) or both (Bolker & Grenfell 1995). Although it has been demonstrated that such models can do a good job of accounting for gross features of the data, there has been less emphasis on requiring these models to account for all features of the data.

Progress on inference methodology for parameter estimation from measles time series has emphasized SEIR-type models for homogeneous populations (Ellner et al. 1998; Bjørnstad et al. 2002; Cauchemez & Ferguson 2008; Keeling & Ross 2008). The study of age structure and spatial effects for measles has placed less emphasis on inferring parameters from data (one exception is Xia et al. 2004). This may be partly because of the additional difficulties of inference for such systems and partly because models based on homogeneous mixing do an impressive job of describing key dynamic features (Earn et al. 2000; Grenfell et al. 2002). Regardless of one's view of the importance of paying further attention to population inhomogeneities, there is a natural methodological question: how applicable are the techniques presented here for larger, more complex models? There is a computational price to pay for the convenience of plug-and-play statistical methods. For the results in table 2, one application of the plug-and-play-iterated filtering inference procedure (based on 50 iterations, with a Monte Carlo size of 10⁴) implemented via the Pomp package in R (R Development Core Team 2006; King et al. 2009) took 5 h to run on a desktop machine. Computational effort scales roughly linearly with the number of parameters plus the number of state variables, so one can see that only a modest amount of additional structure could be included without hitting computational limitations. To increase computational efficiency, however, the iterated filtering algorithm can be implemented in a non-plug-and-play mode, in which the filter is either tailored to the particular model or takes advantage of available analytic properties of the transition densities. Such extensions, which would be required for much larger models, are a topic for future research. One can seek inspiration for future possibilities from numerical climate models, for which filtering operations (the computationally intensive step in the iterated filtering algorithm) have been carried out on systems with 5 × 10⁶ state variables (Anderson & Collins 2007). Filtering in such high-dimensional situations requires the development and evaluation of appropriate approximations for the system under investigation.

In this article, we have carried out the first scientific investigation based on a new framework for continuous-time, discrete-state population dynamics with both demographic and extra-demographic noise (probabilistic and statistical properties of this model class were investigated by Bretó et al. (2009). Extra-demographic stochasticity (interpreted as noise in the rates of a discrete population Markov population model) is equivalent to the possibility of multiple individuals moving simultaneously between compartments (Bretó et al. 2009), and, as such, may be due to social events that affect the behaviour of many individuals (e.g. sporting events) or events that change disease transmissibility, such as variations in temperature and humidity. Variability in the rates gives the model additional flexibility that can also describe model misspecification. Analogously, when carrying out linear regression, it is customary to fit a line to data while understanding that the variation of the data around the line corresponds to unknown and unmodelled deterministic effects as well as to random fluctuations. For linear regression, one typically treats both these sources of uncertainty equally, and certainly all the usual standard errors and test statistics do not discriminate between them. We maintain that the same approach can be applied to dynamic models; in other words, the distinction between model misspecification and process stochasticity should be noted, though it will not usually affect the subsequent analysis.