RLadyBug-An R package for stochastic epidemic models.

RLadyBug is an S4 package for the simulation, visualization and estimation of stochastic epidemic models in R. Maximum likelihood and Bayesian inference can be performed to estimate the parameters in a susceptible-exposed-infectious-recovered (SEIR) model, which is a stochastic model for describing a single outbreak of an infectious disease. The package is thus one step towards statistical software supporting parameter estimation, calculation of confidence intervals and hypothesis testing for transmission models.

Introduction

Understanding the dynamics and spread of infectious diseases is a key component in the design and analysis of defensive strategies. As a consequence, a multitude of epidemiological data on infectious diseases in animal, plant and human communities has been collected to gain insights into the underlying biological and epidemiological processes. The SIR (susceptible-infectious-recovered) model and its variants (S-Exposed-IR, SIS, etc.) are the mathematical tools most commonly used in such analyses.

This paper describes a software program for the statistical analysis of a single outbreak in a small population. Special focus is on the spatial spread between subpopulations arranged on a lattice. In veterinary or plant epidemiology such data arise from so-called transmission experiments, where one or more individuals in a controlled environment are inoculated with the infectious disease pathogen. Subsequently, the course of the epidemic is monitored through visual inspection, clinical testing and other methods. The aim ranges in veterinary epidemiology from quantifying disease transmission (Laevens et al., 1999, Stärk et al., 2000) to determining the effect of a vaccine (Dewulf et al., 2001, Meyns et al., 2004).

Because outbreaks induced by transmission experiments are planned and occur in a controlled environment, the produced outbreak data are especially rich. However, the mathematical setup of course also applies to the analysis of ordinary outbreaks. Interesting is also the application to entirely different areas such as the analysis of computer virus in a network (Wierman and Marchette, 2004), spread of the severe acute respiratory syndrome (SARS) (Donelly and Ghani, 2004) or outbreaks in the wards of a hospital (Grundmann and Hellriegel, 2006).

RLadyBug  is a package implemented in R (R Development Core Team, 2006) providing functionality for the simulation, visualization and estimation in stochastic epidemic models. It enwraps the functionality of the Java program used in Höhle et al. (2005) by S4 classes and adds a volume of methods for the visualization of outbreak data and their estimation results. The aim of this paper is to describe the features of the package in order to make it accessible to statisticians and statistically trained epidemiologists who are looking for software to analyse their infectious disease data. Less attention is thus given to the statistical particulars, which are explained in Höhle et al. (2005).

This paper is organized as follows: Section 2 gives a short introduction to stochastic epidemic models, Section 3 introduces the package and illustrates its use by providing the corresponding R code for analyzing a transmission experiment with classical swine fever virus (CSFV). Section 4 provides a discussion.

Stochastic epidemic models

With focus on the software dimension, only a short introduction to stochastic epidemic models in terms of the SEIR model is given. For a more thorough description see Andersson and Britton (2000).

A closed population P is hosted in k units. Each individual in P can be in one of the states susceptible, exposed, infectious or recovered. A spatial dimension is introduced by assuming that the k units are arranged in a lattice. At the beginning of the outbreak, , the number of susceptibles in each unit is , the number of exposed is and the number of infectious is . At time t, an individual j in unit meets infectious at ratewhere denotes the neighbors of (e.g. in the four compass directions). Furthermore, quantifies the within unit transmission rate, whilst quantifies the transmission rate between neighboring units. If a susceptible meets an infectious it becomes exposed. After being exposed at time an individual j has a gamma-distributed incubation time before becoming infectious, i.e. starting from time , j can infect others. Similarly, the infectious period lasts for after which recovery occurs, hence labels the time of recovery.

Objective of analyzing outbreak data with SEIR models is the estimation of the parameters , which permit the effect of control measures to be studied and to generalize the results to different settings. RLadyBug  is a tool to accomplish this ambition.

The RLadyBug package

Where applicable, RLadyBug  uses the new S4 class system of R implemented in the methods package. This implies a more explicit class system requiring the programmer to define classes, slots and generics explicitly (Chambers, 1998). Fig. 1 shows the object oriented hierarchy of the package: The classes LBExperiment and LBLayout represent the data layer, LBOptions the algorithmic component and LBInference the results. With the help of the rJava package (Urbanek, 2006), this S4 framework is connected to the functionality of an underlying Java program performing the computationally intensive operations.

