Modeling regional disease spread over time using a dynamic spatio-temporal model - With an application to porcine epidemic diarrhea virus data in Iowa, US.

Regional surveillance is important for detecting the incursion of new pathogens and informing disease monitoring and control programs. Modeling disease distribution over time can provide insight into the development of more efficient regional surveillance approaches. Herein we propose a Bayesian spatio-temporal model to describe the distribution of porcine epidemic diarrhea virus (PEDV) in Iowa USA. Model parameters are estimated through a Bayesian spatio-temporal model approach which can account for missing values. For illustration, we apply the proposed model to PEDV test results from the Iowa State University Veterinary Diagnostic Laboratory (ISU-VDL). A simulation study carried out to evaluate the model showed that the proposed model captured the pattern of PEDV distribution and its spatio-temporal dependence.

Introduction

The aim of this project was to explore area surveillance methods for livestock farms, with the long-term objective of developing timely and practical surveillance procedures for detecting the introduction of emerging and/or foreign animal diseases (FAD) into a region. Time-to-detection is a critical feature in disease control because delay allows the pathogen to spread and, thereby, increases the difficulty and cost of achieving control and/or elimination. For example, in a study of livestock farms in California, Carpenter et al. (2011) estimated that a 7 day delay in the detection of foot-and-mouth disease virus (FMDV) would result in a median of 13 infected premises holding 677 animals, whereas a 22 day delay would result in a median of 745 infected premises with 6211 animals. However, the imperative of rapid detection must be balanced by the cost of surveillance. An extreme example is the case of bovine spongiform encephalopathy for which the cost of finding one positive case in slaughter cattle in the European Union was reported as # euro 2.3 million for the year 2003 (Heim and Mumford, 2005). Essentially, the challenge is efficient disease detection at an affordable cost.

As reviewed by Saif et al. (2019), porcine epidemic diarrhea virus (PEDV) is an RNA virus in family Coronaviridae causing enteric infection in pigs and capable of producing extensive mortality in PEDV-susceptible neonatal piglets. Clinical outbreaks of PED were first reported in 1971 in England, but PEDV roused little attention until catastrophic outbreaks were reported in China, Thailand, and Korea beginning in 2007. Previously free of the pathogen, PEDV was detected in the U.S. in 2013 (Stevenson et al., 2013) and is estimated to have caused the deaths of 8 million piglets and economic losses of $481 to $929 million (USD) in 2014 (Paarlberg, 2014). Immediately after its detection, efforts to diagnose PEDV infections in commercial swine herds resulted in the submission of massive numbers of samples for testing to the Iowa State University Veterinary Diagnostic Laboratory (ISU-VDL). Testing results, when combined with precise location data, presented a unique opportunity to evaluate spatial surveillance approaches for the detection of a pathogen introduced into a totally naive population using routine diagnostic testing data. Pragmatically, the use of routine diagnostic submissions in on-going surveillance offers both extensive cost savings (sample collection cost is avoided) and the opportunity for real-time detection.

Typically, contagious infectious diseases spread widely following their initial introduction, reach a peak prevalence, and then establish a cycle in the population that mirrors changes in population immunity over time. To study the spread of infectious diseases, several spatio-temporal Bayesian models have been proposed and evaluated in the literature (Denis et al., 2018, Knorr-Held, 2000, Richardson et al., 2006, Lawson, 2013, Watson et al., 2017). In most of these studies, the goal is to identify areas with high or low disease prevalence and to forecast disease prevalence. For example, a Bayesian spatio-temporal conditional autoregressive model was proposed to analyze and forecast the prevalence of antibodies to Borrelia burgdorferi in domestic dogs in the U.S. (Watson et al., 2017). Bayesian methods provide great advantage in the analysis of complex models with complicated data structure. Given the flexibility and generality of the Bayesian framework, we are allowed to cope with complex problems, especially when we have latent variables, missing data, or multilayered probability specifications (Gelman et al., 2013). Notably, unlike maximum likelihood estimation in the frequentist inference, Bayesian methods have the advantage of avoiding integration over the latent variables or missing values, especially given the complexity of the model there is usually no closed form of the integration. The computation is easier if using Markov chain Monte Carlo (MCMC) methods to draw samples from the posterior distribution (Besag and Mondal, 2013, West and Harrison, 2006, Chap.15). This is relevant to the present study because we propose to use a dynamic model with latent variables and about 27.5% of PEDV test result data from the ISU-VDL are missing.

