Individual-Level Modelling of Infectious Disease Data: EpiILM

In this article, we introduce the R package EpiILM, which provides tools for simulation from, and inference for, discrete-time individual-level models of infectious disease transmission proposed by Deardon et al. (2010). The inference is set in a Bayesian framework and is carried out via Metropolis-Hastings Markov chain Monte Carlo (MCMC). For its fast implementation, key functions are coded in Fortran. Both spatial and contact network models are implemented in the package and can be set in either susceptible-infected (SI) or susceptible-infected-removed (SIR) compartmental frameworks. The use of the package is demonstrated through examples involving both simulated and real data.


Introduction
The task of modelling infectious disease transmission through a population poses a number of challenges. One challenge is that successfully modelling many, if not most, infectious disease systems requires accounting for complex heterogeneities within the population. These heterogeneities may be characterized by individual-level covariates, spatial clustering, or the existence of complex contact networks through which the disease may propagate. A second challenge is that there are inherent dependencies in infection (or event) times.
To model such scenarios, Deardon et al. (2010) introduced a class of discrete time individuallevel models (ILMs), fitting the models to data in a Bayesian Markov chain Monte Carlo (MCMC) framework. They applied spatial ILMs to the UK foot-and-mouth disease (FMD) epidemic of 2001, which accounted for farm-level covariates such as the number and type of animals on each farm. However, the ILM class also allows for the incorporation of contact networks through which disease can spread. Once fitted, such models can be used to predict the course of an epidemic (e.g., O'Reilly et al., 2018) or test the effectiveness of various control strategies (e.g., Tildesley et al., 2006) that can be imposed upon epidemics simulated from the fitted model.
A third challenge when modelling disease systems is that very little software so far has been made available that allows for simulation from, and especially inference for, individual-level models of disease transmission. Most inference for such models is carried out in fast, low-level languages such as Fortran or variants of C, which makes it difficult for researchers (e.g., public health epidemiologists) without a strong background in computational statistics and programming to make use of the models.
A number of R packages have recently been developed for modelling infectious disease systems (e.g., R0 (Boelle and Obadia, 2015), EpiEstim (Cori, 2019), EpiModel (Jenness et al., 2018), and epinet (Groendyke and Welch, 2016)). Most of these packages can be used to carry out epidemic simulation from given models; in addition, R0 or EpiEstim, for example, can be used to calculate the (basic) reproduction number under various scenarios. The EpiModel package allows for the simulation of epidemics from stochastic models, primarily exponential-family random graph models (ERGMs), and provides tools for analyzing simulation output. Functions for carrying out some limited forms of inference are also provided. Another widely used package for monitoring and modelling infectious disease spread through surveillance data is surveillance (Meyer et al., 2017). This package provides for a highly flexible modelling framework for such data. However, the package does not cover mechanistic, individual-level disease transmission models such as those of Deardon et al. (2010).
Here, we detail a novel R statistical software package EpiILM for simulating from, and carrying out Bayesian MCMC-based statistical inference for spatial and/or network-based models in the Deardon et al. (2010) individual-level modelling framework. The package allows for the incorporation of individual-level susceptibility and transmissibility covariates in models, provides various methods of summarizing epidemic data sets, and permits reasonably involved scenarios to be coded up by the user due to its setting in an R framework. The main functions, including for likelihood calculation are coded in Fortran in order to achieve the goal of agile implementation.
The type of spatial and network-based transmission models that EpiILM facilitates can be used to model a wide range of disease systems, as well as other transmissible processes. Human diseases such as influenza, measles or HIV, tend to be transmitted via interactions which can be captured by contact networks. For example, Malik et al. (2014) used a network representing whether two people shared the same household for modelling influenza spread in Hong Kong. Networks can also be used to characterize social or sexual relationships.
In the livestock industries, diseases are often transmitted from farm to farm via supply trucks or animal movements from farm to farm, or from farm to market. For example, ILM's were used by arXiv:2003.04963v1 [stat.AP] 10 Mar 2020 Kwong et al. (2013) to model the spread of porcine reproductive and respiratory syndrome (PRRS) through Ontario swine farms via such mechanisms. Spatial mechanisms are also often important in livestock industries (e.g., Jewell et al., 2009;Deardon et al., 2010;Kwong et al., 2013), as well as for modelling crop diseases (e.g., Pokharel and Deardon, 2016), since airborne spread is often a key factor.
Further, these types of models can also be used to model transmissible processes other than infectious disease spread. For example, Cook et al. (2007) used similar models to model the transmission of alien species through a landscape; specifically, giant hogweed in the UK. In addition, Vrbik et al. (2012) used spatial ILM's to model fire spread. They looked at fire spread under controlled conditions, but such models would likely be useful for modelling the spread of forest fires since important covariates such as vegetation-type could be incorporated into the models.
Data from infectious disease systems are generally 'time-to-event', typically involving multiple states. However, standard survival models (e.g., Cox (1972), Therneau (2015)) or multi-state time-toevent models (e.g., see Jackson (2011)) are not applicable here, because in an infectious disease system individual event times cannot be assumed independent even after conditioning on covariates. That is, my risk of contracting and infectious disease generally depends upon the disease state of other individuals in the population; this is not typically the case for most cancers, for example to which more standard models can be applied.
The remainder of this paper is structured as follows: Section 2 explains the relevant models involved in the package; Section 3 describes the contents of the package along with some illustrative examples; and Section 4 concludes the paper with a brief discussion on future development.

Model
In our EpiILM package, we consider two compartmental frameworks: susceptible-infectious (SI) and susceptible-infectious-removed (SIR). In the former framework, individuals begin in the susceptible state (S) and if/when infected become immediately infectious (I) and remain in that state indefinitely. In the latter framework, individuals once infected remain infectious for some time interval before entering the removed state (R). This final state might represent death, quarantine, or recovery accompanied by immunity. We consider discrete time scenarios so a complete epidemic history is represented by t = 1, 2, . . . , t end , where (typically) t = 1 is the time when the first infection is observed and t end is the time when the epidemic ends. Hence, for a given time point t, an individual i belongs to one, and only one, of the sets S(t) or I(t) if the compartmental framework is SI, and i belongs to one, and only one, of the sets S(t), I(t), or R(t) if the compartmental framework is SIR.
Under either framework, the probability that a susceptible individual i is infected at time point t is given by IP(i, t) as follows: where: Ω S (i) is a susceptibility function that accommodates potential risk factors associated with susceptible individual i contracting the disease; Ω T (j) is a transmissibility function that accommodates potential risk factors associated with infectious individual j contracting the disease; ε is a sparks term which represents infections originating from outside the population being observed or some other unobserved infection mechanism; and κ(i, j) is an infection kernel function that represents the shared risk factors between pairs of infectious and susceptible individuals.
The susceptibility function can incorporate any individual-level covariates of interest, such as age, genetic factors, vaccination status, and so on. In Equation (1), Ω S (i) is treated as a linear function of the covariates, i.e., Ω S (i) = α 0 + α 1 X 1 (i) + α 2 X 2 (i) + · · · + α n s X n s (i), where X 1 (i), . . . , X n s (i) denote n s covariates associated with susceptible individual i, along with susceptibility parameters α 0 , . . . , α n s > 0. Note that, if the model does not contain any susceptibility covariates then Ω S (i) = α 0 is used. In a similar way, the transmissibility function in Equation (1) can incorporate any individual-level covariates of interest associated with infectious individual. Ω T (j) is also treated as a linear function of the covariates, but without the intercept term, i.e., Ω T (j) = φ 1 X 1 (j) + φ 2 X 2 (j) + · · · + φ n t X n t (j), where X 1 (j), . . . , X n t (j) denote the n t covariates associated with infectious individual j, along with transmissibility parameters φ 1 , . . . , φ n t > 0. Also note that if the model does not contain any transmissibility covariates then Ω T (j) = 1 is used.
In this package, we also consider two broad types of ILM models based on the type of the kernel function κ(i, j): spatial and network-based ILMs. In the spatial-based ILMs, the infection kernel function is represented by the power-law function as where β is the spatial parameter that accounts for the varying risk of transmitting disease over the Euclidean distance between individuals i and j, d ij . Whereas in the network-based ILMs, κ(i, j) can be represented by one or more contact network matrices and is written as ij denotes the (i, j) th element of what we term the contact matrix of a given contact network; in graph theory this is more typically referred to as a (weighted) adjacency matrix. The corresponding β (.) 's represent the effect of each of the n networks on transmission risk. In each contact network, each individual in the population is denoted by a node and is connected by lines or edges. These connections represent potential transmission routes through which disease can spread between individuals in the population. If the network is unweighted, the contact matrix is treated as binary (0 or 1). If the edges have weights assigned to them, then C ij ∈ [0, 1] are typically used. These weights can be used to allow for different infection potential between different pairs of individuals. If the network is undirected, the contact matrix will be symmetric; if directed, it can be non-symmetric. Finally, the C ii (diagonal terms) are not used in the models and are typically set to C ii = 0, ∀ i.
Note that IP(i, t) gives the probability that susceptible individual i is infected at time point t, representing some interval in continuous time (e.g., a day or week), but they actually become infectious at time t + 1.
Following Deardon et al. (2015), the likelihood function for the ILMs (1) is given by and where, θ is the vector of unknown parameters, I(t + 1)\I(t) denotes all new infections observed at t + 1 in the infectious state at time t, and t max ≤ t end is the last time point at which data are observed or being simulated.

Contents of EpiILM
The EpiILM package makes use of Fortran code that is called from within R. This package can be used to carry out simulation of epidemics, calculate the basic reproduction number, plot various epidemic summary graphics, calculate the log-likelihood, and carry out Bayesian inference using Metropolis-Hastings MCMC for a given data set and model. The functions involved in the package are summarized in Table 1.

Simulation of epidemics
The function epidata() allows the user to simulate epidemics under different models and scenarios. One can use the argument type to select the compartmental framework (SI or SIR) and population size through the argument n. If the compartmental framework is SIR, the infectious period is passed through the argument infperiod. Depending on whether a spatial or network model is being considered, the user can pass the arguments: x, y for location and contact for contact networks. Users can also control the susceptibility function Ω S (i) through the Sformula argument, with individuallevel covariate information passable through this argument. If there is no covariate information, Sformula is null. An expression of the form Sformula = ∼ model is used to specify the covariate information, separated by + andoperators similar to the R generic function formula(). For example, Ω S (i) = α 0 + α 1 X(i), i = 1, . . . , n can be passed through the argument Sformula as Sformula = ∼ 1 + X. In a similar way, the user can control the transmissibility function Ω T (i) through the Tformula argument. Note that, the Tformula must not include the intercept term to avoid model identifiability issues, i.e., for a model with one transmissibility covariate (X), the Tformula becomes Tformula = ∼ -1 + X. The spatial (or network), susceptibility, transmissibility, and spark (if any) parameters are passed through arguments beta, alpha, phi, and spark, respectively.
The argument tmin helps to fix the initial infection time while generating an epidemic. By default, tmin is set as time t = 1. We can also specify the initial infective or infectives using the argument inftime. For example, in a population of 10 individuals, we could choose, say, the third individual to become infected at time point 1, using the option inftime = c(0,0,1,0,0,0,0,0,0,0). We could  also infect more than one individual and they could be infected at different time points. This allows simulation from a model conditional on, say, data already observed, if we set the tmin option at the maximum value of inftime.
The output of the function epidata() is formed as class of epidata object. This epidata object contains a list that consist of type (the compartmental framework), XYcoordinates(the XY coordinates of individual for spatial model) or contact (the contact network matrix for the network model), inftime (the infection times) and remtime (the removal times). Other functions such as plot.epidata and epimcmc involved in the package use this object class as an input argument.

Descriptive analyses
We introduce an S3 method plot function to graphically summarize, and allow for a descriptive analyses of, epidemic data. The function plot.epidata() illustrate the spread of the epidemic over time. One of the key input arguments (x) of this function has to be an epidata object. The other key argument plottype has two options: curve and spatial. Specifying the first option produces various epidemic curves, while the latter show the epidemic propagation over time and space when the model is set to spatial-based. When the plottype = curve, an additional argument needs to be passed through the function through curvetype. This has four options: curvetype = "complete" produces curves of the number of susceptible, infected, and removed individuals over time (when type = "SIR ); "susceptible" gives a single curve for the susceptible individuals over time; "totalinfect" gives the cumulative number of infected individuals over time; and "newinfect" produces a curve of the number of newly infected individuals at each time point.
The plot functions plot.epimcmc() and plot.pred.epi() can be used to illustrate inference results (see Bayesian inference section). Detailed explanation is provided in the corresponding subsections.

Example: spatial model
Suppose we want to simulate an epidemic from the model (1) using a spatial kernel with type SI, Ω s (i) = α, and no transmissibility covariates, Ω T (j) = 1. Choosing the infectivity parameter α = 0.3, spatial parameter β = 5.0, sparks parameter ε = 0, and t max = 15, the model is given by: where d ij is the Euclidean distance between individuals i and j, with their locations specified through x and y.
First, we install the EpiILM package and call the library.

Example: network model
To illustrate simulation from a contact network-based model, we consider a disease system in which disease transmission can occur through a single, directed, binary network over a population of n = 100 individuals. The elements in the contact matrix represent the existence or non-existence of a directed connection through which disease can transmit between two individuals in the population. Each individual within the population is represented by a row and column within the matrix. Specifically, an element C ij in the contact matrix is given by: if a directed edge exists between i and j 0 otherwise (5) We also consider the inclusion of a binary susceptibility covariate Z in the model. This can be thought to represent, say, treatment or vaccination status. The infection model is then given by where α 0 is the baseline susceptibility and α 1 is the binary treatment effect. The parameters in the model are set to be (α 0 , α 1 ) = (0.1, 0.05) and ε = 0, and we also set t max = 15. We simulate a directed network using the following code: R> contact <-matrix(rbinom(10000, 1, 0.1), nrow = 100, ncol = 100) R> diag(contact[, ]) <-0 Various packages are available in R for network visualization, such as igraph, ergm, etc and as an example, we use the igraph package for the network display as shown in Figure 3.

Bayesian Inference
In EpiILM, ILMs can be fitted to observed data within a Bayesian framework. A Metropolis-Hastings MCMC algorithm is provided which can be used to estimate the posterior distribution of the parameters. The function epimcmc() provides three choices for the marginal prior distribution of each parameter: the gamma, half-normal, and uniform distributions. The parameters are assumed to be a priori independent. The proposal used is a Gaussian random walk. Again, users can use the Sformula and Tformula arguments to specify any individual-level susceptibility and transmissibility covariates. Users can also control the number of MCMC simulations, initial values, and proposal variances of the parameters to be estimated. Note that in case of fixing one parameter and updating other parameters, users can do it by setting the proposal variance of fixed parameter to zero. This is usually the case to avoid identifiability issue when the model has both susceptibility and transmissibility covariates without intercept terms. Again, spatial/network, susceptibility parameters and transmissibility parameters are passed through arguments beta, sus.par and trans.par, respectively. One can specify the spark parameter using the spark argument, but by default its value is 0. epimcmc() can also call the adaptive MCMC method of inference facilitated by the adaptMCMC package. Specifically, users can pass the argument adapt = TRUE along with acc.rate to run the adapt MCMC algorithm.
We can use the S3 method functions summary.epimcmc() and plot.epimcmc() available in the package for output analysis and diagnostics. Both are dependent upon the coda package. The argument plottype in the function plot.epimcmc, has two options to specify which samples are to be plotted: (1) "parameter" is used to produce trace plots (time series plots) of the posterior distributions of the model parameters, and (2) "loglik" to produce trace plots of the log likelihood values of the model parameter samples. Other options that are used in the coda package can be used in the plot.epimcmc() as well; e.g., start, end, thin, density, etc.
Spatial or network-based ILMs can be fitted to data as shown in the following examples.

Example: spatial model
Suppose we are interested in modelling the spread of a highly transmissible disease through a series of farms, say n = 100, along with an assumption that the spatial locations of the farms and the number of animals on each farm are known. In this situation, it is reasonable to treat the farms themselves as individual units. If we treat the extent of infection from outside the observed population of farms as negligible (ε = 0), we can write the ILM model as where the susceptibility covariate A represents the number of animals on each farm, α 0 is the baseline susceptibility, α 1 is the number of animals effect, and β is the spatial parameter. Let us use the same (simulated) spatial locations from the previous spatial model example and set the parameters (α 0 , α 1 ) = (0.2, 0.1) and β = 5. We also set t max = 50. Considering an SI compartmental framework for this situation, the epidemic is simulated using the following command: R> A <-round(rexp(100,1/50)) R> SI.dis.cov <-epidata(type = "SI", n = 100, tmax = 50, x = x, y = y, + Sformula =~A, sus.par = c(0.2, 0.1), beta = 5) We can now refit the generating model to this simulated data and consider the posterior estimates of the model parameters. We can do this using the following code: R> t_end <-max(SI.dis.cov$inftime) R> unif_range <-matrix(c(0, 0, 10000, 10000), nrow = 2, ncol = 2) R> mcmcout_Model7 <-epimcmc(SI.dis.cov, Sformula =~A, tmax = t_end, niter = 50000, + sus.par.ini = c(0.001, 0.001), beta.ini = 0.01, + pro.sus.var = c(0.01, 0.01), pro.beta.var = 0.5, + prior.sus.dist = c("uniform","uniform"), prior.sus.par = unif_range, + prior.beta.dist = "uniform", prior.beta.par = c(0, 10000)) where niter denotes the number of MCMC iterations and sus.par.ini and beta.ini are the initial values of the parameters to be estimated. The proposal variances for (α 0 , α 1 ) and β are set to (0.01, 0.01) and 0.5, respectively (after tuning). Vague uniform prior distributions are used for all three parameters, i.e., we choose U(0, 10000) for α 0 , α 1 , and β. As the locations x and y are specified, a spatial ILM is fitted rather than a network-based ILM. Note that the full data set to which the model is being fitted consists of the spatial locations (x, y) and the infection times. Figure 5 displays the MCMC traceplot after 10000 burn-in using the command R> plot(mcmcout_Model7, partype = "parameter", start = 10001, density = FALSE) The posterior means and 95% credible intervals (CI) of the parameters, calculated as the 2.5% and 97.5% percentiles of 50000 MCMC draws after a burn-in of 10000 iterations has been removed, can be obtained using the following code: R> summary(mcmcout_Model7, start = 10001)

Example: network model
Now consider modelling the spread of an animal infectious disease through a series of farms (n = 500) in a region. We once again consider individuals in the population to be the farms, rather than animals themselves, and a single binary contact network to represent connections between farms. We assume that there are two species of animals in farms that play an important role in spreading the disease; let us say, cows and sheep. Thus, we include the number of cows and sheep on the farm as susceptibility and transmissibility covariates in the model. Thus, the probability of susceptible farm i to be infected at time t becomes: with susceptibility parameters, (α 1 , α 2 ), transmissibility parameters, (φ 1 , φ 2 ), and X 1 and X 2 represent the number of sheep and cows in farms, respectively. Note, when we have a single network, there is no network parameter to estimate; this is to avoid problems of non-identifiability.

Conclusion
This paper discusses the implementation of the R software package EpiILM. Other than this package, there does not appear to be any R software that offers spatial and network-based individual-level modelling for infectious disease systems. Thus, this package will be helpful to many researchers and students in epidemiology as well as in statistics. These models can be used to model disease systems of humans (e.g., Malik et al. (2014)), animals (e.g., Kwong et al. (2013)), or plants (e.g., Pokharel and Deardon (2016)), as well as other transmission-based systems such as invasive species (e.g., Cook et al. (2007)) or fire spread (Vrbik et al., 2012). The EpiILM package continues to exist as a work in progress. We hope to implement additional models and options in the future, that might be useful for other researchers in their research and teaching. Such additions may include the incorporation of time-varying networks, covariates and/or spatial kernels (e.g., Vrbik et al. (2012)), models that allow for a joint spatial and network-based infection kernel, uncertainty in the times of transitions between disease states (e.g., Malik et al. (2016)), unknown covariates (e.g., Deeth and Deardon (2013), and extensions to other compartmental frameworks such as SEIR and SIRS. Finally, we hope to extend the package to allow for the modelling of disease systems with multiple interacting strains or pathogens (Romanescu and Deardon, 2016).