Fig. 1

The object hierarchy of the RLadyBug S4 classes.

Simulation and visualization

To get an understanding of the SEIR model specified in Section 2, simulation based on the Sellke construction (Andersson and Britton, 2000) is utilized to generate data from the model. Code wise, this is accomplished by creating an object of class LBExperiment using the simulate function. Several visualizations of the experiment data can then be generated by plot, e.g. the type=state time position argument shows and as a function of time for each unit u in the lattice. The below stated code creates a simulated epidemic in a lattice, see Fig. 2 . Initially and . Further views of LBExperiment objects can be generated using type=state time or type=state 1 position as arguments in plot. The former shows and , with e.g. , the latter creates an animation illustrating the course of the epidemic by plotting the state given position for a set of fixed time points. As exemplification, the code below uses a lattice with 15 individuals in each unit (a realistic setup for a pig farm) and creates a 20-picture animation of the epidemic. Fig. 3 shows the eighth picture of this sequence: Each stacked bar shows the current ( days) percentage of susceptible, exposed, infectious and recovered in the unit.

Fig. 2

, , and for the three units in a lattice layout.

Fig. 3

Snapshot of the animation illustrating the outbreak in a lattice.

Test data

Besides the ability to simulate data, the package contains several data sets from human and veterinary epidemiology. Amongst others are the Smallpox Epidemic in Abakaliki, Nigeria, analyzed e.g. in Andersson and Britton (2000) or O’Neill and Becker (2001), and the data from CSFV experiments by Laevens et al. (1999) and Dewulf et al. (2001). The next section describes a statistical analysis of the experiment by Laevens et al. (1999) using RLadyBug. In this experiment the spread of CSFV was investigated in a layout with and slaughter pigs. Every second day all pigs still alive were investigated using a virus isolation test based on blood plasma.

Analysis of an CSFV transmission experiment

For full data, i.e. with known event times for all individuals, RLadyBug  provides likelihood or Bayesian inference for . Details about the underlying equations leading to the log-likelihood and posterior distribution can be found in Höhle et al. (2005).

In practice, however, is unobservable. A typical assumption is thus to assume a fixed and known incubation time c and hence compute as . This was also done in the CSFV experiment by assuming . With this assumption all event times are known and can be illustrated as in Fig. 4 : The individual 12:01 was inoculated at time . By maximization of the log-likelihood the maximum likelihood estimators are readily determined using RLadyBug. Nonetheless, assuming a fixed and known incubation time is not very realistic. Höhle et al. (2005) therefore use a Bayesian framework to handle the missing but gamma-distributed exposure times. Unknown (or censored) waiting times are imputed and updated through a Gibbs-within-Metropolis–Hastings Markov Chain Monte Carlo (MCMC) algorithm.

Fig. 4

Infectious period of each individual () in the CSFV example. Crosses denote the time of exposure (in days), lines connect the and events.

A Bayesian-analysis of the CSFV data with unknown exposure times could be conducted by the following code: Instead of creating the necessary LBOptionsMCMC object by hand, the call to data also loads an appropriate object laevens.opts for MCMC estimation. For example: shows the requested number of samples to draw from the posterior together with the burn-in and thinning rate; a total of samples*thin+burnin samples are generated. The results are as follows. The results of the MCMC inference are provided as realizations from the Markov chain having the posterior distribution of as stationary distribution. Access to the samples is obtained through the samplePaths method—this makes allowance for further processing using e.g. the coda (convergence diagnostic and output analysis) or boa (Bayesian output analysis) packages from the Comprehensive R Archive Network (Plummer et al., 2006, Smith, 2005). To exemplify, the sampling paths and the marginal posterior of are inspected in Fig. 5 by A matter of particular interest in the CSFV experiment is the relationship between and —especially whether there is a difference in spread within the pen and between neighboring pens. The following code uses the plot method for LBInferenceMCMC objects to generate Fig. 6 and to display the posterior mean together with lower and upper boundaries of a 95%-highest posterior density (HPD) interval for . An important epidemiological quantity of an infectious disease is the basic reproduction ratio , i.e. the expected number of new infections generated by a single infectious individual in a large susceptible population. In a multitype setup with k groups this can be computed as the largest eigenvalue of the matrix containing the mean number of infectious contacts between all units (Andersson and Britton, 2000, Section 6.2). Samples from the posterior distribution of are generated by performing this computation for each posterior sample of . The code stated below uses the method R0 to retrieve the samples for and thus to compute the posterior median together with a symmetric 95% credibility region.