Iowa is a highly agricultural state located in the north-central U.S. The state is 150,930 square kilometers in size and divided into 99 counties of relatively uniform size (Fig. 1). Federal Orders (June 5, 2014, January 4, 2016) entitled “Reporting, Herd Monitoring and Management of Novel Swine Enteric Coronavirus Diseases” required that swine producers report cases of swine coronavirus infections and provide premises identification number, date of sample collection, type of unit sampled, the diagnostic test used, and the results of testing. This provided an extensive dataset on PEDV in Iowa swine farms. The initial dataset included a total of 30,843 PEDV PCR test results from samples collected from premises located in the state of Iowa and tested at the ISU-VDL between May 2014 and March 2017. A data query was performed to retrieve all PEDV test results. In order to maintain data integrity, the results were processed to remove any result that did not have an associated valid address representing a swine premises. Premises identification numbers (PIN) and premises submission level identifiers were entered into Google Earth and if a swine premises was located at the entered address the result was maintained. If the address or PIN provided did not pertain to a swine premises, the result was removed from the data set. Binary (positive/non-positive) results were used for summary and analyses. In our study, we selected the 6-month period, August 2016 to January 2017, with the most complete data to address the spread of PEDV.

Fig. 1

Monthly observed proportion of PEDV among counties in Iowa from August 2016 to January 2017 in the Veterinary Diagnostic Laboratory at Iowa State University (ISU-VDL) testing result.

In this paper, we propose a Bayesian spatio-temporal model for the distribution of PEDV at the county level in Iowa that accounts for the effect of spatial distance among counties. As illustrated in Fig. 1, the observed proportion of PEDV in each county from August 2016 to January 2017 shows a clear spatio-temporal dependence. There are a few counties with no tested samples at certain time points. The top two reasons that can address the majority of missing are farmers may not submit samples for testing if the disease status is not a concern and the farm lost all of the animals due to acute outbreak. There might be a small portion of missing due to reasons such as producers may have elected to submit samples for testing to a laboratory other than the ISU-VDL. Our goal is to build a model reflects the major trend of missing. Thus to account for the missing values, we added constraints based on the two major reasons described earlier to the Bayesian spatio-temporal model.

The paper is organized as follows: Section 2 presents the proposed spatio-temporal model and the corresponding Bayesian inference including MCMC algorithm and choice of prior. Section 3 shows the simulation results. Section 4 reports the results of PEDV data analysis. Section 5 provides the discussion of the results and evaluation of the model performance. Section 6 concludes the paper with a summary and some future works.

Modeling procedure

Bayesian spatio-temporal model for PEDV distribution

In this section, we propose a Bayesian spatio-temporal model to describe changes in the distribution of PEDV over time. Suppose there are counties in our study, and the time ranges from 0 to with each time point representing one month. Let be the total number of samples submitted for testing from county at time , and is the number of samples tested positive, and . follows a binomial distribution:where is the probability of samples in county at time being tested PEDV positive, and is its logit transformation. Here are latent variables which have spatio-temporal dependence. That is, is influenced by its previous time value and the status of neighbor counties of at time . To account for this spatio-temporal dependence among counties, we model the increment of , , to be a linear function of the probabilities of being tested as PEDV positive from other counties. In addition, the magnitude of the impact of neighbor counties on county is also related to their distance. Therefore, we propose the model in the following way:Here is the intercept which accounts for the average increment of across time without considering any spatial interactions. is the distance between county and county , and is the probability of samples from county at time being tested PEDV positive. In this model, the spatial impact of county on county is quantified as , which decreases exponentially as the distance increases. Thus for any county, its neighbor counties have larger influence than other counties that are far away. The parameter controls the overall magnitude of the spatial impacts on from all other counties at the previous time point. The structure of the model is autoregressive and thus has the Markov property. There are also some other covariates can potentially affect , such as the number of farms in the county and the transportation infrastructure (interstate, US and state highways) in the county. The latter was considered because pigs are routinely moved between farms and counties in the production cycle. The presence of more farms and more developed transportation infrastructure may have affected the distribution of the virus. Thus, we include covariates in the model:for and .

