Evaluating the Impact of Environmental Temperature on Global Highly Pathogenic Avian Influenza (HPAI) H5N1 Outbreaks in Domestic Poultry

The emergence and spread of highly pathogenic avian influenza (HPAI) A virus subtype H5N1 in Asia, Europe and Africa has had an enormously socioeconomic impact and presents an important threat to human health because of its efficient animal-to-human transmission. Many factors contribute to the occurrence and transmission of HPAI H5N1 virus, but the role of environmental temperature remains poorly understood. Based on an approach of integrating a Bayesian Cox proportional hazards model and a Besag-York-Mollié (BYM) model, we examined the specific impact of environmental temperature on HPAI H5N1 outbreaks in domestic poultry around the globe during the period from 1 December 2003 to 31 December 2009. The results showed that higher environmental temperature was a significant risk factor for earlier occurrence of HPAI H5N1 outbreaks in domestic poultry, especially for a temperature of 25 °C. Its impact varied with epidemic waves (EWs), and the magnitude of the impact tended to increase over EWs.


The Integrated Approach of Bayesian Cox Proportional Hazards Model and Besag-York-Mollié (BYM) Model
Clayton [1] formulated the Cox proportional hazards model using the counting process introduced by Andersen and Gill [2] and also discussed the approach of Bayesian estimation using Markov Chain Monte Carlo (MCMC) method. Say we have n HPAI H5N1 outbreaks. Let N i (t) (i = 1,..., n) count the number of outbreaks which occurred up to time t. Then N i (t) can be regarded as a counting process with random intensity process [1,2]: where, is the familiar Cox proportional hazards model: represents the underlying baseline hazard function, is the unknown regression coefficients, and is the observed or measured covariates for the ith outbreak; is a process taking the value 1 or 0 according to whether or not outbreak i is observed at time t.
For HPAI H5N1 outbreaks, we have data D= N i (t), Y i (t), environmental temperature , unknown regression parameter β and: (2) Then the joint posterior distribution of Bayesian estimation for model (1) is: Assuming non-informative censoring, the likelihood of the data ( ) is proportional to , where: This is as if the counting process increments in the time interval [t, t + dt) are independent Poisson random variables with means : Following formula (1), can written as: where, a vague Normal prior distribution was used for regression coefficient  (~N(0, 0.0001)): is the increment of the integrated baseline hazard function occurring during the time interval [t, t + dt).
Since the conjugate prior for the Poisson mean is the Gamma distribution, it would be convenient if is a process in which the increment is distributed as Gamma distribution. So, the increment was assigned the vague Gamma prior distribution as suggested by [3], namely: where, is a prior guess of ; c represents the degree of confidence in that guess and small values of c correspond to weak prior beliefs. In this study, we set (9) (r is a guess at the failure rate per unit time; dt is the size of the time interval). The vague Gamma prior distribution (G (0.01, 0.01)) was used for both c and r. In our study, there are many factors that can affect HPAI H5N1 outbreaks, but only environmental temperature was obtained. To efficiently explore the impact of environmental temperature on HPAI H5N1 outbreaks, the effects from the other unobserved and unknown variables that are related with HPAI H5N1 outbreaks should be controlled. However, taking into account all the risk factors in a practical analysis is always impossible either because of the difficulties of obtaining data or our limited knowledge on disease occurrence and transmission. Hence, we adopted the idea of the Besag-York-Mollié (BYM) model to represent all the unmeasured and unknown risk factors using two latent variables.
The BYM model was first introduced by Clayton and Kaldor [4] and further developed by Besag et al. [5] It decomposed the position-specific random effects of all the unmeasured and unknown variables into two components: one component is the clustering or correlated heterogeneity, which represents the variables that, if observed, would display substantial spatial structure in that the values for a pair of contiguous zones would be generally much more alike than for two separate zones (e.g., precipitation and elevation) and another component is the uncorrelated heterogeneity that models those intra-area effects that vary in an unstructured way between areas (e.g., gender, age, control strategy). That is, the total effects can be divided as: where, represents the total effects; is the impact from the baseline level; is the measured covariates that are known or suspected to be related with the disease (e.g., environmental temperature in this study); and are the surrogates (latent variables) for the clustered and uncorrelated heterogeneities from all the unknown and unobserved covariates, respectively.
For the uncorrelated heterogeneity, the independent normal prior distribution was used: (11).
For the correlated component , a spatial correlation structure should be used, where the risk estimation in any area depended on the neighboring areas. The conditional autoregressive (CAR) model proposed by Besag et al. [5] was adopted Lawson et al. [6]: where: if i and j outbreaks are adjacent, and 0 if they are not.
The vague Gamma prior distribution (G (0.01, 0.01)) was used for both and , as suggested by Bernardinelli et al. [7]. By combining the above two ideas of Bayesian Cox proportional hazards model and BYM model, the integrated approach to analyze the impacts of environmental temperature on HPAI H5N1 outbreaks was defined as: The meanings of the parameters ( , , , , , , , ) and corresponding prior specifications ( , , , , ) are the same as the above mentioned notations.
WinBUGS version 1.4 (Imperial College and MRC, London, UK) was used to fit the integrated model. One chain Markov Chain Monte Carlo simulation was used for parameter estimation. Model convergence was assessed by visual examination of time series plots of samples, and the inference of the parameters was based on a burn-in of 10,000 iterations after 30,000 iterations.