Fig. 5

Sample path and kernel density of the -samples generated using coda.

Fig. 6

Left panel: 2D-Kernel estimated posterior density surface of . Right panel: 1D-Kernel estimate of the posterior density of together with posterior mean, posterior median and a 95%-HPD interval.

Discussion

In this paper we have illustrated the use of RLadyBug  for the simulation, visualization and estimation of infectious disease outbreak data. To our knowledge, the package is the first publicly available software for the estimation in stochastic epidemic models. By providing such specialist functionality within a standard software package as R we hope to make a thorough statistical analysis of infectious disease data a bit more routine.

Many extensions towards more complex and realistic models than provided by the package are imaginable. For example, the general multitype SEIR model (Andersson and Britton, 2000, Chapter 6) allows for more complicated neighbor dependencies: Letting be a matrix given as a function of the parameter vector , the transmission rates are given by . Nearest neighborhood transmission as in (1) thus corresponds to and hence . An alternative would be to let transmission from neighbors be a function of distance: . Modifications of the Java code to handle the simulation and maximum likelihood estimation in such models should be feasible. However, obtaining MCMC estimates in case of missing data would require a substantial amount of work. The same comment applies if one wants to extend beyond the currently implemented waiting time distributions (gamma-distributed, exponential-distributed and constant).

The multitype model could also be applied to the handling of heterogeneous units exemplified by veterinary experiments quantifying the effect of a vaccine. Here, all individuals in specific units are vaccinated, thus having different parameters than individuals in non-vaccinated units. A Monte-Carlo approach could be employed to calculate sample sizes necessary to detect a certain difference in parameters with a given accuracy.

In addition, we are currently working on an implementation of the logistic-regression approach in Klinkenberg et al. (2002). This extension illustrates the benefits of providing a flexible and extensible package: the code is purely R-based exploiting the class structure of the package, while using the optimization routines of R for inference.

Sources, binaries and documentation of RLadyBug  are available for download from the Comprehensive R Archive Network http://cran.r-project.org/ under the GNU Public License. Once installed, the analyses of this article can be reproduced using .

>	library("RLadyBug")
>	layout <- new("LBLayout", S0 = matrix(c(10, 9, 10), 1, 3),
+			E0 = matrix(c(0, 1, 0), 1, 3))
>	options <- new("LBOptions", initBeta = list(init = 0.125),
+		initBetaN = list(init = 0.018), initIncu = list(g = 6.697, d = 0.84),
+		initInf = list(g = 1.772, d = 0.123))
>	plot(simulate(options, layout = layout), type = state^∼time|position)

>E0 <- matrix(0, 8, 16)

>E0[4, 8] <- 1

>S0 <- matrix(15, 8, 16) - E0

>exp <- simulate(options, layout = new("LBLayout", S0 = S0, E0 = E0))

>plot(exp, type = state^∼1|position, options = list(noOfPics = 20))

>data("laevensML")

>plot(laevensML, type = individual^∼time|position)

>data("laevens")

>inf.mcmc <- seir(laevens, laevens.opts)

>algo(laevens.opts)

samples	thin	burnin
2500	25	50000

>inf.mcmc

An object of class LBInferenceMCMC

Parameter Estimations (posterior mean from 2500 samples):
Parameter:
beta	betaN	gammaE	deltaE	gammaI	deltaI
0.03706	0.02837	56.82000	9.37400	2.16200	0.25640
StandardErrors (posterior std.dev. from 2500 samples):
beta	betaN	gammaE	deltaE	gammaI	deltaI
0.018500	0.009481	45.510000	7.761000	0.738100	0.097760

>samples <- mcmc(samplePaths(inf.mcmc))

>plot(samples[, "beta"])

>betabetaN <- plot(inf.mcmc, which = "betabetaN")

>c(mean = betabetaN$mean, betabetaN$hpd)

mean	LB 95% HPD	UB 95% HPD
1.4740612	0.1245430	3.3158915

>quantile(R0(inf.mcmc, laevens), c(0.025, 0.5, 0.975))

2.5%	50%	97.5%
0.3245193	0.6370434	1.2426485