Here is the dimensional vector of constant one and the covariates that affect the distribution of PEDV. is the corresponding vector of regression parameters. 's are the random errors accounting for other unmeasured factors, which are independently identically distributed as . The initial states, , for , are modeled as

Considering missing values

Within the considered time period, there are some counties that have no samples submitted from farms for testing in at least one month, i.e., missing values. Presumably, farms do not submit samples for testing if there are no clinical signs of disease among piglets. That is, we assume that farms submit samples for testing only if the proportion of clinically-affected piglets is relatively high. Then one assumption we can make in the analysis is that the decision not to submit samples means that the disease status at that time point is not a concern, thus a low proportion (less than 0.1) of PEDV can be assumed. Therefore, if a county at time has missing, we put constraints on as to make the proportion of PEDV to be less than 0.1.

In addition, if a farm has submitted samples at time but none at time , i.e. is observed but is missing, it is very likely that all the piglets in that farm died. PEDV infection remains endemic, and it is impossible for the infection to disappear naturally within a short period of time. In many cases during the acute outbreak, farms were unable to produce healthy piglets for weeks or even a few months. These cases do not fit in our model, since the model is only intended to describe the spatio-temporal trends in PEDV. Therefore, the spatio-temporal dependency will be cut for such cases, which results the corresponding to be no longer influenced by others and we propose a noninformative distribution on for , i.e. .

To incorporate these constraints into the model, we first define the following five sets:

For and all county , we have which is the same as the model described in Section 2.1. For and county , we have the following model:

For , this model sets density of as a constant, which means its distribution is noninformative. For , the model is the same as the Bayesian spatio-temporal model shown in Eq. (1).

Joint probability density function (pdf) of the PEDV outcome

Here we develop the joint pdf of the logit transformed probability of tested positive , where represents a vector of logit transformed probability of tested positive at time for all counties. The joint pdf of can be derived based on the spatio-temporal model described above. Given that there are two different model settings, one is without missing value described in Section 2.1 and the other is under missing values with constraints described in Section 2.2, we have two different joint pdf forms.

Joint pdf under complete data setting

For , given 's are independently and identically distributed, has the density as:

For , based on the Markov property, the joint density of has the following form

Thus the joint pdf for is derived as:

Let be the observed outcomes for all counties at all time, where each is a vector of observed outcome for all counties at time point .

Finally, the joint density of and has the following form:with and as shown in equation (3).

Joint pdf under missing value with constraints setting

The joint pdf of is then derived as:where is the number of elements in set and is an indicator function.

Finally, the joint density of and has the following form:with .

Bayesian inference

Bayesian posterior inference will be used to estimate the parameters of interest include , , , and . The prior distribution for each of the parameters are as follows:

The parameterization of the Gamma distribution prior for has the mean of and variance , and has the mean of and variance . The conjugate non-informative priors are used for all model parameters, where , and take the value 100, and . Sensitivity analysis has been performed for different choices of prior parameters (, or 0.0001). The results showed that different values of prior distribution parameters only have little impact on the estimates. Given the joint pdf shown in Eqs. (4) and (6), the joint posterior distributions can be derived as:

Markov chain Monte Carlo (MCMC) with a combination of Gibbs and Metropolis-Hastings (MH) steps is used to draw samples from the posterior distributions of the parameters sequentially. Based on the draws from the posterior distributions we can make inferences for the model parameters. The detailed posterior sampling algorithm is described in section Appendix: MCMC Algorithm, and the computing program is written in R.

