A comparison of several theories of epidemic dynamics using AIC

Craig W. McCarty¹, Kenneth P. Burnham², and Michael W. Miller³

¹Graduate Degree Program in Ecology,

Colorado State University, Fort Collins, CO 80523, U. S. A.

²Colorado Cooperative Fish and Wildlife Research Unit,

Colorado State University, Fort Collins, CO 80523, U. S. A.

³Colorado Division of Wildlife, Wildlife Research Center,

317 West Prospect Road, Fort Collins, CO 80526, U. S. A.

We define and compare 4 classes of epidemic models (41 total) using 2 data sets (measles, common cold) from the literature and model selection based on Akaike's information criteria. The 4 classes include multinomial, structured multinomial, nested logistic, and biological-mechanistic. Although model selection results are not entirely consistent, one model emerges as superior: a mechanistic model addressing heterogeneity but denying the classical assumption that risk of infection increases in proportion to the number of infectious individuals in the population (i.e., the mass action assumption). Our findings underscore the need for replicated data and experimental evidence to support selection of any particular epidemic model to explain or predict host/parasite dynamics.

Key words: Greenwood, Reed-Frost, mass action assumption, binomial chain, epidemic modeling, model selection, statistical inference, AIC, measles, common cold

    Models of specific epidemics are abundant in the literature. Several general theories of epidemic dynamics exist and have allowed authors a variety of approaches to specific problems. For obvious reasons, experimental challenges to these various epidemic theories are rare and most analyses are based on opportunistic data (D'Amico et al., 1996). Given the well known difficulties of such data sets, we are surprised that (1) so few authors have compared alternative models from within a statistical framework, (2) of these, even fewer have applied the modern paradigm of information theory based model fitting and selection, and (3) that no authors have entertained more than 4 candidate models. Here we define and compare 4 classes of epidemic models (41 total) using 2 classic observational data sets from the literature and model selection based on Akaike's information criteria corrected for small sample bias (AICc; Burnham and Anderson, 1998).

    We divide models of epidemic dynamics into two general categories: "statistical models" which attempt only to describe the pattern of the data, and "mechanistic models" which attempt to mimic the processes assumed to have generated the data. The pivotal issue of mechanistic models of epidemic dynamics is how to mathematically describe the transition from the susceptible to the infected state throughout the course of the epidemic. Authors have rationalized a variety of approaches, modeling this transition as dependent on various combinations of constants, random variates, and functions of some combination of the numbers of infectious, susceptible, and resistant individuals in the population. Although often neglected, many other phenomena are potential sources of significant effects. First, does infectiousness or susceptibility vary by individual? Second, does the mechanism of transmission vary by situation or pathogen? Third, what is the consequence of the population being open to births/deaths and emigration/immigration during the epidemic? Fourth, does infectiousness and/or susceptibility change through time independently of the population's composition in terms of the number of infectious and/or susceptible individuals present during the epidemic. Finally, does the spatial arrangement and movements of members of the population affect transmission?

2. Data

    A literature search revealed two data sets useful in addressing these questions (Table 1). Both are frequency tabulations of epidemic "binomial chains" in which the progress of successive viral infections was tracked over time within human family households (e.g., the chain 1-2-0 implies 1 initial infection which caused 2 new infections which caused no new infections). The first data set, "common cold," is a tabulation of binomial chains from 664 households of size 5 with 1 initial infection in each household (Schenzle, 1982; Becker, 1998). The second data set, "measles," is a tabulation of binomial chains from 434 households of sizes 3 and 4 with 1 initial infection in each household (Bailey, 1975).

    Together these 2 data sets display variation in several of the components of the epidemic problem. Each household of each data set represents a replicate of a specific epidemic. Both data sets exhibit variation in the number of infectious and susceptible individuals during the course of the epidemic. The measles data, while not open to births/deaths or emigration/immigration, shows variation in household size, thereby simulating this effect via comparison between households of different sizes within the same data set. In both sets, the spatial component to the epidemic process is largely constant because, when compared to entire communities, houses are approximately the same size and are likely to share similar daily activity dynamics. Further, households are most likely to meet the "homogeneity of mixing" assumption of most mechanistic models of epidemic dynamics and necessary for many of the models compared here. Unfortunately, additional within-household covariate information was not collected by the investigators (e.g., age and gender composition of each household).

3. Models

3.1 Rationale

    For convenience we divide our set of 41 candidate models into 4 classes (multinomial, structured multinomial, nested logistic, biological-mechanistic) based on our interpretation of each individual model as "mechanistic" or "statistical" and whether or not individual models are nested within the global model of each class. For example, model class 2 contains 14 models nested under one global model and is considered to be more mechanistic than model class 1, but less so than model class 3.

3.2 Model class 1, multinomial (1 member)

    The multinomial model, MULT, recognizes no generating process beyond the simple observation of different categories of chains, each with different cell probabilities. For example, under this model the likelihood of the measles data is

All subsequent models add structure to these underlying multinomial models by recognizing that each cell represents a specific binomial chain (e.g., chain 1-1-0 corresponds to p₂ in the above multinomial) and equating each cell probability with a function of the population status whose form is based on mechanistic and/or statistical assumptions about how that specific chain came about.

3.3 Model class 2, structured multinomial (15 members)

    Global model NTSI suggests that the probability of transition from the susceptible to the infected state is conditional on the population size (N), the sequence step (T), the number of available susceptible individuals (S), and the number of infectious individuals (I). The number recovered (R), is not included since R=N-S-I. Thus for example, the probability of the chain 1-1-0 from households of size 5 with 1 initial infection under model NTSI is 4p_{N=5,T=1,S=4,I=1}q³_{N=5,T=1,S=4,I=1}q³_{N=5,T=2,S=3,I=1}. Model class 2 contains 14 nested sub-models (N, T, S, I, NT, NS, NI, TI, TS, SI, NTS, NTI, NSI, TSI), some of which are redundant when applied to these 2 data sets because N = 3 or 4 is small.

3.4 Model class 3, nested logistic (16 members)

    Nested logistic models, global model logitNTSI, considers the probability of transition from the susceptible to the infected state as a logit transform of a linear combination of N,T,S, and I. Included in this class are the obvious 14 sub-models and an additional model representing our interpretation of the popular "mass action assumption." Unfortunately, the "mass action assumption," has no unique interpretation as a stochastic model or in discrete time (Bailey, 1975; Heesterbeek and Roberts, 1995): "The I infectives make on average I potentially infective contacts per unit of time. Only a fraction S/N of those contacts is with a susceptible. So, IS/N is the number of new cases arising per unit of time." We chose to represent the "mass action assumption" as a logit transform of I/N giving model logitI/N.

3.5 Model class 4, biological-mechanistic (10 members)

    For heuristic purposes, all members of this class are considered as modifications of the McCarty-Miller model, MM (McCarty and Miller, 1998; McCarty, 1999):

This model recognizes two possible biological mechanisms of disease transmission: (1) a random apportionment of Ib (the product of I infectious individuals and b infectious contacts per infectious individual) total infectious contacts among the N individuals of the entire population at each time step, and (2) a possible constant level of background risk (1- a).

    Model MM simplifies to the Reed-Frost model, RF, if N is held constant and a=1 (i.e., no background risk) and the Greenwood model, G (Bailey, 1975), if b=0 and N is held constant. Consequently, the Greenwood model has a mechanistic interpretation as a single point source of infection, making the probability of transition from the susceptible to the infected state a constant regardless of population composition. Expanding model G by replacing a with a beta random variate gives model BG which accounts for variation in the degree of risk to infection from a point source specific to a population (Bailey, 1975); allowing separate beta random variates for each population size gives model BGN. Model MM can be further modified to give the beta-Reed-Frost model, BRF (Bailey, 1975; Becker, 1989), if a=1, and the term (1 - 1/N)^b is replaced by a beta random variate; conditioning on N gives Model BRFN. Finally, included in this class is the propagating point source model, PP, which considers each new infectious individual as a new point source of risk to infection:

4. Model fitting and selection

    Model selection based on AICc requires 4 steps (Burnham and Anderson, 1998): selection of candidate models for each data set, maximizing the resulting likelihoods, calculation of an AICc value for each model-data set combination, and comparison of the resulting AICc values with selection leaning towards models with the lowest AICc values. Models differing in AICc value from the best model by more than 4-7 units are generally poor candidates. Alternatively, Akaike weights provide a useful scaled metric (0-1) of model selection uncertainty.

    Because previous authors have questioned the binomial chain interpretation of these data sets, suggesting that chains may have been incorrectly classified, we considered underlying multinomial models of each data set based on both the chain interpretation and on the final size of each epidemic (e.g., chains 1-1-1 and 1-2 both give a final size of 3; Longini and Koopman, 1982; Becker, 1989). Standard approaches give simple closed-form solutions for maximum likelihood estimates of the parameters for all models in classes 1 and 2 (Bailey, 1975; Becker, 1989). We used numerical optimization in cases where no simple estimator exists.

5. Results and discussion

5.1 Complications

    Although model BG is the best overall candidate among our analyses, model selection results are not entirely consistent between data sets and among analyses (Tables 2, 3). Curiously, in opposition to traditional tenets of epidemic dynamics, neither data set suggests that the number of infectious or susceptible individuals, population size, or time strongly affected epidemic dynamics in these examples. While our analyses raise several questions about the validity and utility of traditional mechanistic theories of epidemic dynamics, three issues must be considered as we attempt to address these concerns.

    First, were the data collected precisely? Misclassification of chains could have occurred in numerous ways. For example, some instances of the chain 1-2-1 could have been incorrectly classified as 1-1-1-1: analyses based on final size attempt to address this issue (Bailey, 1975). Also, we have no guarantee that new infections were not introduced into some households from outside sources during the period of data collection (Schenzle, 1982; Longini et al., 1988): models incorporating heterogeneity attempt to address this issue.

    Second, did heterogenous factors affect outcomes? Perhaps, households of size 4 contained more young children than households of size 3, and children are more or less susceptible or infectious than adolescents and adults. Also, were resistant individuals present in any households prior to the start of the study? This seems more likely in epidemics of measles than the common cold because, unlike the common cold, prior infections of measles confer long-term immunity.

    Third, were the true generating processes of these data complicated enough to mask the effects of an underlying general mechanism like the "mass action assumption?" Since AIC selects models which best balance the tradeoff between bias and variance (Burnham and Anderson, 1998) it is virtually certain that different models would be selected given larger sample sizes and "cleaner" data. However, if the popular mechanistic models of epidemic dynamics are useful concepts, such mechanisms should at least strongly effect outcomes (McCarty, 1999).

5.2 Evidence for I, S, T, and N effects

    Four observations illustrate a lack of conformity between current theories of epidemic dynamics and these data sets. Neither data set suggests a strong I effect; that risk to infection increases or decreases when the number infectious individuals in the population exceeds 1. All model-data set combinations incorporating an I effect gave individual Akaike weights of <0.14 and commutative weights for this effect among models within each analysis were <0.2 in all cases. Further, while we recognize that the measles data set shows variation in I for only households of size 4, separate analyses by household size resulted in no remarkable differences in the model selection results; model BG is still strongly favored even when the data are partitioned by household size and no evidence for an I effect emerges. However, it is important to note that the highly infectious nature of these pathogens and the relatively large proportion of the population initially infected may have overwhelmed the I effect. Second, the measles data only hint at a possible N effect as model BGN ranked second in one analysis (likelihood based on chains using household of sizes 3 and 4) but with a low relative weight, making it 5 times less likely than the selected model. The 30% increase in population size is apparently not large enough to detect an N effect with this sample size (i.e., N_{size 3} = 334, N_{size 4} = 100) if it is present in these data. Third, while a T effect seems biologically reasonable, only the common cold data offer any support for this notion with the relatively high total weights of models logitT, logitTI, and logitTS. It is difficult to argue that even a weak T effect could have been masked by heterogeneity in the common cold data, given the relatively low rank of model BG and the extremely good fit of virtually all candidate models. Fourth, reasonable evidence is suggested for a monotonic S effect in the common cold data (i.e., a cumulative weight of 0.55 among models incorporating this effect when likelihoods were based on chains) but clearly not in the measles data. We can give no biological explanation for this outcome.

5.3 Evidence for simple generating mechanisms

    Selection of the beta-Greenwood model, BG, has a pleasing mechanistic interpretation as constant background risk varying in intensity among households. Explanatory hypotheses are plentiful. Perhaps infections were spread to susceptibilities through contact with common infected objects within each home such as food, towels, sinks, or eating/drinking utensils; or perhaps the infectious agent was uniformly distributed in the air throughout each home. These conditions could produce constant risk; the addition of new infectious individuals may have little or no further effect other than replenishing and maintaining risk through time.

5.4 Conclusion

    "A discipline is a science only if the theories it entertains are universal and falsifiable" (Popper, 1972). Thus, although we recognize the difficulty of falsifying theories expressed as probability statements (Dolby, 1982), we are concerned with the proliferation of models of specific epidemics in the absence of a well-supported underlying theory of epidemic dynamics. Even a cursory literature search reveals hundreds of such models published since 1930 (McCarty, 1999). Further, although it may be unrealistic to expect that a single unifying model of epidemic dynamics can adequately represent the myriad of host/parasite relationships encountered in nature, we are troubled by the lack of significant supporting data sets, replication, and experimental evidence supporting any particular epidemic model. While the data analyzed here are admittedly imperfect, the results of our analyses raise questions about the propriety of blindly invoking popular epidemic models in attempts to explain or predict specific host/parasite dynamics.

ACKNOWLEDGMENTS

Support from the Colorado Division of Wildlife, W-153-R is gratefully acknowledged. Richard Bartman, Sara Oyler-McCance and Mary Conner provided useful reviews of earlier drafts.

REFERENCES

Bailey, N. T. J. (1975). The mathematical theory of infectious diseases and its applications, second edition. London: Charles Griffin and Company.

Becker, N. (1989). Analysis of infectious disease data. New York: Chapman and Hall.

Burnham, K. P. and Anderson, D. A. (1998). Model selection and inference: A practical information theoretic approach. New York: Springer-Verlag.

D'Amico, V., Elkinton, J. S., Dwyer, G. and Burand, J. P. (1996). Virus transmission in gypsy moths is not a simple mass action process. Ecology 77(1), 201-206.

Dolby, G. R. (1982). The role of statistics in the methodology of the life sciences. Biometerics 38, 1069-1083.

Heesterbeek, J. A. P. and Roberts, M. G. (1995). Mathematical models for microparasites of wildlife. In Ecology of Infectious Disease in Natural Populations, B. T. Grenfell and A. P. Dobson (eds), 90-122. Cambridge: Cambridge University Press.

Longini, I. M. Jr. and Koopman, J. S. (1982). Household and community transmission parameters from final distributions of infections in households. Biometrics 38, 115-126.

Longini, I. M. Jr., Koopman, J. S., Monto, A. S. and Fox, J. P. (1982). Estimating household and community transmission parameters for influenza. American Journal of Epidemiology 115(1), 736-751.

McCarty, C. W. (1999). Contributions to the mathematical theory of epidemic dynamics. PhD dissertation. Colorado State University.

McCarty, C. W. and Miller, M. W. (1998). A versatile model of disease transmission applied to forecasting bovine tuberculosis dynamics in white-tailed deer. Journal of Wildlife Diseases 34, 722-730.

Popper, K. R. (1972). The Logic of Scientific Discovery. London: Huchinson.

Rampey, A. H. Jr., Longini, I. M. Jr., Haber, M. and Monto, A. S. (1992). A discrete-time model for the statistical analysis of infectious disease incidence data. Biometrics 48, 117-128.

Schenzle, D. (1982). Problems in drawing epidemiological inferences by fitting epidemic chain models to lumped data. Biometrics 38, 843-47.