The vast majority of epidemiological models, such as the well-known susceptible-infectedrecovered (SIR) model, rely on compartmentalizing individuals according to their disease status, meaning that individual level information is lost (Kermack and McKendrick, 1927) .
Classical compartmental models adopt the "random mixing" assumption under which an infectious individual can spread the disease homogeneously to any susceptible individual (Kermack and McKendrick, 1927; Bailey et al., 1975; Anderson and May, 1992) . In the last two decades, an alternative assumption-that the disease is transmitted through links in a contact network-has gradually gained popularity. As a result, there has been an increase in the relevance and importance of networks in transmittable disease research.
Networks can flexibly and informatively represent human interactions, contact patterns, and social structures, which can fundamentally impact the behavior of epidemic processes (Wallinga et al., 1999; Edmunds et al., 1997 Edmunds et al., , 2006 Mossong et al., 2008; Melegaro et al., 2011) . Also, there is evidence that the population network structure can be influenced by infection events as well (Bell et al., 2006; Funk et al., 2010; Eames et al., 2010; Van Kerckhove et al., 2013) . Numerous recent epidemiological methods consider both disease transmission and the contact network as dynamic processes (Bansal et al., 2010; Pastor-Satorras et al., 2015; Masuda and Holme, 2013; Enright and Kao, 2018; Masuda and Holme, 2017) , and even as coupled processes: disease spread is impacted by the evolving network structure, which also adapts to contagion progression in that individuals with different health statuses adopt different social behaviors (Kiss et al., 2012; Ogura and Preciado, 2017; Van Segbroeck et al., 2010; Tunc et al., 2013; Group, 2006) .
In this paper, we propose a stochastic generative model that jointly describes SIR-type epidemic processes on temporal networks and dynamics of adaptive networks, and then develop a likelihood-based Bayesian inference scheme that accommodates missingness in individual epidemic recovery times. We choose stochastic models over deterministic models widely used in existing literature to account for the randomness in small-scale epidemic outbreaks and allow for measures of uncertainty in estimation (Britton, 2010) . The model formulation is also motivated by the increasing need for inference based on real-world epi-demiological data that contain time-resolved social interactions data. Almost all previous methods for adaptive network epidemic processes focus on theoretical derivation or numerical simulation due to the long-standing difficulty in dynamically tracking social contacts.
However, thanks to technological advancements in mobile devices, several recent observational studies have collected high-resolution data on social contacts (Vanhems et al., 2013; Barrat et al., 2014; Voirin et al., 2015; Kiti et al., 2016; Aiello et al., 2016; Ozella et al., 2018) , which enables inference on the interplay between epidemic and network processes.
Moreover, real-world epidemiological data almost always contain some missingness, particularly in the exact times of individual disease episodes. This is a major challenge for inference, even without the additional complexity of network dynamics (Ho et al., 2018) .
To address this challenge, we design an inference scheme that includes a novel data augmentation algorithm that leverages both the network dynamics and the generative process implied by the model. This algorithm is, in principle, similar to existing agent-based data augmentation methods for fitting stochastic epidemic models that propose individual disease episodes compatible with the observations (Cauchemez et al., 2006; Hoti et al., 2009; Fintzi et al., 2017) . It also bears some resemblance to simulation-based methods (Andrieu et al., 2010; Pooley et al., 2015) that use the underlying model to propose epidemic paths (yet at the population-level) as the foundation for inference. However, these methods are inadequate in coping with the intricate constraints imposed by time-varying social contacts, or utilizing network information to reduce the difficulty in evaluating the exact likelihood function. Moreover, these methods are computationally expensive, especially those based on simulations. In comparison, the data augmentation method introduced in this work is more effective and efficient. This paper is organized as follows: Section 2 provides some background knowledge on epidemic models (the stochastic SIR model in particular), basic network concepts, and dynamic and adaptive network processes. Section 3 formulates the generative model and derives maximum likelihood estimators as well as Bayesian posterior distributions based on the complete data likelihood. Section 4 describes a Bayesian inference scheme that deals with incomplete observations on individual recovery times. Section 5 and 6 present experiment results on simulated datasets and a real-world dataset. Finally Section 7 gives further discussions.
The vast majority of epidemiological models are based on compartmentalizing individuals into non-over lapping subsets according to their disease statuses. In classical models, the sizes of those subsets are described by differential equations (Hethcote, 2000) . One famous and widely used model is the susceptible-infected-recovered (SIR) model, which assumes three disease statuses-susceptible (S), infected (or infectious, I), and recovered (or removed, R). On a closed population of N individuals (with N sufficiently large), the deterministic dynamics of the SIR model can be expressed as
where S(t) and I(t) refer to the number of susceptible and infected individuals at time t,
respectively, and the number of recovered individuals satisfies R(t) = N − [S(t) + I(t)]. Here β is the rate of disease transmission per contact between an S individual and an I individual, and γ is the rate of recovery for an I individual.
By setting the growth rate of infection to be proportional to S(t)I(t), the model in (1) implicitly assumes that any two members can interact with each other. This assumption is easily violated in reality, where an individual only maintains contact with a limited number of others. Moreover, the differential equations can only account for the average, expected behavior of the process, but the transmission of an infectious disease exhibits randomness and uncertainty by nature.
To account for the underlying network structure of a population as well as the random nature of an epidemic process, we adopt an individual-level, stochastic variation of the SIR model (similar to that used in Auranen et al. (2000) ). An individual of status S (susceptible) at time t (> 0) changes disease status to I (infected/infectious) at time t + h (h is an infinitesimal quantity) with a probability that is dependent on both the infection rate β and his/her contacts at time t. An infected individual at time t becomes a member of the R (recovered) sub-population at time t + h with a probability determined by the recovery rate γ. Specifically, for any susceptible individual p 1 and infected individual p 2 in the population at t, conditioned on the current overall state of the process, Z t , then P r(p 1 gets infected by p 2 at t + h|Z t ) = βh + o(h)
if p 1 and p 2 are in contact at t, and P r(p 2 recovers at t + h|Z t ) = γh + o(h).
(3)
A network, or a graph, is a two-component set, G = {V, E}, where V is the set of N nodes and E is the set of links. A network can be represented by its "adjacency matrix", A, where A ij = 1 indicates there is a link from node i to j. Since most infectious diseases can be transmitted in both directions through a contact, we assume that the adjacency matrix is
A special network structure is the fully connected network (or the complete graph), K N , and its adjacency matrix A satisfies A ij = 1 for any i = j. This network structure corresponds to the widely adopted "random mixing" assumption in epidemiological models, which, as stated before, is unrealistic and restrictive. Therefore, in the rest of the paper, we consider arbitrary network structures on a population instead.
Interactions between individuals are dynamic in nature, and such dynamics is important when modeling epidemic processes (Masuda and Holme (2017) ; also as demonstrated later in Section 5.1). We consider a continuous-time link-Markovian model for temporal networks (Clementi et al., 2010; Ogura and Preciado, 2016) . For two individuals i and j who are not in contact at time t, they form a link at time t + h (h 1) with probability αh, where α is the link activation rate. Similarly, if there is an edge between i and j at time t, then the edge is deleted at time t + h with probability ωh, where ω is some link termination rate.
If, instead, individuals establish and dissolve their social links with rates that vary according to their disease statuses, then the evolution of the network is coupled to the epidemic process and thus becomes adaptive. This mechanism can be described by instantaneous rates of single-link changes. For any two individuals i and j, their corresponding entry in the adjacency matrix can be thought of as a {0, 1}-valued Markov process, A ij (t), t > 0. Suppose that at time t, i is of status A, j is of status B, 1 then for an infinitesimal quantity h,
3 Adaptive Network Epidemic Processes and Inference with Complete Data
In this subsection, we lay out a generative model for the joint evolution of an individualized SIR process on a networked population and the dynamics of the contact network. The key feature of the model is the interplay between epidemic progression and network adaptation:
on one hand, transmission of infection depends on the existence of susceptible-infected links, which may change through time; on the other hand, network links temporally update in a manner that in turn depends on individual disease status.
We formulate this complex process as a continuous-time Markov process that comprises all the competing Poisson processes with exponentially-distributed wait times at the individual level, as described in (2)-(5). Throughout the process, the system is updated stochastically with one event (for one individual or a pair) at a time, and each event is of one of the four types:
• Infection: The disease is transmitted through a link between an S (susceptible) and
an I (infected) individual (S-I link) with rate β;
• Recovery: Each I individual recovers with rate γ independently;
• This model formulation allows for joint inference of both disease spread and network evolution, and as illustrated in the next sub-section, inference is straightforward when all the events are fully observed. Furthermore, this formulation corresponds to a relatively simple generative process at the population level too. By the superposition property, conditioned on the current state of the process Z t at time t (> 0), the very next event of the entire process is the earliest event that occurs among the four competing processes:
• Infection: An infection occurs with rate βSI(t), where SI(t) is the number of S-I links at time t;
• Recovery: A recovery occurs with rate γI(t), where I(t) is the number of infected individuals at time t;
• Link activation:
This generative model is a generalization of two simpler models. If we set α AB ≡ α and ω AB ≡ ω for any status A and B, the coupled process is reduced to a decoupled process, where network evolution is independent of individual disease status. Moreover, if we fix α ≡ ω ≡ 0, the process is further reduced to an SIR process over a static network.
Here we assume that the population size N is fixed, and both the initial network structure G 0 and the number of infection cases at onset, I(0), are known. We summarize a list of important parameters and notation used for model formulation and inference in Table 1 . 
Derivation of complete data likelihood Let i k be the infection time for individual k (k = 1, 2, . . . , n E ), r k be the k th observed recovery time (k = 1, 2, . . . , n R ), and without loss of generality, set i 1 = 0. Recall that the widely used "random mixing" assumption in classical epidemiological models is equivalent to assuming that the contact network is a complete graph, K N , and the complete data likelihood under this assumption is
Instead, assume that the network is an arbitrary static network G. Then when writing down the complete data likelihood, we have to explicitly account for the number of infected contacts for each individual at the time of infection as well as the total number of S-I links in the system:
Here I k (t) denotes the number of infectious individuals that are connected to person k at time t, which remains fixed for a static network.
If the network is not static, but the entire network process {G t : 0 < t < T max } is given (or all the network changes during the contagion process are fully observed), then the form of the data likelihood remains unchanged, conditioned on {G t }. Note that the dynamic nature of the network can be implicitly subsumed into the terms I k (i k )'s and SI(u).
To incorporate network dynamics, we begin with the simpler decoupled process in which the network evolves independently of the epidemic process with edge activation rate α and deletion rate ω. Given the initial network G 0 , the likelihood of the network process alone can be easily written down as
Here s is the time of the th network event, and D = 1 if this event is a link termination and otherwise D = 0. Therefore, when the epidemic process and network process are decoupled, the complete data likelihood is simply a product of the terms in (6) and (7):
Then we consider the coupled process with an adaptive network, where link activation and termination are dependent on individual disease status. Define g(p, t) as the indicator function of infectiousness, i.e. g(p, t) = 1 if person p is infected at time t and g(p, t) = 0 otherwise. We further assume that the S and R populations behave identically in the network process:
We shall refer to them as the "H" (healthy) population; let
denote the number of healthy individuals at time t. Naturally the term "H-H link" represent an S-S link, an S-R link, or an R-R link, and the term "H-I link" represent either an S-I link or an R-I link.
Combine the epidemic events and network events and re-denote all events as {e j = (t j , p j1 , p j2 )} n j=1 , n = n E + n R + n N . Here t j is the event time (set t 1 = 0, the infection time of the first patient), with t 1 < t 2 < . . . < t n . If e j is a network event, p j1 , p j2 are the two individuals getting connected or disconnected, and if e j is an epidemic event, let p j1 be the person getting infected or recovered and set p j2 = 0. Furthermore let event type indicators F j , C j , D j take the value 1 only if e j is an infection, a link activation, and a link deletion, respectively, and otherwise take the value 0.
The contribution of all network events to the complete data likelihood is in essence of the same form as (7), except that for every activation or termination event the link type has to be considered. Then the likelihood component of the adaptive network process is
Therefore, given the initial network structure G 0 and one infected case at time 0, the complete data likelihood can be expressed as
Inference Given Complete Event Data The derived likelihood function in (14) indicates that inference should be relatively straightforward when all the changes in the process are fully observed. Given the complete event data {e j } n j=1 and the initial conditions of the process G 0 and I(0), the only unknown quantities in (14) are the model parameters Θ = {β, γ, α SS , α SI , α II , ω SS , ω SI , ω II }. In the following theorems, we state results on maximum likelihood estimation as well as Bayesian estimation.
Theorem 3.1 (Maximum likelihood estimation). Following the likelihood function in (14),
given G 0 and complete event data {e j }, the MLEs of the model parameters are given as
.
The above results can be directly obtained by setting all partial derivatives of the loglikelihood to zero. The detailed proof is provided in Supplement S2.
Theorem 3.2 (Bayesian inference with conjugate priors). If Gamma priors are adopted for all the parameters:
Then following the likelihood function in (14) and given G 0 as well as the complete data {e j }, the posterior distributions of the parameters are given by
whereβ,γ,α SS ,α SI ,α II ,ω SS ,ω SI ,ω II are the MLEs for the parameters defined in Theorem 3.1.
Note that (14) implies that each parameter is essentially the rate of Exponential wait times between consecutive events of a continuous-time Markov process. Applying the Gamma-Exponential conjugacy leads to the posterior distributions in (15).
Relaxing the closed population assumption Throughout this section, the host population is implicitly assumed to be closed. If, instead, the observed population is in fact a subset (with its size N still fixed) of a larger yet unobserved population, then it is possible for an individual to get infected by an external source that is not under study. This external infection force cannot be reflected by the per S-I link infection rate β, so we introduce an "external infection" rate ξ as the rate for each susceptible individual to contract the disease from an external source. ξ can also be thought of as the constant rate for every S individual to spontaneously enter status I, independently of the influence from the infected members inside the population. In this scenario the complete data likelihood becomes
The MLEs for {γ,α,ω} remain unchanged, and although there isn't a closed-form solution for the MLEs of β and ξ, numerical solutions can be easily obtained.
Moreover, if we have information on which infection cases are caused by internal sources population-wide prevalence counts (Fintzi et al., 2017) , and observations on a structured (but static) population (Neal and Roberts, 2004; O'Neill, 2009; Tsang et al., 2019 ), but has not yet been designed for an epidemic process coupled with a dynamic network. The timevarying nature of social interactions imposes complex constraints on the implementation of data augmentation, which already involves difficult conditioned sampling, so previous techniques can be quite ineffective or inefficient. Yet, network dynamics also provides information on possible infection sources and transmission routes that should be exploited.
In this section, we derive a data augmentation method specifically designed to enable inference under missing recovery times. The algorithm utilizes the information presented by the dynamic contact structure. In contrast to existing methods (such as Fintzi et al. (2017) and Hoti et al. (2009) ), it is able to efficiently impute unobserved event times in parallel instead of updating individual trajectories one by one.
Problem setting Our goal is to conduct inference despite the absence of exact times of recovery events in the observed data. Throughout the observation period (0, T max ], suppose
is the collection of disjoint time intervals in which a certain number of recoveries occur, but the exact times of those recoveries are unknown.
That is, for each = 1, 2, . . . , L, some individuals are reported as infected up to time u , and they are reported as healthy again starting from time v . Within one particular interval (u , v ] , let n ( ) E be the number of infections, and n ( ) R the number of recoveries for which the exact times are known, so the number of unknown recovery times for this interval is
Label the R individuals with unknown recovery times by {k ,1 , . . . , k ,R }, and denote their recovery times by latent variables {r ,1 , . . . , r ,R }.
Further assume that we have a health status report (indicating ill or healthy) of each individual periodically during (0, T max ]. Access to such information is usually granted in epidemiological studies where every study subject gives updates on health statuses through regular surveys (for example, weekly surveys). Let S be the maximum number of iterations, then for s = 1 : S, do: 
In the inference scheme stated above, the data augmentation step (step 1) is challenging. This is because (16) describes the distribution of missing recovery times conditioned on both historical events and future events, and thus sampling from (16) is in fact sampling unobserved event times from a continuous-time Markov process with a series of fixed endpoints (Hobolth and Stone, 2009 ). Even though (3) suggests that, in forward simulations, the time it takes for an infected person to recover only depends on the recovery rate γ and is completely independent of the individual's social links or the epidemic history of any other individual, when recovery times need to be inferred retrospectively, there are additional constraints imposed by the observed data. First of all, an individual q cannot recover before a certain time point t if it is observed that at time t the person is still ill. Second, and more subtly, if another individual p gets infected during his contact with q, then q shouldn't have a recovery time that leaves p without any possible infection source. The first condition is easy to satisfy. The second constraint, in comparison, is much more complicated due to the network dynamics, which a simple forward simulation approach would fail to effectively accommodate.
We tackle the challenge in data augmentation by first simplifying the expression of (16) and then stating an efficient sampling algorithm.
Lemma 4.1. (16) can be simplified into the following expression:
where Z t is the state of the process at time t, including the epidemic status of each individual and the social network structure.
Proof. Consider the joint density of the complete data given parameter values Θ (s−1)
.
Examining all terms concerning {r ,i } i=1:R for each leads to
The lemma above suggests that imputation of missing recovery times inside an interval (u, v] only depends on the events that occur in (u, v] , the state of the process at the start of the interval, Z u , and the value of recovery rate γ.
Now consider sampling recovery times within any interval (u, v] . Let Q denote the group of individuals who recover at unknown times during (u, v] , and for each q ∈ Q, let q's exact recovery time be r q (∈ (u, v]). Similarly, let P denote the group of individuals who get infected during (u, v] ; for p ∈ P, let p's infection time be i p , N p (i p ) be the set of p's contacts at time i p , and I(i p ) be the set of known infected individuals at time i p (that is, I(i p ) excludes any individual who may have recovered before i p ). Proposition 4.2 (Data augmentation regulated by contact information (DARCI)). Following the notation stated above, given a recovery rate γ, the state of the process at time u, Z u , and all the observed events in the interval u,v] , Z u in the following steps:
1. Initialize a vector LB of length |Q| with LB q = u for every q ∈ Q;
2. Arrange the set P in the order of {p 1 , p 2 , . . . , p |P| } such that i p 1 < i p 2 < . . . < i p |P| , and for each p ∈ P (chosen in the arranged order), examine the set
If I p ⊂ Q, then randomly and uniformly select one q ∈ I p , and set LB q = i p .
3. Draw recovery times r q ind ∼ TEXP(γ, LB q , v), where TEXP(γ, s, t) is a truncated Exponential distribution with rate γ and truncated on the interval (s, t).
The intuition is to draw a sample of recovery times that are "consistent with" the observed data. To achieve this goal, an imputed recovery cannot occur in a way that leaves a tobe-infected individual without any infectious neighbor at the time of infection, nor take place before the corresponding individual gets infected. Effectively there is a "lower bound"
for each missing recovery time conditioned on the observed data, particularly the dynamic contact structure. An illustration of the DARCI algorithm is provided in Figure 1 .
Combining Lemma 4.1 and Proposition 4.2 enables exact sampling from the conditional distribution (17) in the data augmentation step: for each = 1, 2, . . . , L, applying the DARCI algorithm to the interval (u , v ] gives an updated set of missing recovery times,
This makes the MCMC sampling scheme simply a Gibss sampler.
In this section we present results of a series of experiments with simulated datasets. In all experiments, we employ a forward simulation procedure that can be seen as a variation of Gillespie's algorithm (Gillespie, 1976) in order to generate realizations of the dynamic network epidemic process. The input consists of the parameter values Θ = {β, γ,α,ω} 2 , an arbitrary initial social network G 0 , the number of infectious cases at onset I(0), and the maximum observation time T max . The output is the complete collection of all the events that occur within the time interval (0, T max ], {e j = (t j , p j1 , p j2 , F j , C j , D j )}, and for each event e j , there is a timestamp (t j ), labels of the individuals involved (p j1 , p j2 ), and the type of the
The steps of the simulation procedure are detailed as follows:
1. Initialization. Randomly select I(0) individuals to be the infected/infectious (then the rest of the population are all susceptible). Set t cur = 0.
2. Iterative update. While t cur < T max , do:
(a) Bookkeeping. Summarize the following statistics at t cur : 1) SI(t cur ), the number of S-I links in the population; 2) M max (t cur ), the possible number of links of each type defined in (12); 3) M(t cur ), the number of existing links of each type defined in (13).
Compute the instantaneous rate of the occurrence of any event, Λ(t cur ) = βSI(t cur ) + γI(t cur ) +α T M d (t cur ) +ω T M(t cur ), and draw ∆t ∼ Exponential(Λ(t cur )).
(c) Next event type. Sample Z ∼ Multinomial(λ(t cur )), wherẽ
Then do one of the following based on the value of Z: proportional toω•M(t cur ), and uniformly pick one existing "Y link" to terminate.
(d) Replace t cur by t cur + ∆t, record relevant information about the sampled event, and repeat from (a).
In
Step 2 (c), "•" refers to the Hadamard product (entrywise product) for two vectors.
In this subsection, we first demonstrate the insufficiency of analyzing disease spread without considering the network structure or its dynamics. Then we validate our claims on maximum likelihood estimation and Bayesian inference given complete event data (Theorems 3.1 and 3.2). Finally, we show that the model estimators can detect simpler models such as the decoupled process and the static network process.
Unless otherwise stated, we set the initial network G 0 as a random Erdős-Rényi graph (undirected) with edge probability p = 0.1, let I(0) = 1 individual to be infected at onset, and choose the ground-truth parameters as β = 0.03, γ = 0.12;α T = (0.005, 0.001, 0.005),ω T = (0.05, 0.1, 0.05).
For Bayesian inference, we adopt the following Gamma priors for the parameters 3 :
For each parameter, 1000 posterior samples are drawn after a 200-iteration burn-in period.
The danger of neglecting networks or network dynamics Adopting the "random mixing" assumption about an actually networked population can lead to severe underestimation of the infection rate. Erroneous estimation can also happen if contacts are in fact dynamic but are mistaken as static during inference. Table 2 These results make clear that neglecting the constraints of social networks, even when the quantity of interest is the disease transmission rate, is dangerously misleading, producing estimates that are not at all close to the truth. While incorporating the initial network structure statically throughout the process helps (note that the 95% confidence interval does include the truth), disregarding the evolution of the network over time is still a significant model misspecification and leads to biases. Assessing model flexibility Our proposed model is a generalization of temporal network epidemic processes where the network evolves independently of the contagion (the "decoupled" process), which in turn is a generalization of epidemic processes on static networks (the "static network" process). Therefore, if a set of events is actually generated from the decoupled process, then we expect all the link activation and termination rates to be estimated as the same. Likewise, if the events are generated from the static network process, then we expect all the link rates to be estimated as zero.
To confirm this, experiments are conducted on complete event datasets generated using the two simpler models. Here we only show select results of Bayesian inference on datasets generated from static network epidemic processes ( Figure 4 ) and delegate other results to Supplement S4. We can see that information about a moderate number of events is sufficient to accurately estimate the epidemic-related parameters (β and γ) and uncover the static nature of the network; note how quickly the posterior credible bands for α SI and ω SI shrink toward zero. 
Upon validating the model and inference framework in the previous subsection, we are now ready to assess the performance of our proposed inference scheme in the more realistic setting where epidemic observations are incomplete. In this subsection, we first verify that the MCMC sampling scheme in Section 4 is able to retrieve the parameter values despite the unavailability of exact recovery times in the observed data. Then we compare our DARCI algorithm (Prop. 4.2) with two baselines and show that it produces posterior samples of higher quality and with higher efficiency.
Simulating partially observed data We first generate complete event data using the simulation procedure stated earlier in this section, and then randomly take out η × 100% of the recovery times and treat them as unknown. Meanwhile, a disease status report (as described in Section 4.1) is produced every 7 time units throughout the entire process. Note that the scale of time is insignificant here-the physical length τ of 1 time unit can be thought of as a second, an hour, or a day, since τ only determines how we interpret the parameters (which govern how fast the system evolves per unit time) but not the realization of the process. Therefore, we can regard 1 time unit as a day, and the periodical disease status report as a weekly report.
Efficacy of the inference scheme We demonstrate the efficacy of the inference scheme outlined in Section 4.1 through experiments on an example dataset, where the settings and parameters are the same as those in (18) First set η = 50, that is, randomly select 50% of recovery times to be taken as missing. We can see that for every parameter, the 95% sample credible interval covers the true parameter value, suggesting that the proposed inference scheme is able to estimate parameters from incomplete data reasonably well.
Then set η = 100, which implies that all the recovery times are unknown. Figure 6 presents outcomes of the inference algorithm in this case. Understandably, parameter estimation is affected by the total unavailability of exact recovery times, but the drop in accuracy is marginal. Moreover, the credible bands are slightly wider, suggesting an increased uncertainty with more missingness. 
which equals to 1 when the proposed {r * ,i } i=1:R are consistent with the observed event data in (u , v ] and 0 otherwise. If the proposal is not accepted, then set {r
,i } i=1:R . We label this method by "MH".
Note that, the "MH" method employs the same principle as existing agent-based data augmentation methods (Cauchemez et al., 2006; Hoti et al., 2009; Fintzi et al., 2017 ) that propose candidates of individual disease histories and accept them with probabilities computed through evaluating the likelihood (or an approximation of it) and the proposal density.
But in "MH" the implementation is actually simpler and computationally lighter, because the proposal is conditioned on known infection times and the current posterior draw of γ.
The acceptance step is reduced to inspecting compatibility with known data, thus avoiding the intensive computation of likelihood evaluation.
Although the three methods give the same inference results since they all sample from the same posterior distributions, our data augmentation algorithm (labeled by "DARCI")
is more efficient than the other two in two aspects. First, since we draw a new sample of recovery times from the conditional distribution in (16) in every iteration, there is lower autocorrelation in the resulting Markov chain, which leads to better mixing and fewer iterations needed to achieve a certain effective sample size (especially so when compared with "MH").
Second, the DARCI algorithm is generally much less time consuming when drawing from the conditional distribution in (16) compared to "Reject", because it parses out a configuration of "lower bounds" for imputing missing recovery times while accounting for the constraints imposed by contagion spreading and the dynamics of social links.
Three MCMC samplers are run using the three methods respectively on the dataset showcased above. In each sampler, 1000 consecutive samples are retained for each parameter after a 200-iteration burn-in period. For the resulting chain of each parameter, we calculate the effective sample size (ESS), the Geweke Z-score (Geweke et al., 1991) , and the twosided p-value for the Z-score. The results are presented in Table 3 . Among the three methods, "MH" suffers the most from the correlation between two successive samples, while "DARCI" seems to produce high-quality MCMC samples. Table 3 : MCMC diagnostics for three data augmentation sampling methods, labelled as "DARCI", "Reject", and "MH". "ESS" stands for "effective sample size". The "Z-score" is the test statistic for MCMC convergence proposed by Geweke et al. (1991) , and the two-sided p-value for each standard Z-score is also computed. Samples acquired by MH tend to have higher auto-correlations and thus smaller effective sample sizes. We then compare "DARCI" and "Reject" in their running times (see Table 4 ). A dataset is simulated where there are different numbers of recoveries with unknown times within 5 time intervals. The two sampling methods are applied to draw a set of recovery times for each of those 5 intervals, and over multiple runs, the minimum and median times they take are recorded. Although the two methods draw samples from the same conditional distribution, "DARCI" tends to take less time than "Reject" in one iteration. 
In this section, we apply the proposed model and inference scheme to a real-world dataset on the transmission of influenza-like illnesses among students on a university campus.
The data we analyze in this section were collected in a 10-week network-based epidemiological study, eX-FLU (Aiello et al., 2016) . The study was originally designed to investigate the effect of social intervention on respiratory infection transmission. 590 university students enrolled in the study and were asked to respond to weekly surveys on influenza-like illness (ILI) symptoms and social interactions. 103 individuals further participated in a sub-study in which each study subject was provided a smartphone equipped with an application, iEpi.
The application pairs smartphones with other nearby study devices via Bluetooth, recording individual-level social interactions at five-minute intervals.
The sub-study using iEpi was carried out from January 28, 2013 to April 15, 2013 (from week 2 until after week 10). Between weeks 6 and 7, there was a one-week spring break (March 1 to March 7), during which the volume of recorded social contacts dropped noticeably. In our experiments, we use data collected on the N = 103 sub-study participants from January 28 to April 4 (week 2 to week 10), and treat the two periods before and after the spring break as two separate and independent observation periods (T max = 31 days for period 1 and T max = 28 days for period 2).
We take an infection event as the occurrence of a positive ILI case, and define a recovery event as the end of an individual's flu-like symptoms (with the exact time unknown). A link between two study subjects is activated when the iEpi application registers the initial pairing of their smartphones, and the link is terminated when the Bluetooth detection stops.
Altogether, the data contain 24 infection events in total (the exact times are known), with 14 before the spring break week and 10 after, as well as 45,760 social link activation and termination events. The weekly health status (healthy or ill) of every participant is acquired from the weekly surveys, so we know which individuals recovered during each particular week, but the exact recovery times are unknown. Summary statistics of the data are provided in Table 5 . Overall, infection instance counts peaked in the middle of each observation period and dropped at the end, and the dynamic social network was quite sparse; more activity (in both the epidemic process and network process) was observed in the weeks before the spring break. Further details on data cleaning and pre-processing are provided in Supplement S5.
Since the 103 individuals are sub-sampled from the 590 study participants, which are also sub-sampled from the entire population on the university campus, the real data are actually
observed on an open population. Following the parameterization introduced at the end of Section 3.2, we include the parameter ξ to denote the rate of infection from an external source for each susceptible individual. Every infected individual that came into contact with any infectives within 3 days prior to the onset of symptoms is regarded as an internal case (governed by parameter β); otherwise the infection is labeled as an external case (governed by parameter ξ). This enables the inference procedure stated at the end of Section 4.
The data collected during the two observation periods are considered as separate and independent realizations of the same adaptive network epidemic processes. For each param- No new infection cases took place in week 2, but two participants were already ill at the beginning of the week. The dynamic network remained sparse throughout the duration of the sub-study, except for one instance in week 3-the unusually high network density only occurred on the night of February 4, possibly due to a large-scale on-campus social event.
Week that an external ILI case is expected to occur every other three days. This is a reasonable estimate, since we observed 9 external infection cases within 28 days during the second observation period.
The inferred link rates reflect an interesting pattern in social interactions in this particular population: individuals are reluctant to establish contact and active contacts are broken off quickly-an average pair of healthy people initiate/restart their interaction after waiting 20 days and then end it after spending less than 40 minutes together. Moreover, it seems that on average a healthy-ill link is activated more frequently than a healthy-healthy link, but the former is also terminated faster-this might be because those students who fell ill in the duration of the study happen to be more socially proactive, but once their healthy social contacts realize they are sick and thus potentially infectious, the contact is cut short to avoid disease contraction.
It is also notable that the sample 95% credible intervals for β, α SI and ω SI are relatively wide, indicating a high level of uncertainty in the estimation for these parameters. It is challenging to estimate the internal infection rate β because dataset contains only 6 cases of internal infection in total (5 in period 1, 1 in period 2), providing limited information on the rate of transmission. Similar issues are present for the estimation of α SI and ω SI ; since there were no more than 5 infected individuals at any given time, network events related to them were few and far between. Moreover, since their exact recovery times are unknown, there is additional uncertainty associated with their exact disease statuses when they activated or terminated social links. Such measure of uncertainty, readily available through stochastic modeling and Bayesian inference, provides valuable insights into the amount of information the data contain and the level of confidence we possess when making conclusions and interpretations. The inference outcomes imply that, for example, the real data sufficiently inform the contact patterns among healthy individuals in this population but are limited toward understanding how long a healthy person and a symptomatic person typically maintain their contact.
This paper has focused on enabling inference for partially observed epidemic processes on dynamic and adaptive networks. We formulated a continuous-time Markov process model to describe the epidemic-network interplay and derived its complete data likelihood. This lead to the design of conditional sampling techniques that enable data augmented inference methods to accommodate missingness in individual recovery times.
The proposed framework for dynamic network epidemiology is a generalization and extension of previous methods. Certain concepts, such as infection and recovery rates, disease prevalence, and network density, are analogous to those in existing literature, and their interpretations are relatively straightforward from our generative, mechanistic model. Other concepts, however, are not readily transferable from more traditional settings to the adaptive network setting. One quantity of particular interest in most epidemiological models is the basic reproductive number, denoted by R 0 . It denotes the average number of secondary infections caused by a typical infectious individual at the early stage of the epidemic. Here, "early stage" refers to the initial period of the process during which only a negligible fraction of the entire population is infected, such that any person that comes into contact with an infected individual can be assumed to be susceptible. R 0 is intimately tied with the initial behavior and the final epidemic size (the number of individuals that ever get infected) in SIR models. Loosely speaking, if R 0 < 1, then the epidemic will quickly die down and eventually only a tiny fraction of the population will be affected; if R 0 > 1, then the epidemic will take off and a considerable proportion of the population will be affected by the disease.
For the deterministic SIR model in (1), the basic reproductive number can be written as R 0 = βN/γ, whereas for the stochastic SIR model in (2) and (3), R 0 = βd/γ, whered is the average number of contacts for individuals in the population at the early stage of the epidemic. In temporal network epidemic processes, however, the basic reproductive number does not necessarily carry a clear meaning, nor does it directly relate to the final size of the epidemic (Holme and Masuda, 2015) . For our proposed generative model where the epidemic and network jointly evolve, the definition and expression of R 0 become even more complicated. Here, R 0 may be vaguely expressed as (Tunc et al., 2013; Van Segbroeck et al., 2010) , have investigated 5 two simplified, extreme scenarios where the early stage epidemic behavior can be easily described and R 0 has an explicit expression:
1. Fast network dynamics: when the rate of link changes is significantly faster than that of the infection (and recovery), it can be assumed that the network process has reached a steady state before the disease starts spreading. Then the basic reproductive number is R (F ) 0 = β γ N ψ SI , where ψ SI = α SI /(α SI + ω SI ) is the expected "density" of S-I links at the network steady state.
2. Slow network dynamics: when the rate of link changes is much slower than that of the infection, it can be assumed that the network structure remains the same (as the initial configuration G 0 ) when the epidemic process begins. Thus the basic reproductive
The expression of R 0 is unclear when the network and the epidemic evolve at comparable rates, and the concept "basic reproductive number" may not even be meaningful with an adaptive network, since the system can exhibit complex behaviors (Kiss et al., 2017) with features and outcomes that cannot be described by a single quantity like R 0 . One possible direction of future work may be to investigate, under the coupled process of epidemics and adaptive networks, the general concept and formula of R 0 , or other metrics that quantitatively describe the systematic behavior. 
For an SIS-type infectious disease, the complete data likelihood can be derived following the same steps in Section 3.2. Alternatively, one can slightly modify (14) to arrive at the complete likelihood for an SIS-type contagion. Since an individual doesn't acquire immunity upon recovery, it is equivalent to setting H(t) ≡ S(t) at any time t. Thus the complete data likelihood is
Proof for Theorem 3.1 From (14), we can obtain the log-likelihood:
(β, γ,α,ω|G 0 ) = log L(β, γ,α,ω|G 0 ) = n j=2 logM (t j ) + F j log I p j1 (t j ) + n R log γ + (n E − 1) log β
Taking partial derivatives of the right hand side of (21) with respect to the parameters and setting them to zero yield the results above.
Suppose the observed population is not fully closed, but is a subset of a larger yet unobserved population. Then it is possible for an individual to get infected by an outsider. Let ξ be the "external infection" rate, the rate for any susceptible individual to be infected by any external infectious source, then the complete data likelihood is L(β, ξ, γ,α,ω|G 0 ) = p(epidemic events, network events|β, ξ, γ, α, ω, G 0 )
MLEs for {γ,α,ω} remain unchanged, but estimating β and ξ is less straightforward.
Let (β, ξ, γ,α,ω|G 0 ) be the log-likelihood, then the partial derivatives of the log-likelihood w.r.t. β and ξ are
which do not directly lead to closed-form solutions.
Reparameterizing by ξ = κβ leads to the following partially derivatives
which are slightly more straightforward in form, and can be solved numerically to obtain the MLEs.
If, somehow, we have information on which infection cases are caused by internal sources and which are caused by external sources, then we can directly obtain the MLEs and Bayesian posterior distributions for all the parameters. For an infection event e j (with F j = 1), let Int j = 1 if it is "internal" and let Int j = 0 otherwise. Then the complete data likelihood can be re-written as
where n int E and n ext E are the total numbers of internal and external infection events, respectively.
Estimation for all parameters remains unchanged except for β and ξ. Their MLEs arê
and with Gamma priors β ∼ Ga(a β , b β ) and ξ ∼ Ga(a ξ , b ξ ), their posterior distributions are
When there is missingness in recovery times, the Bayesian inference procedure described in Section 4 can still be carried out, with two slight modifications. First, in the data augmentation step, when drawing missing recovery times in an interval (u, v] , the DARCI algorithm ( Prop. 4.2) inspects I p only for each p ∈ P int , where P int is the group of individuals who get internally infected during (u, v] . Second, in each iteration, parameter values are drawn from the posterior distributions specified in (15) except for β and ξ, for which the posteriors are stated in (27) and (28), respectively.
Supplement for "inference from complete event data" Figure S1 and S2 complement Experiments on larger networks Figure S3 shows MLEs and 95% confidence bands for parameters with complete data generated on a network with N = 500 individuals. Other experimental settings are the same as those in Section 5.1. With a larger population, there tends to be more events available for inference, so the accuracy is in fact improved.
Experiments on different initial network configurations Still set population size N = 100, but instead of a random Erdős-Rényi graph as G 0 , the initial network is a "hubnet": one individual (the "hub") is connected to everyone else in the population while the others form an ER(N − 1, p) random graph, with edge probability p = 0.1. Figure S4 summarizes results of Bayesian inference carried out on complete event data generated in this setting.
Supplement for "Assessing model flexibility" Estimate parameters Θ of the full model on datasets generated from 1) the decoupled temporal network epidemic process with type-independent edge rates, and 2) the static network epidemic process where the network remains unchanged. For both simpler models, fix β = 0.03 and γ = 0.12, and for the former model, let link activation rate α = 0.005 and termination rate ω = 0.05. Still, set population size N = 100 and let the initial network be a random Erdős-Rényi graph with edge probability p = 0.1.
We present, in Figure S5 , the results of Bayesian inference on datasets generated from the decoupled process model. Across four different realizations, it can be observed that, the posterior samples of link activation rates (α SS , α SI , α II ) concentrate around the same mean, and uncertainty is reduced with more events available for inference. Same can be said about the link termination rates, ω SS , ω SI , ω II . This verifies that the proposed model is indeed a generalization of the aforementioned two simpler processes, and the inference method is able to recover the truth under mild model misspecification. 
All infection events and weekly health statuses of all N = 103 individuals are extracted from the weekly surveys. In every survey, study participants were asked if they ever felt ill at all in the past week, if they ever experienced certain symptoms, and, if there were symptoms, when the approximate illness onset time was. We take an "infection" as a positive ILI (influenza-like illness) case, which, following the protocol in Aiello et al. (2016) , is defined as a cough plus at least one of the following symptoms: fever or feverishness, chills, or body aches. We further examine each ILI case and only accept one as a positive infection if the individual also indicated that they "felt ill" in the past week, thus eliminating a small number of reoccurring ILI cases for the same participants 6 . Moreover, since an individual may start exhibiting symptoms at most 3 days after getting infected and becoming infectious, for each infection event, we set the "real" infection time as the reported onset time minus a random "delay time" uniformly sampled between 0 and 3 days.
Social link activation and termination events are obtained from the iEpi Bluetooth contact records. Each time two study devices were paired, the iEpi application recorded the unique identifiers of the devices, a timestamp, and a received signal strength indicator (RSSI). Since
Bluetooth detection can be activated whenever two devices are within a few meters of each other while the two users may not actually be in contact, we only keep those Bluetooth records with relatively strong signals (high values of RSSIs) 7 . If two consecutive Bluetooth records for one pair of devices are no more than 7.5 minutes apart in time 8 , then the two records are considered to belong to one single continuous contact; a social link between two individuals is activated at the time of the first Bluetooth detection record in a series of consecutive records that belong to a single contact, and the link is terminated at a random 6 One particular individual had positive ILI cases and felt ill in week 2, 3, and 5, but not in week 4. We therefore treat his/her illness as an extended one, starting in week 2 and lasting till week 5. 7 The RSSIs range from -109 to 6, and we set the threshold as -90, so only those records with RSSIs larger than -90 are kept. 8 We choose 7.5 minutes as a threshold instead of 5 minutes to accommodate potential lapses in Bluetooth detection.
time point between 1 and 6 minutes after the last Bluetooth detection of a continuous contact.
The resulting processed data contain 24 infection events in total, with 14 before the spring break week and 10 after, as well as 45,760 social link activation and termination events. The weekly disease status (healthy or ill) of every participant can be acquired from the weekly surveys, so we know, for example, if an individual recovered sometime after day 7 and before day 14, but the exact times of all recoveries are unknown.
Instead of assuming the knowledge of which infection cases are internal and which are external, we directly estimate all the parameters based on the likelihood function in (22), solving (23) and (24) for the MLEs of β and ξ.
However, the real data are incomplete, with the exact times of all the recoveries unobserved. We resolve this issue using a naive imputation method-for each recovery, an event time is randomly sampled from a uniform distribution between the time of infection and the earliest time point the individual no longer felt ill (in response to the weekly surveys).
Such imputation, of course, is subject to a considerable level of uncertainty, so we randomly generate 10 differently imputed datasets , obtain the MLEs from every dataset, and then report the averages over the 10 runs (see Table S1 ).
We can see that the MLEs acquired in this manner generally agree with the Bayesian estimates in Section 6.2. 