The parameters , , , and have closed form of full conditional posterior distributions from which we can directly sample. However, the conditional posterior distribution of given other parameters does not have a closed form. We use random walk MH to update one by one, which means the proposal distribution is . Here, how to scale is of importance: too small values will make the chain move very slowly; too large values will result in very low acceptance rate of the proposals. To choose in an appropriate way, we use the idea of adaptive MCMC in Haario et al. (2001). At the beginning, we use pre-specified such as , to carry out a MCMC simulation. After an initial period, for each , a good choice of would be the standard derivation of the posterior samples of in this MCMC. After that, a new MCMC is implemented with these new 's.

Simulation study

A simulation study was performed to evaluate the performance of both models (complete and missing values). Under each model, there are two different settings. In both settings we set the number of counties and the total period of time to be the same as the real data. In the first setting, we have the number of samples tested from each county at each time, , fixed as 25 for all counties at all times, since the averaged number of samples tested in the PEDV test data over all counties is about 25. In the second setting, we fix as 100 for all counties at all times. Posterior means of parameters are used as the true parameters in the simulation. Each setting has 100 simulations. For each simulation, the data is generated from the true parameters, and then we applied the proposed models to estimate the parameters. The parameter estimations from the 100 simulations are compared with the true parameters to evaluate the models performance.

In the simulation under complete data model, we first generate , for , from , where and . Then for , the is generated as equation (1), where model parameters are set to be the posterior means from real data analysis, i.e. , , and . Given the generated , the number of samples tested positive, , can then be generated from a binomial distribution given size of or 100 and probability of . The result based on 100 simulation for model with complete values is shown in Table 1.

Table 1

Average and standard error of posterior means for simulation studies of the model with complete data (based on 100 simulations).

N_s,t	Parameter	True value	Averaged posterior mean	SE of posterior mean
25	β_0	−0.0942	−0.0933	0.1490
	α	0.9705	0.9589	0.0907
	μ	−3.9164	−3.9675	0.3851
	σ^2	2.0935	2.3269	1.0320
	τ^2	2.1679	2.3425	0.4663
100	β_0	−0.0942	−0.0907	0.1299
	α	0.9705	0.9690	0.0710
	μ	−3.9164	−3.9858	0.2980
	σ^2	2.0935	2.2916	0.8586
	τ^2	2.1679	2.2044	0.3194

In the simulation under missing data model, we first generate , for , from , where and . Then for , is generated as equation (1), where model parameters are set to be the posterior means from real data analysis, i.e. , , and . Given that we assume the proportion of PEDV for samples that are missing in the PEDV data are less than 0.1 for all and . We calculate the proportion of missing among samples with proportion of PEDV less than 0.1 in the PEDV data as . In the simulated data for 's with values of , we randomly set of them as missing. The simulation result for model with missing values constraints is shown in Table 2.

Table 2

Average and standard error of posterior means for simulation studies of the model with missing values. (based on 100 simulations).

N_s,t	Parameter	True value	Averaged posterior mean	SE of posterior mean
25	β_0	−0.7324	−0.6220	0.1433
	α	0.9188	0.9065	0.0176
	μ	−3.7945	−3.7387	0.1665
	σ^2	0.9721	0.7542	0.3117
	τ^2	3.0639	3.1134	0.4806
100	β_0	−0.7324	−0.6425	0.1248
	α	0.9188	0.9073	0.0161
	μ	−3.7945	−3.8084	0.1345
	σ^2	0.9721	0.8769	0.2261
	τ^2	3.0639	3.1577	0.3665

Results of the PEDV data

The proposed spatio-temporal model is applied to the PEDV test results from the ISU-VDL with testing dates from August 2016 to January 2017. The testing results are at the county level to protect client confidentiality. We are interested in modeling the monthly spatio-temporal pattern of PEDV. In Fig. 1, we can see that the averaged observed proportion of PEDV over all counties is very low, with an average of 0.023. Over time, the proportion of PEDV continues to increase until the end of the sixth month. Samples for testing were received monthly from 85 out of 99 counties during this 6-month period. We proposed two models, one is for the complete data setting which will be applied in Section 4.1, and the other is for the incomplete data setting and will be applied in Section 4.2.

Results under complete data setting

After deleting the counties with missing values to obtain a complete dataset, only 31 counties are left. The number of samples submitted for testing from each county at each month ranges from 1 to 331. We have considered six covariates in the model in order to achieve a more accurate estimation of the spatio-temporal PEDV trend. These six covariates are season, number of interstate highways, US highways, and state highways in each county, area of each county in square miles and the number of farms in each county divided by one thousand. However, none of these six covariates has significant contribution to the model (Table 3). Thus, our final model only includes the intercept, .

Table 3

Posterior results for complete data model with all six covariates analyzed with PEDV data.

	Posterior mean	Posterior SD	LB 95% CI	UB 95% CI
β_0	−0.7775	0.9739	−2.6839	1.1574
β_I−HW	−0.0036	0.2970	−0.5853	0.5800
β_US−HW	−0.0643	0.2191	−0.4938	0.3640
β_ST−HW	−0.0390	0.0994	−0.2348	0.1534
β_area	0.0013	0.0020	−0.0027	0.0053
β_farmnum	1.1584	1.1540	−1.0956	3.4993
β_season	−0.0491	0.4094	−0.8427	0.7807
α	1.0039	0.2859	0.4458	1.5909
μ	−4.1299	0.5359	−5.2853	−3.1673
σ^2	2.4397	0.6305	1.4174	3.8664
τ^2	2.4416	1.3383	0.7108	5.8476

To start the MCMC iteration, a common choice of the initial value of is the logit of the observed proportion of PEDV with . However, given the number of samples tested in some counties are small (less than 3), it is easy to have a taking value of 0 or , which makes the corresponding to be either 0 or 1. In such case, is no longer valid. To avoid setting as an arbitrary number such as 0.001 or 0.999, we borrowed the idea of Agresti and Coull (1998) and used the following weighted average to approximate :where is the quantile of the standard normal. Here we choose and then . This makes fall into and shrinks the observed proportion of PEDV to 0.5. When goes up, the weight assigned to increases and the shrinkage becomes less. We implement this correction for all 's, and then use as the initial value of .

The estimation of parameters were based on 15,000 MCMC iterations after disregarding the first 10,000 burn-in period. The Gelman-Rubin diagnostic test was performed to check the convergence, and we ran three different chains with starting values as 50, 0, and 50. The scale reduction factor from the Gelman–Rubin diagnostic test is 1 means that between variance and within chain variances are equal, which is an indication of convergence. The multivariate effective sample sizes is 3581 for all parameters of interest, based on 15,000 MCMC with 10,000 burn-in. Table 4 shows the posterior means, posterior standard deviations and 95% Bayesian credible intervals (CI) from the model without missing values.

Table 4

Posterior results for complete data model without covariates analyzed with PEDV data.

	Posterior mean	Posterior SD	LB 95% CI	UB 95% CI
β_0	−0.0942	0.1542	−0.3959	0.2071
α	0.9705	0.2556	0.4748	1.4862
μ	−3.9164	0.4285	−4.7596	−3.1277
σ^2	2.0935	1.0108	0.7055	4.5333
τ^2	2.1679	0.5347	1.2651	3.3686

Model under missing value with constraints setting

The spatio-temporal model under missing value with constraints setting is applied to the same PEDV dataset. We are keeping all testing results from those 85 counties regardless if they have any values of missing or not during the 6 months period. Notice that there are 99 counties in total in Iowa, but farms tend to be concentrated in certain counties (farms are not evenly distributed across counties) and not all farms will have had clinical PEDV during this six-month period. Among the 85 counties over the six-month period, there are about 27.5% values missing. This missing rate is calculated as dividing the total number of samples that are observed by . Fig. 2 shows the observed proportion of PEDV for the 85 counties that with at least one observation over 6 months, where each line represents a time series for one county. Again, similar as under the complete data setting, we first considered six covariates in the model but none of these covariates have significant contribution to the model (Table 5). The final model only includes the intercept with result shown in Table 6.

Fig. 2

Plot of observed proportion of PEDV over 6 months with missing values.

Table 5

Posterior results for missing values model with all six covariates analyzed with PEDV data.

	Posterior mean	Posterior SD	LB 95% CI	UB 95% CI
β_0	−1.4744	0.5956	−2.6594	−0.3215
β_I−HW	0.0172	0.1804	−0.3342	0.3790
β_US−HW	−0.1283	0.1505	−0.4328	0.1610
β_ST−HW	−0.0286	0.0685	−0.1626	0.1052
β_area	0.0019	0.0011	−0.0004	0.0042
β_farmnum	0.8997	0.8757	−0.8077	2.6383
β_season	−0.4391	0.3368	−1.1014	0.2153
α	0.9506	0.0856	0.7788	1.1141
μ	−3.9780	0.2683	−4.5710	−3.5207
σ^2	1.1178	0.3958	0.5244	2.0535
τ^2	2.9549	0.5472	2.0550	4.2167

Table 6

Posterior results for missing value model without covariates analyzed with PEDV data

	Posterior mean	Posterior SD	LB 95% CI	UB 95% CI
β_0	−0.7324	0.1338	−1.0068	−0.4806
α	0.9188	0.0886	0.7493	1.0951
μ	−3.7945	0.2272	−4.2546	−3.3776
σ^2	0.9721	0.3640	0.4433	1.8664
τ^2	3.0639	0.5012	2.2247	4.1546

Discussion

The result of the simulation study with complete values is shown in Table 1, and the simulation result for model with missing values constraints is shown in Table 2. We can see for each parameter the average of the 100 posterior means from the 100 simulations is close to the true value, compared with the standard error of these 100 posterior means, for both models. Both simulation results support the validity of the proposed models under complete data setting and missing data settings, respectively.

For the PEDV data analysis result for complete data model with only the intercept, Table 4 shows the posterior means, posterior standard deviations and 95% Bayesian credible intervals (CI) from the model without missing values. Except the intercept all of the parameters are significant, given their 95% credible intervals do not cover 0. indicates that the proportion of PEDV of a county is influenced by other county , where . If the distance between two counties is larger, the influence on will be smaller. Also, if the probability of tested positive, , of county is smaller, its influence on will be smaller. The result for missing data model with only the intercept is shown in Table 6. The parameter estimates are very similar to the model with complete data except that the estimates for is smaller and has significant contribution to the model (with 95% Bayesian credible interval doesn’t cover 0). This means that without considering any spatial interactions, the average increment of the proportion of PEDV across time under the missing value model is larger than under the complete data model. Also notice that the posterior standard deviation for the parameters for model with missing values (Table 6) are all smaller than those without missing values (Table 4). Together with the consistent posterior mean estimates between the two models, it is shown that the model with missing values provide more precise estimates by incorporating more data. Thus when there are missing values, the analysis with missing value data model is preferred to the analyses with complete data only.

Conclusions and future work

In this project, we used a spatio-temporal model to account for the changes in the spatio-temporal pattern of PEDV. The PEDV test result data we used were derived from routine diagnostic testing, which presents a huge advantage in terms of lower cost and timeliness (real-time surveillance). Applying the proposed model to PEDV test results from the ISU-VDL, the parameters estimation results showed strong evidence of distance-dependency at a regional (county) level. Given the importance of efficient regional surveillance, the proposed model can be used to predict the changes in the distribution of PEDV over time and thus lead to better guidelines for detecting the incursion of new pathogens and informing disease monitoring and control programs. For future work, the proposed model could inform the development of sampling guidelines that could be more timely and sensitive in detecting emerging diseases than surveillance based on simple random sampling or passive methods.