‘flré\OV—ln§1'. '. «‘9; fl <0. a. j... . k .3. f .13... LT: (xx.- . ravmuwflmfi... ” .. ,....w.% .wmfiw... . 11350941. r was... % wk fir . . .110». .mv! 31!. :3. $51.11} 3.9...bs . . , 3.24.. . .5: h.” 5., 3.1K . P- A 1 _ I. II ....1....xua.....aru.. any...» 1m." a...“ .3. uréais. .3 3c .. .. 5.333 .. ... )- ..\ ‘5'. «or... .ufi‘ 5.. u. k. :uijawflntafi fiwfliifl .5. .flhfifi. $.31" .ir. a: I Each . .i {33.7 i. I. .113)... w £1535 .. :75? i a... a. .......93.. v RI........P¢..?L .52; 2.2, 1 tutti"; II it}?! .25. ‘ l. 3.3191). .Ia.‘ .12.}!1. . 23...}: .9135! ll‘ . r . 2A9 .... :, ... 3!. zt‘aglvtsiuihflna. .f i 5...: ‘ Jai. 3.83:2.14! .3): is...‘ . tari7‘ltrfl ‘14.... affix... Amhivfluwmhuflfim | . Sara; 3 LIBRARIES ‘ V MICHIGAN STATE UNIVERSITY - > 5. EAST LANSING, MICH 48824-1048 935.1'75255 This is to certify that the dissertation entitled SPACE-TIME MODELING AND APPLICATION TO EMERGING INFECTIOUS DISEASES presented by CHENG-YU LEE has been accepted towards fulfillment of the requirements for the Ph.D. degree in Forestry and EEBB program —.-.-.-.-.-.-.-.-.-._.-.-.- —— 4 - - fijor ProtessE/Jr’é Signature Mflt/i /0 .7005 Date MSU is an Affirmative Action/Equal Opportunity Institution PLACE IN RETURN BOX to remove this checkout from your record. TO AVOID FINES return on or before date due. MAY BE RECALLED with earlier due date if requested. DATE DUE DATE DUE DATE DUE 2/05 C‘EIRC/Dateouejndd-pJS SPACE-TIME MODELING AND APPLICATION TO EMERGING INFECTIOUS DISEASES By Cheng—Yu Lee A DISSERTATION Submitted to Michigan State University in partial fulfillment of the requirements for the degree of DOCTOR OF PHILOSOPHY Department of Forestry; Ecology, Evolutionary Biology, and Behavior Program 2005 ABSTRACT SPACE-TIME MODELING AND APPLICATION TO EMERGING INFECTIOUS DISEASES By Cheng-Yu Lee A refined modeling framework for space-time analyses, specifically developed for general Space-Time Autoregressive Moving Average (STARMA) models, is proposed. To enhance accuracy and performance of these analyses, statistical tools and algo- rithms were extended from univariate case to space-time case, including space-time extensions of the Hannan—Rissanen algorithm, the bias-corrected Akaike information criterion, and the Bayesian information criterion. Methods for assessing statistical significance of model parameters are also presented. A general-purpose statistical software, called Integrated Enviromuent for Analyzing STARMA models (IEAST), is developed for space-time analyses in this research. As an empirical example. the framework and these space-time modeling methods are then applied to investigate the spreading dynamics of West Nile virus (WNV) epidemic in crews and humans in the Detroit Metro area in 2002. Both datasets of dead crews and human cases fit very closely to those expected from a purely STAR (Space-Time Autoregressive) pro- cess having low spatial and temporal orders. The use of the STARMA model allows estimation of the rate of spread of W NV at different spatial scales and thus charac- terization of the expected spatial and temporal scales. In addition, a space-time cross correlation analysis between crew and human cases is conducted. The result shows that there exists high cross correlation between dead crow and human cases at spe- cific spatial and temporal lags. This evidence provides a foundation for the control of human VVNV epidemics by using dead crows as a sensitive indicator variable. Statis- tical inferences from a biological point of view based on these analyses can be used to formulate the prevention and control policies for WNV. The determination of spatial- temporal autoregressive parameters using STARMA holds considerable promise for characterizing emerging infectious diseases. To my parents and my wife iv ACKNOWLEDGEMENTS First, I would like to thank my dissertation adviser Professor Bryan K. Epperson for his enthusiasm and support during the years. It would not have been possible to complete this dissertation without Professor Epperson‘s valuable. help. Further, Professors Joseph P. Messina, Kim T. Scribner. and Edward D. Walker are greatly acknowledged for their guidance and inspirational cooperation during my study. I would also like to acknowledge data collation efforts of Dr. Erik Foster and Amy Williams. Dr. Jonathan S. Patterson (Diagnostic Center for Population and Ani- mal Health, l\-*Iichigan State University) performed iInmunohistoclremistry analyses to confirm West Nile virus infection on a subsample of dead crows. This study was supported by the Michigan Agricultural Experiment Station and the Center for Emerging Infectious Diseases, Michigan State University; and a coop- erative agreement with the Centers for Disease Control and Prevention, US Public Health Services. V’ TABLE OF CONTENTS LIST OF FIGURES viii LIST OF TABLES x 1 Introduction 1 1.1 Space-Time Processes and Space-Time Dependency in Biology and Ecology .................................. 1 1.2 Statistical Modeling and Various Space-Time Models ......... 3 1.3 Scope of the Work ............................ 6 1.4 Organization of the Thesis ........................ 8 Space-Time ARMA Models 11 2.1 Historical Background .......................... 11 2.2 Model Definitions ............................. 13 2.3 Spatial Correlation Structure and Weight Matrices .......... 16 2.4 Box-Jenkins Modeling Approach ..................... 17 2.5 Space-time Autocorrelation and Partial Autocorrelation Functions . . 18 2.6 Assumptions ................................ 19 2.7 Finest Lag Structure and Selection of Spatial Correlation Structure . 23 Development of the Proposed Method for STARMA Modeling, and IEAST 28 3.1 Specification of Spatial Correlation Structures and Derivation of Weight Matrices .............................. 29 3.2 Identification Stage ............................ 35 3.3 Estimation Stage ............................. 36 3.4 Diagnostic Checking Stage ........................ 41 3.5 A Refined STARMA Modeling Procedure ................ 46 3.6 Cross Correlations and Space-Time Regression Models ........ 47 3.7 The Development of IEAST ....................... 49 Application to the Spreading of West Nile Virus in Detroit Metro Area 61 4.1 West Nile Virus .............................. 61 vi 4.2 Description of the Datasets and Diagnoses ............... 65 4.3 Assumptions and the Specification of Spatial Correlation Structure . 67 4.4 Space-Time Trends and De-Trending .................. 68 4.5 Results for the Analysis and Modeling of the Dead Crow Dataset . . 71 4.6 Results for the Analysis and Modeling of the Human Case Dataset . . 78 4.7 Summary of Analyses for the Dead Crow Dataset and the Human Case Dataset .................................. 83 4.8 Results for the Cross Analysis Between Human Case and Dead Crow Datasets .................................. 85 4.9 Discussion ................................. 89 5 Conclusion 107 BIBLIOGRAPHY 1 10 vii 1.1 2.1 2.2 3.1 3.2 3.3 3.4 3.5 3.6 3.7 4.1 4.2 4.3 4.4 LIST OF FIGURES Relationship between data-generating and statistical modeling proce- dures .................................... An example of the definition of spatial orders / lags in systems with two spatial dimensions. ............................ Box-Jenkins modeling approach ...................... STACF and STPACF for a simulated pure STAR dataset with maxi- mum temporal order 2 and maximum spatial order I: Zt = 0.5Zt11 + 0.3W(1)Zt_1 + 0.15Z¢_2 + E; ........................ STACF and STPACF for a pure STMA dataset with maximum tem- poral order 1 and maximum spatial order 1: Z = Et — (—0.6)et_1 — (——0.4)W(1)et_1 ............................... The refined STARMA modeling procedure. The area enclosed by the dotted line is for the case where model type cannot be decided from the STACF/STPACF. AICC / BIC of every possible model type is eval- uated and compared, in order to select the best one ........... The linear relationship between space-time variable Z and independent variable Y under the assumption of space-time regression (STR) model. . Inputs and outputs of IEAST ....................... System structure of the software IEAST. The gray blocks are the com- ponents of IEAST which consist of three major parts: menu—driven interface, interpreter, and core functions. ................ The first two levels of the menu system in IEAST ........... Dead crows reported in the Detroit Metro area. The rectangular area enclosed by thick line is the area retained after truncation. ...... Dead crow dataset in the Detroit Metro area with 10 x 10 cells in space and spanning 28 weeks in time. ..................... Human cases reported in the Detroit Metro area. The rectangular area enclosed by thick line is the area retained after truncation. ...... Human cases in the Detroit Metro area with 10x10 cells in space and spanning 13 weeks in time ......................... viii 10 26 27 CH CH 1;“ ‘1 g“! 00 59 60 92 93 94 4.5 4.6 4.7 4.8 4.9 4.10 4.11 4.12 4.13 4.14 4.15 The spatial trends (averaging over time) of dead crow (upper) and human case data. The areas cover by these two datasets are not aligned geographically ................................ The spatial and temporal trends removed from the dead crow dataset. Space-time autocorrelation analysis of the de-trended dead crow dataset. Space-time autocorrelation analysis for the modeling residuals of dead crow dataset ................................ 96 97 98 99 The spatial and temporal trends removed from the human case dataset. 100 Space-time autocorrelation analysis for the de—trended human case dataset. .................................. Space-time autocorrelation analysis for the modeling residuals of hu- man case dataset ............................. Dead crows reported and area for cross analysis ............. Human cases reported and area for cross analysis. ........... Space-time cross correlation between the human case data and the dead crow data ............................... Space-time cross correlation between the filtered human case. data. and the filtered dead crow data. ....................... ix 195 3.1 4.1 4.2 4.3 4.4 4.5 LIST OF TABLES Aecuracies of model type selection using variances, AICC. BIC. and -AICC*BIC based on 150 Monte Carlo simulated datasets. ...... The estimates and significance levels of the model parameters (STAR(MaxT=3, MaxS=4)). ...................... The maximum spatial lags and the corresponding distances between cell midpoints for that lag, for three models using different cell numbers (larger numbers, smaller areas). ..................... The maximum spatial lags for positive autoregressive effects and the corresponding distances between cell midpoints for that lag, for three models using different cell numbers .................... Initial estimates and significance levels of the model parameters (STAR(MaxT=6, MaxS=6)) for the human dataset ........... Final estimates and significance levels of the model parameters (STAR(MaxT=4, MaxS=6)) for the human dataset ........... 76 80 81 CHAPTER 1 Introduction 1.1 Space-Time Processes and Space-Time Depen- dency in Biology and Ecology Spatio—temporal models have gained widespread popularity for the last decade. One major reason for this is an abundance of new challenging applications arising in the environmental sciences and epidemiology. Typical examples include forecasting of global climate change, infectious disease mapping, and their inter-relationship. Space- time datasets are usually very large and, therefore, require substantial computing power for modeling. The major improvement in computational power, especially personal computing, is another sigI'Iificant cause for the recent surge in using the models. Biological and ecological data are often organized by units of time as well as by geographic locations. The procesSes that produce, such data may have strong corre- lations not only in time but also in space. It is not always reasonable to study such processes by considering space and time separately. With increasing accessibility and accuracy in remote sensing technology, large scale analyses (especially in space) of space-time data become possible. This trend highlights the importance and necessity of space-time analysis in various disciplines. In particular, some problems with high dependency in space and time, such as precipitation / temperature forecasting, pollu- tion spreading, population succession of organisms. and epidemics, require space-time analysis and modeling to reveal or even forecast process dynamics. In the natural world, biological or ecological systems are neither spatially nor tem- porally homogeneous. Rather, the processes in these systems are dynamically (and often systemically) self- or inter-correlated in space and in time. Analyses considering only time or only space may produce misleading results and not reveal the dynamical behavior of the system. In fact, processes such as epidemics, succession, competi- tion, evolution, interactions, and population dynamics, assume that the elements of an ecosystem close to one another in space or in time are more likely to be affected by the same generating process. The lack of spatio-temporally explicit. analytical framework is considered to be a major obstacle to understanding the fundamental mechanisms of such processes. 1.2 Statistical Modeling and Various Space-Time Models In biological and ecological sciences, we usually have incomplete knowledge of the physical mechanisms responsible for the. dynamics, and thus must resort to statistical estimation to represent the dynamic behavior of a system. Statistical modeling is an important procedure for analyzing stochastic systems because it allows us to connect the data to theoretical processes. Furthermore, in principle theoretical models can be used to understand, forecast, or control the behavior of a realistic ecosystem. Fig- ure 1.1 shows how a model is used to represent the data-generating process. The data together with knowledge of the underlying process are used to identify a statistical model. VVlIile modeling, a dataset is used to fit this hypothetical model, following determination of some of the characteristic features (e.g. model types and orders) of the process. The parameters of the model are estimated during the fitting process. Next, adequacy of this statistical model can be evaluated. Finally, the statistical model can be used for forecasting, analysis, control, or other purposes. The objective. is to find a ”parsimonious” model with the smallest number of parameters needed to adequately fit the patterns in the data. Many space-time models have been developed since the initial development of basic theory of space—time models in the 1970s. The following gives an intrtxlnction to some popular space-time models. Kalman filtering has been applied in areas as diverse as aerospace, marine navi- gation, signal processing, demographic modeling, etc. Kalman filtering estimates the state of a system from observations containing random errors. Due to high success, the Kalman filter quickly became an essential role in control theory. In the mid 19908, Cressie suggested that a Kalman filter incorporating space and time can be a powerful tool for modeling space-time processes [17]. Two years later, Huang and Cressie showed the details of the space-time Kalman filter (STKF) [31]. The (.lynamic space—time structure is separable in space and time. In the late 19908, Mardia and Goodall et al. suggested an approach cmnbining Krigging and the Kalman filter, called Krigged Kalman filter (KKF), to model space- time processes [38]. This model uses a time-varying linear combination of spatial fields to represent variations in space and time. The most important feature of this approach is that it reduces the number of parameters to a small number of key parameters. This can be a good method for forecasting; however, these parameters do not provide much information about the underlying mechanisms. Vector Autoregressive Moving Average (VARMA) models (Lutkerpohl [36, 37]) are used to determine the dynamic interactions among multiple time series. VARMA models can be viewed as a subclass of the state-space models that are widely used in physics and engineering. VARMA processes are featured by their autoregressive. and moving average orders. A p—th autoregressive order and q-th moving average order VARMA process can be expressed as P q 2,, :2 Z¢*’Z“k — 20,451-), + e, (1.1) k=1 k=1 where Z(t) = [Zl(t), Z2(t), - ' - , ZN(t)]’ is a kf-dIIIICIlSiOIlal time series vector, EU) is a random noise vector where €(t) = [51(t), 52(t), - - - ,€N(t)]’, (pk and 0,, are unknown coefficient matrices to be estimated. As shown in Equation 1.1, no structured spatial dependence is defined. Instead, the spatial dependence structure is distrilmted among the two sets of parameters, (15,, and 0k. For most ecological and biological space—time variables, spatial correlation (dependence) structures usually can be regularly defined. Thus, the disadvantage of VARMA is that it does not. provide a systematic structure of spatial dependence. If processes or data reveal systematic dependencies between observations in neigh- borhoods, the multi-variate time series models can be refined [27]. Weight matrices can be used to reflect spatial correlation (or spatial dependence) for a given spatial configuration and can be incorporated into a VARMA model to result in Space-Time Autoregressive Moving Average (STARMA) models. In fact, STARMA models can be viewed as special cases of the VARMA models. The extension of univariate ARMA models into the spatial—temporal domain re- sults in a general class of models known as Space-Time AutoRegressive Moving Av- erage (STARMA) models [12, 40]. STARMA models can be used to represent a wide range of stochastic processes that are space—time correlated. In 1980, Pfeifer and Deutsch culminated the collective efforts to extend the Box-Jenkins approach [9] for time series modeling to STARMA modeling [44, 48]. These studies also provided a computational basis for STARMA modeling and analyses. CI! 1.3 Scope of the Work Space-time ARMA modeling potentially has very wide applications. Its primary theory and basic modeling framework have been established since the mid 1970s. However, because of its computational complexity, paucity of appropriate algoritl’nns and appropriate modeling framework, there have been few applications and S[)a(f(-‘,-l.llll(;‘. analyses in the literature. Thanks to modern computing power and the sophistication of statistical algorithms, it is now possible to represent explicitly spatitytemporal dynamics of natural processes in mathematical models and to realize the space—time model in computation. In this thesis, a computationally efficient and widely applicable method is intro— duced. The difficulties encountered while applying the traditional modeling frame- work are reduced by using the proposed algoritlnns and method with the space-time. extensions of the statistical criteria. The thesis will introduce basic STARMA theory, describe how the improvements can be achieved, and give examples of the statistical inferences that can be made using the proposed methods and STARMA modeling theory. The main aim of this thesis is to introduce the improved methods, to establish general space-time ARMA models for empirical data and to show how to establish stochastic space-time models for biological or ecological processes. Because of the presence of noises and errors in natural processes and datasets, it is often not easy to identify and estimate a candidate model from a real dataset using existing methods. To solve these problems, several statistical tools and algorithms were developed in this research for increasing model accuracy and improving computing performance. The secondary purpose of this thesis is to make the stochastic space—time mod— eling theory more accessible to scientists. Since space-time modeling is a relatively unfamiliar field for most scientists, a refined STARMA modeling frammvork is pro- vided for general purpose modeling. Furthermore, to support 011' analysis and to increase the accessibility to space-time modeling for scientists, we implemented a new method and general-purpose software called Integrated Environment for Analyzing STARMA models (or IEAST, please refer to http://fried.for.msu.edu/"ieast), developed explicitly for space-time analysis of processes with two-dimensional lattice data. This is the first general—purpose space-time analysis and STARMA modeling software. Since this work is one of the essential parts of the study, this software will be described and used for data analysis and modeling in this thesis. Finally, a thorough application of STARMA modeling, including statistical space— time modeling and cross analysis of real data for the infectious disease, West Nile virus, is conducted. This application illustrates how the proposed 111ethodology can be tailored to model the spread of infectious diseases. In addition to its potential use for short—term forecasting, this model class contributes to the understanding of the spatio—temporal evolution of disease spreading process, since it can be used to estimate how changes in spreading patterns in some specific locations are propagated to the remaining of the spatial locations. 1.4 Organization of the Thesis This thesis consists of theoretical and application components, according to the fol- lowing chapters: . Chapter 1 ( this chapter) Introduction: An overview of space-time processes and space-time dependence in biological and ecological sciences, space-time statis— tical modeling and various space—time models, and the coverage of the thesis is provided. 0 Chapter 2 Space- Time ARMA Models: The space-time modeling in this re- search is based upon a. well—defined model. Space—Time ARMA (STARMA). This chapter gives a historical background of the development of STARlVIA modeling, the general definition of STARMA, spatial correlation structures, limitations, assumptions, and conventional modeling method/ statistical tools. 0 Chapter 5’ Development of the Proposed Method for S TARMA Modeling, and IEAS T: This chapter explores the details of both the Box-Jenkins modeling approach and the proposed modeling method. Statistical algorithms and tools for improving modeling performance are developed and introduced in this chap- ter. In addition, a. software enviromnent based on the developed Inethod for modeling STARMA models, called IEAST, is implemented and introduced. 0 Chapter 4 Application to the Spreading of West Nile Virus in Detroit Metro Area: This chapter is another important. part of the thesis. Based on the two real datasets, human cases and dead crow data for WNV in Detroit, Michigan in 2002, the dynamic behaviors of both datasets are analyzed and their space- time models are established, using the method and tools discussed in Chapter 3. Statistical inferences are given based on the analysis results. lVIoreover, the cross analysis between these two datasets were examined for statistical inferences on the space-time relationship between human and dead crow variables. Conclusion: Future research and limitations for this method are discussed in this chapter. Stochastic Process Real World Statistical Modeling Statistical Generating Model Figure 1.1. Relationship between data-generating and statistical modeling procedures Modeling oooooooooooooooooooooooooooooooooooooooooooooooo Procedures 10 > Forecasting/ Analyzing/ Control CHAPTER 2 Space-Time ARMA Models 2.1 Historical Background Research in statistical / ecological models that describe the spatio—temporal evolution of a single variable or multi—variable relationships in space and time started in the mid 19708. It has increased significantly during the last twenty years, because it is closely related to progress in computer technology and the prevalence of remote sensing for collecting large space-time datasets. Cliff and 0rd were among the first to establish a model for interactions in space and time [14, 15], and since then several teclmiques have been developed corresponding to different inferential needs and data types. The STARMA model class was first presented in the literature in the late 1970‘s by Bennett [6]. Since then, it has been applied to spatial time series data from a wide variety of disciplines such as river flow, traffic control, and spatial ecemnnetrics. A series of studies on the basic space-time modeling theories was completed during the beginning of 19808. Areian defined tit-dimensional time series, derived properties, and 11 applied to the examples of space-time series [4]. In following papers, univariate time series autoregressive models [3], moving average models [55], and mixed models [42], i.e. autoregressive moving average models, were extended to viz-dimensional cases (STARMA is a special case of m-dimensional time series). Furthermore, some. of the statistical properties of these models were investigated by Pfefier and Deutsch [46, 47]. Recently. Giacomini and Granger pointed out that the STARMA class of models can be derived through a transformation of the Vector Autoregressive Moving Average (VARMA) models [27]. In fact, the transformation is a restriction related to the systematic spatial correlation structure as revealed by a set of weight matrices [27]. STARMA models can represent well most applications at large spatial scales, but for the applications involving only a few observations across space it may be too parsimonious. For the related univariate ARMA modeling methodology, Box and Jenkins in 1970 approached the identification problem for time series by using temporal autocorre- lation and partial autocorrelation analysis [9]. This method consists of identifying the characteristic behaviors of the correlograms for temporal autoregressive, moving average, and mixed processes. Martin and Oeppen then extended the identification procedure in the Box-Jenkins approach to the space-time case, also using space-time correlation functions [40]. At the same time, Bennett proposed an alternative proce- dure for model identification for the spatial time series [5]. Various aspects of identi- fication have been addressed by Bennett [5] and l\-*Iartin and ()eppen [40]. Pfeifer and Deutsch in 1980 provided a comparison of accuracy and precision of various estima- tion methods for the parameters in space-time autoregressive models, using computer 12 simulations [45]. The three-stage iterative model-building philosophy conunonly referred to as the Box-Jenkins approach [9] for building univariate time series models has been adapted for use with STARMA models. Pfeifer and Deutsch (1980) [44, 45, 48] were among the first to develop space-time modeling techniques for lattice spaces in the context of STARMA models, and they illustrated the model-building details for the identifica- tion, estimation, and diagnostic checking of the STARMA model, using an iterative three—stage procedure. The extension of univariate ARMA time series models into the space-tin'ie domain results in the general model class of STARMA models [12, 14]. These models apply to a single random variable. observable at N fixed sites or locations in space at discrete points or periods of time, t = 1, 2. T. They are of value for descriptive and fore- casting purposes when the observed system exhibits spatial autocorrelation defined by Cliff and 0rd [13]: ’If the presence of some quality in a county of a country makes it presence in neighboring counties more or less likely, we say that the phenomenon exhibits spatial autocorrelation.” 2.2 Model Definitions STARMA is a space—time extension of the ARMA model [9]. The STARMA model class is characterized by linear dependence lagged in space as well as time. We summarize some of the results of Pfeifer [44], using essentially the same notation. STARMA expresses Z”, the observation of the space-time random variable at site 13 2' (i=1,2,...,N} and time t, as a weighted linear combination of past observations and random noise inputs, which may be lagged both in space and time. Let Z, = [Z1,_,, Zg‘t, - - - , Z N,,]' be the N x 1 vector of observations at time t, where N is the total number of sites in space. The concept of a spatially lagged variable requires an ordering of the neighbors of each site. A definition of spatial order appropriate for regularly spaced systems is found in [7]. To formulate the spatial relations among sites, spatial weight matrices W“) are used. Assume weight matrix W“) has elements mg) that are the wr-iightiug contributions of site j to site i, and which are nonzero if and only if site i and j are l-th order neighbors in space. Then, the general STARMA model can be expressed in the following form: 1' q 3 WWW.-. -— Z Z auw<”e,_k + e. (2.1) =0 k=1 (:0 p z,=Z k=ll where p and r are respectively the maximum autoregressive temporal and spatial orders, q and s are respectively the maximum moving average temporal and spatial orders, 96“ and 0“ are respectively the autoregressive and moving average. paranwters at temporal lag k and spatial lag l, Wm is the N x N weight matrix for spatial order I, and at = [51,t352,ta - - - ,8N‘1]/ is the random noise vector at time t. The weights w]? should reflect an ordering or proximity of spatial neighbors. Figure 2.1 illustrates an example of a spatial order definition. The first order neighbors (corresponding to w(1)) are those that are closest. to a. given site. (the central gray circles). The 2nd order neighbors are farther away from than the 1”" order neighlmrs, 14 but closer than the 3"! order neighbors. There are three major mode] types (STAR, STMA, and mixed models) defined for general STARMA models. A process is said to be a Space-Time AutoRegressive process of temporal order p and spatial order 7‘ if q = 0 (named as STAR(p. 1)). A STAR process can be expressed as P M: Zr = Z Z ¢sz(l)Zt-k + 51- (22) [:0 k2] A Space-Time Moving Average process is of temporal order q and spatial order s if p = 0 (named as STMA(q, 5-)). A STMA process can be expressed as q "7k z, = e, — Z Z 9,,,w<’)e._,,. (2.3) k=1 [=0 The mixed model combines both autoregressive and moving average effects (if p > 0 and q > 0), and is named as l\t"lixed(p,q,'r,s) (its mathematical form is given by Equation 2.1). STAR and STMA model are popularly used in practice. Not only can many practical stochastic processes be simply attributed to either STAR or STMA. but also there exist primary statistical differences between them. 15 2.3 Spatial Correlation Structure and Weight Ma- trices As shown in Equations 2.1-2.3, spatial correlation structure is specified by a. set of weight matrices (Wm). There are ma.r(r, s) + 1 weight matrices, i.e. W(0),W(1), - - - ,W("’"""(""“)), for correlation structures of each spatial lag for a Mixed(p,q,r,s) model; r + 1 weight matrices for STAR(p,r); s + 1 weight matri- ces for STMA(q, s). Specifically, Wm) is a unit matrix IN. If the total number of sites in space is N (i.e. the length of the observation vector Z), the size of all weight matrices for the configuration is N x N. The nv-th row of the weight matrix for spatial lag l is the distribution of the relative contributions of the sites in the neighborhood of spatial lag l to the site 77. at some temporal lags. In (1) other words, if the weights ll’ij are the elements in W”), it is nonzero only if sites 2' and site j are neighbors for spatial lag l with i NW—i . d +‘”>0f '111' =1 21111,]. _ an. a}. _ or a l,j,.. J: Weights mg.) for each spatial lag 1 can be selected to reflect physical connections and relationships among sites (or locations) of the system under study. The weight matrices can be specified by model builders based on prior knowledge of the data observed, or can be assigned to be isotropic and generated systemati- cally. Because weight matrices define the relationship between each site in a general sense, STARMA model can be easily adapted to various spatial dimensions (eg. 16 one-dimensional application, as in modeling the dispersal of species along a river, or three-dimensional application, like modeling air conduction and circulation) or special spatial arrangements (e.g. anisotropic weighting for directional spatial corre- lation, like modeling spatial tendency of migration of a. species). 2.4 Box-Jenkins Modeling Approach The classic (1970s) time series analysis uses a Box-Jenkins approach which is a general procedure for modeling and forecasting stationary autoregressive and moving average processes. The main output from such an approach is a regression model explaining current values of the series in terms of past values. The coefficients in the. model can then be used to forecast the series into the future. This approach involves identifying an appropriate ARMA model, fitting it to the data, and using the fitted. model for forecasting. One of the attractive features of Box-Jenkins approach for forecasting is that ARMA processes are a very rich class of possible models. It is usually possible to find a model which provides an adequate description for the data. The Box-Jenkins approach consists of iterative steps of model identification, parameter estimation, diagnostic checking and forecasting (as shown in Figure 2.2). Iterations of these steps are then used to find increasingly better solutions. The method for STARMA modeling in our research is based on a space-time extension of the Box-Jenkins approach in [48]. In chapter 3, improvements and details of STARMA modeling methodology are described. 17 2.5 Space-time Autocorrelation and Partial Auto- correlation Functions There are two questions that arise in modeling procedures. First, what type (STAR. STMA, or mixed) of the space-time model should be used? Secondly, what are the spatial and temporal orders of the model? Because there are so many potential candi- dates, it is advantageous to have a computationally efficient method for identifying the type and orders. As shown in Box and Jenkins [9], the popular identification method for time series is to use the estimated autocorrelation and partial autocorrelation functions. In an analogous way, the following statistical tools, space-time covariance. autocorrelation, and partial autocorrelation can provide a systematic method for re- vealing the characteristics (type and orders) of the underlying process for a given dataset. The space-time covariance between (”I and A?“ spatial order neighbors at time lag s can be defined as in [44] [WmZtl’IW‘k’Zm] WINS) = Er. N (2.4) Based on the definition of space-time covariance in Equation 2.4, the Space-Time Autocorrelation Function (STACF) between 1’” and Ir” spatial order neighbors at time lag s can be defined as 71143) . t. : . 2.5 [)M (9) 7H(O)’7kk(0) ( ) 18 Having the space—time autocorrelation function from Equation 2.5, to fully imple— ment the identification stage. in modeling procedure the corresponding partial auto- correlations have to be estimated. The Space-Time Partial Autocorrelation Function (STPACF) can be found in the following way. Given temporal lag m and spatial lag n, the best linear predictor of Z, can be expressed as 171 2(t)lm.n = Z Z ¢le(I)Zt—k- (2.6) k=l 1:0 For every m (=1...maximum temporal lag) and n (=0...maximum spatial lag), a set of coefficients {mu 2 {am/c = 1...?7‘z.,l = 0...n} can be found. For given temporal lag m and spatial lag n, the STPACF is the element an", in the set Em". En,” that. is the solution of the space~time Yule—\Valker equations for maximum temporal lag m and spatial lag n [40, 48]. In other words, the space-time partial correlation function can be found by successively fitting STAR(MaxT=m,MaxS=n) models for m =1, 2. 3, and for n =0, l, 2, and picking out the estimates of the last coefficient from this sequence of models. Computationally, the partial autocorrelations can be obtained by solving the space-time Yule-Walker equation for each spatial order 0, 1, 2, and for each temporal order 1, 2, using the recursive method shown by Durbin [19]. 2.6 Assumptions For space-time modeling, there are some assumptions made, including spatial reg- ularity [25], stationarity [47], invertibility [47], constant correlation structure, and 19 normality. A spatial system has spatial regularity if the spatial structure is regularly defined over entire space, e.g. as in a lattice. Without spatial regularity, it is not meaningful to define spatial lag structures [25]. Thus, we generally will not obtain spatial station— arity without spatial regularity. In practice, it is extremely difficult to define a. spatial lag structure so that STARMA models are stationary and correlation functions well defined without the assumption of regular distribution of locations. The most important assumption in STARMA modeling is the validity of stationar- ity. A space-time process is called stationary if its statistical properties do not change with time or locations. However, it is often impractical to assume that the observa- tion data is generated by a stationary process, especially in space. In the natural world, seasonality, trends, and clusters can cause non-stationarity. Generally, some. preprocessing or transformations such as differencing, de-seasonalizing, de—trending, or logarithmic transformation must be performed in advance to obtain a stationary series before a STARMA model can be fitted. It makes no sense to try to fit a data generating from strongly non-stationary process to a stationary model like STARMA. Without stationarity or approximation of stationarity we cannot make detailed sta- tistical inferences about the properties of the underlying processes. Furthermore, by ensuring the property of spatial rerflication, spatial stationarity allows us to gain more statistical power [25]. Using a stationary model to fit a process without stationarity is just like trying to characterize a system in which the characteristics of the system are constantly changing. There is no way to look inside or to characterize such a system from its inputs and outputs. 20 With the assumption of normality, only weak stationarity (second-order) is needed for STARMA modeling to achieve strict stationarity. Weak stationarity is defined as follows: Given a space-time process {Zml t E Z, .7: E R"}, it is weakly stationary if: 1) second moment of Z” is finite for all t and :17: 2) first moment of ZN is indepen- dent of t and at: 3) cross moment of ZN and ZIHHh depends only on s and h. Weak stationarity plus an assumption of normality are sufficient to produce strict station— arity. For a strictly stationary process, we have the same probability distribution for all time and locations. If the probability distribution is not independent of time and locations, it is difficult to define statistical properties of the process, e.g. level, variance, covariance, correlation, and even hypothesis test. Among these properties, the critical ones are correlations and hypothesis tests, which are necessary during STARMA modeling and general space-time analysis. To be stationary and having a unique solution, as in the univariate case, the causality and invertibility conditions of the process nmst be satisfied. If p = 0 in Equation 2.1, the model is a pure STMA model (as in Equation 2.3) and Z, is al— ways stationary. When a STARMA model has p autoregressive terms (p > 0), Z) is stationary if and only if all roots of PM (1615 I —' Z: Z ¢k1W(I)IIT—k = 0 k=1 [:0 are inside the unit circle, that is. [r] < 1. If q = 0 in Equation 2.1, the model is a. pure STAR model (as in Equation 2.3) and Z, is always invertible. When a STARMA 21 model has q moving average terms (q > 0), Z, is invertible if and only if all roots of m 1: (1(21‘. I— 6HW(1).‘IT_k =0 (1 k=1 (:0 are inside the unit circle, that is, [:1] < 1 [48]. However, stationarity does not guarantee constant spatial correlation structure. Processes that can be modeled by STARMA are the stochastic processes whose prop- erties, such as correlative structures, do not vary with location and time. Hence, we need to further assume that the spatial correlation structure is constant, not varying with either time or locations. It is assumed that each observation of the space-time series has the same expecta- tion function, standard deviation, and probability distribution function. We furtlwr assumed a noise error component, which is a sequence of uncorrelated random vec- tors with a constant distribution (Gaussian), constant variance, and zero mean. The random noise vector at is normally distributed at time t and satisfies the following characteristics: at N1V(0,0‘2IN) 0213;. s = 0 ’ _ El€(€(+sl _ 0, otherwise This assumption of normality assures that weak stationarity is sufficient for a uniform probability distribution for all times and locations. Although teclmically the normality assumption is required, little to nothing is known about how strict 22 violations may affect inferences. 2.7 Finest Lag Structure and Selection of Spatial Correlation Structure Although spacetime ARMA appears to be a direct extension of time series ARMA, several properties of time series models do not apply to STARMA models. In par- ticular, some important characteristics of STARMA models are only valid for the finest lag structure [30]. A finest spatial lag structure is a special case of corre- lation structures when, for each pair of site i and site j, there exists s such that the neighborhood for spatial lag s of site i contains only site j. Strictly speaking, properties such as stationarity, causality, invertibility, autocorrelation functions, and Yule—\Nalker equations, are only valid for finest lag structure [30]. Spatial correla- tion structure significantly affects the a1V1t0(.'.()rrelation function, which may further influence the determination of model types and orders during model identification. The following points should be kept in mind while defining spatial correlatirm structures. First, choosing a particular lag structure other than finest lag structure limits the generality of the models, so care must be taken in defining spatial correlation structure to ensure that useful models are not eliminated from consideration. Second, spatial correlation structure should be defined in a uniform manner across space to ensure that the STARMA models are stationary and the correlations are well- defined [30]. It may often be desirable to use a coarser lag structure in order to reduce the number of parameters that must be estimated. However, having a coarser lag struc- ture restricts the generality of the class of STARMA models. Therefore, the choice of an appropriate definition of lag structure is an important first step in the. process of selecting an adequate but parsimonious model. Knowledge about the process un- derlying the data is helpful in deciding on an appropriate lag structure or correlation structure. For instance, if the interactions among locations, for the process being studied, depend only on the distance not the direction (i.e. isotropic), then the lag structure like the following can be a reasonable choice. 32123 3210123 32123 If both distance and direction influence the interactions among locations, the following definition (considering the distance and the direction in south, north, west, and east) might be desirable. 24 11 12 12 12 12 12 12 OCOCOOOCOONI Orb-##CJJN OWN CI! CECEGCD C}! 10 10 10 10 10 10 OOOOOOO OOOOOOO OOOOOOO OOOOOOO OOOOOOO OOOOOOO OOOOOOO 00.0000 0000000 0000000 OOOOOOO OOOOOOO OOO®OOO OOO®OOO OOO®OOO OOO®OOO 0000000 0000000 OOOOOOO OOOOOOO OOOOOOO OOOOOOO OOOOOOO OOOOOOO OOOOOOO OOOOOOO OOOOOOO OOOOOOO 1st order 2nd order 3rd order 4th order Figure 2.1. An example of the definition of spatial orders/ lags in systems with two spatial dimensions. 26 <1Dataset > Identification f; l [ MEN f l Diagnostic Checking [Modify Model] No Good? Yes [ Forecasting ] Figure 2.2. Box-Jenkins modeling approach. 27 CHAPTER 3 Development of the Proposed Method for STARMA Modeling, and IEAST In STARMA modeling, there are some problems that have. impeded applications since the theory was developed, including: no robust modeling method; no appropriate statistical tools/ algorithms; and no flexible or general purpose software. This chapter describes the proposed STARMA modeling method with the related extension of statistical tools and algorithms. This chapter will provide an introduction to the conventional method as well as detailed description of the refined method. At the end of this chapter, the general-purpose software IEAST for STARMA modeling will also be introduced. 28 3.1 Specification of Spatial Correlation Structures and Derivation of Weight Matrices The first step of STARMA modeling is the specification of spatial correlation struc- ture (i.e. assignment of spatial weight matrices, or SVVM for short). A systematic. procedure used in this research to specify spatial correlation structure. is introduced in this section. The procedure for specifying spatial weigl'it matrices consists of three steps: (1) specify Spatial Order Definition (SOD) matrix; (2) specify Spatial Relation Ma.- trix (SRM); (3) based on (1) and (2), weight matrices (SW M) can therefore be determined and generated. In summary, the specification follows the sequence of SOD —> SRM —> SVVM. SOD is a matrix defining the spatial order of each ’related’ neighboring cells relative to the central cell (i.e. the cell which is under consideration for spatial correlation with neighboring cells). The matrix size depends upon how large. the correlated spatial neighboring area must be as needed. For example, if we define. the first spatial order cells as the cells which are at a single ”rook‘s move" away; the second spatial order cells as the cells which are at a single ”bishop‘s move” away: and the third spatial order cells as the four cells which are exactly two cells away from the central cell in the direction of east, west, south, and north; and so on, up to fourth order, we can specify the SOD matrix as 29 -_ 4 3 4 —- 4 212 4 SOD: 3 1 0 1 3 . (3.1) 4 212 4 -_ 4 3 4 ‘, The number of cells that can be defined for a specific spatial order is arbitrary. The central element of SOD represents spatial order zero, that is, it is the cell (under consideration) itself. In this example, four neighboring cells are. defined in each of the first three orders of the central cell, but. eight neighboring cells are defined in the fourth order. The definition of SOD needs not be symmetric. and its arrangement could depend on the spatial correlation structures of the applications. Once the SOD is specified, the weighting distribution in a given spatial order can be defined. The weighting distributions of each spatial order are. defined by a. matrix called Spatial Relation Matrix (SRM). The sum of the weight distribution of every specific spatial order should be unity. The coefficients in SRM represent the importance (in the sense of spatial correlation to the central cell) of every neighboring cell in each specific spatial order. For uniform weights (isotropic), a SRM for the SOD defined above can be defined as — — 125 .25 125 — l 125 25 .25 25 125 81m: .25 .25 0 .25 .25 . (3-2) 125 25 .25 25 125 _ — 125 .25 125 — ] 30 '3 I]; By comparing SOD and SRM, we. can definitely associate the spatial order of any given cell with its corresponding weighting (in that spatial order) in the. correlated neighborhood. The above is an isotropic (directionless) example. We can also use a directional spatial correlation structure, i.e. anisotropic correlation structure. The following is an anisotropic example (a biased distribution towards the north in the first and second spatial orders) of SRM for the SOD given in Equation 3.1. — .125 .25 .125 — .125 .35 .40 .35 .125 SEN: .25 .25 0 .25 .25 (3.3) .125 .15 .10 .15 .125 — .125 .25 .125 — Every two—dimensional spatial matrix can be rearranged into a linear vector (as the. column vector Z, in Equation 2.1) for mathematical manipulations and computations. For example, an observation in 4 X 4 spatial matrix can be rearranged into a linear vector Z as following. a ' T a b c a' b c f h c ” -—» = z (3.4) I J A? I i - m n o 1) J I- p d 31 Weight matrix W“) is actually a. linear transformation matrix for spatial lag l (or spatial order I). Consider WmZ: the observation vector Z is spatially-weighted and transformed into vector Z”) with the weighting distribution at spatial order I by multiplying Z with the weight matrix of spatial order I, W"). The sum of the transformed Z for all spatial orders can be written as Z Z“). t [- T " ‘ '7 ' " (l (l b b" 125 35 .40 125 (I) C Cf (I) w Z _ = = Z (3 5) .- _[ h p d - pi! — In the above example, the n-th row of the weight matrix for spatial order I (i.e.. w(l)) is a rearranged linear row vector of the relative contributions of the cells in the neighborhood of spatial order I to the cell it. Hence, the weight matrices can be generated systematically in the following way. For the given SOD (Equation 3.1), SRM (Equation 3.3), and space dimension of 5 x 5, we can now construct the weight matrix row by row for the first spatial order weight matrix W“). Let us look at a simple example. After removing the coefficients for the spatial orders other than the first order, the matrix SRM can be rewritten as 32 , _ 0 .40 0 SRM = .25 0 .25 . (3-6) L 0 .10 0 . In order to use a weight matrix based on the above first order SRM to transform an observation (say vector Z) for each cell (from the 1st to the 25th element) in a space with dimension of 5 x 5, we need to have a weight matrix W“) with 25 rows, and it can be obtained as follows. To obtain the first row of W“), put. the SRM into an empty (filled with zeros) 5 x 5 matrix by aligning the central cell (bold Us in the following equations) of SRM to the first cell of the 5 X 5 matrix. Everything outside the 5 x 5 matrix is discarded. Thus, we have. 0.25000 .10 0 000 Row1-> 0 0000 —*[0.25000.1000---0] 0 0 000 L0 0 000- For the second row of W“), the SRM is shifted for a cell to the right of the. first cell so that the central cell of SRM aligns with the second cell of the 5 x 5 matrix. Keep doing this until the 25th row and have 33 p.25 0 .2500- 0 .100 00 Row2 —+ 0 0 0 0 0 —+ [.25 0 .25 0 0 0 .10 0 0] 0 0 0 00 _0 0 0 00. ’0 25 0 25 l 0 0 .100 Row3 —» 0 0 0 0 0 —> [0.250.25000.10-..0] 0 0 0 0 0 _0 0 0 0 0_ [000 0 01 000 0 0 Row25_+ 000 0 0 —> [00-~0.40000.250]. 000 0 .40 _000 .25 0‘ Integrating all these 25 rearranged row vectors above, the weight matrix W“) is thus obtained. The same procedure can be also applied to spatial order two, three, and so forth, to generate Wm, Wm, etc. Finally, the whole. set of weight matrices W”) (l = 1,2, ~ - - , max spatial order) for the given SOD and SRM can be generated. This procedure has be automated in IEAST. 34 3.2 Identification Stage The first stage of the modeling procedure is the identification stage, in which the initial most probable underlying model type (STAR, STMA or Mixed) and its maxi- mum spatial/ temporal orders are determined. There are diagnostic differences in the STACF and STPACF among pure STAR, pure STMA, and mixed models. The model type is determined based on whether the curves of STACF or STPACF drop abruptly, or if both of them decay in an exponential manner along the time (or spatial) axis (e.g. Pfeifer and Deutsch [48]). In the following sections, we call the abrupt dropping behavior "cut-off" and the exponential decaying behavior "tail-off". Technically. we can distinguish these two behaviors by their first order derivatives of the autocor- relation functions. The derivatives of tailing-off functions generally are exponential shape, and the derivatives of cutting-off functions have abrupt change around the cut-off point. Conventionally, the method of identifying the process type was simply a qualitative comparison of observations to theoretical behaviors. In the. following sections, I have developed some quantitative metrics and methods to set up a sys- tematic identification process that significantly improves computational performance. Furthermore, the spatial and temporal order can be determined by investigating the lags (in space or in time) at which the STACF or STPACF cuts-off. In a manner analogous to that of the univariate ARMA, the subclasses of STARMA family are characterized by their distinct STACF and STPACF behav- iors. A pure STAR model exhibits a STACF that tails-off in both space and time, and STPACF that cuts-off after certain lags in space and time. In contrast, a pure 35 STMA model is characterized by a STACF that cuts-off after some temporal and spatial lags and a STPACF that tails-off spatially and temporally. A mixed model exhibits the spatial and temporal tailing-off of both the STACF and STPACF. Two examples shown in Figure 3.1 and 3.2 are to illustrate the cut-off/tail-off behaviors of STACF/STPACF of two typical simulated data. These data were. generated by picking the last 100 generations from a total of 5000 simulated generations. The two-dimensional space was assumed to be 8 x 8, i.e. 64 cells in total. Uniform weight. matrices (or isotropic spatial correlatirm structure) were used for all spatial orders (from first to fourth). The random noise vector was 5, ~ N (0, o2 = 0.00125). In summary, for STAR(2.1), the STPACF cuts-off at temporal lag two and at. spatial lag one; STPACF tails-off. On the other hand, for STIV’IA(1.1), the STACF cuts-off at temporal lag one and at spatial lag one and the envelope for STPACF tails-off. These two examples illustrated how the. model type is selected and how the temporal / spatial orders are determined by observing the behavior of STACF and STPACF. 3.3 Estimation Stage The second stage is to estimate the parameters in the model identified in terms of the model type and orders. In general. parameter estimation is to minimize the following sum of squared error function (or maximize the likelihood function) (Equation 3.7) to find a set of maximum likelihood estimates. The maximum likelihood function to 36 be optimized is 2 T p A1,. q 171),. 5(3) = Z Zt " Z Z $k1WU)Zt—k + 6kIW(l)Et—k (3.7) t=l k=II=0 k=ll=0 where T is the number of observations in time, Z, is the observation vector at time t, 6t is the random error vector at time t, and fl = [$10~$11~"°a$pxaé10~é1u”Ham],- This is a quadratic nonlinear optimization problem. For linear models (i.e. STAR for STARMA case), these maximum likelihood es- timates can be found by applying linear estimators. For example, in a model of STAR(maxT=2,maxS=1), all observations can be substituted into the model equa- tion and arranged as in Equation 3.8. . - F . - - 2h () O () O f q 61 1 $10 22 Z1 W( )Z1 0 0 45 52 23 = z2 w<1>z2 z1 w<1>z1 ” + 53 (3.8) $20 $21 27‘ ZT—‘l W(I)ZT—1 ZT—2 WmZ-T—2 ‘ ‘ ET The equation above can be simply expressed in linear model form as Z = XQ +E, and the best linear unbiased estimator is Q = (X’ X)‘1X’ Z, where the observation vector Z = [Z1.Z2, - - - ,ZT]’, model parameter vector Q = [tbm.(,5n,0§20.962.]’. error vector E = [51.52, - - ' .51], and X is the large matrix in Equation 3.8 [48]. Because of the nonlinear nature of STMA and mixed models, the linear estimator above is not appropriate. In this research, the parameters for nonlinear processes (i.e. 37 STMA and mixed) are estimated by maximizing Equation 3.7 and searching through a quadratic surface. A quadratic optimization algorithm, l\larquardt’s algorithm [39]. is used to maximize the likelihood function (or minimizing sum of least squared errors) and to find the parameter estimates. lVIarquardt’s algorithm is a Gauss-Newton based algorithm for least squares optimization. However, as one may commonly encounter in such optimization problems, it is important to locate an appropriate starting point for the optimization process. This is especially critical for multi-variate nonlinear optimization problems. During this research, it was found that the optimization pro- cess for parameter estimates, in most cases, either converges to local optima or does not converge. Thus. to avoid reaching a local maxima and to reduce the number of iterations needed during the optimization process of the maximum likelihood func— tion, a preliminary stage, pre-estz'matz'on. is implemented to calculate an appropriate starting point for Marquardt‘s algorithm. In this study, I have extended and imple- mented the Hannan-Rissanen algorithm from univariate case [29] to space-time case for calculating the initial points for STMA and mixed processes as shown below. The space-time extension of Hannan-Rissanen algorithm has three iteratim steps. 1. A high order STAR model is fitted to the data using the space-time Yule-Walker equation. Then, we have the following approximate model Z; = i : 7),,1W(”Zf_k + E; k=1 [=0 where {fikllk = 1 - - - u,l = 0 - - - v} are the Yule-Walker estimates. 38 2. The estimated random noise vectors can be computed as U 1’ e __ t t 2 :2 :_~ I t . 1:21 1:0 3. Once the estimated random noise vectors é,.t = m + 1,---,T have been found from Equation 3.9. pre—estimates of the model parameters, 61 = A A I I O I [cp',1p]’ are deternnned by least squares linear regress1on of Z, onto the space {Zt—l. Zt_2. - ' ’ , Z¢_,,. 5-1. €143. ' - - ,é¢_q}.t = m + 1, . - - , T. By minimizing the sum of square errors T ,\ 2 P k q mk A 5(5)) = E Z: — E $L~1szt_k+ E E ’l/u-iwméz—k (:171.-+-l k=1 [:0 k=1 (:0 with respect to d, we can obtain the space-time extension of the Harman- Rissanen estimator a = (X’xrl x’e where 39 A 'A A A A! 7! 7v I a = [3010,99113' ' 399k191$10~ (+11. " Nil-kl] s Q : [Zm+1~Zm+2s"'~.ZT]I~ l' Zm Zm-l-l From simulation Using the pro—estimate d as the. initial point for the optimization, the final least A ’ A I I U Q 0 t O O I 0 [Q , G ]’ to uunnnlze Equation 3.7. The maxnnum likellhood estlmate of the residual variance (32 is thus W(1)Zm W(1)Zm+l (3,2 W(r) z"! Zm—l w(r)zm+l Zm improves the robustness of the modeling method. TN 40 5(, e) ooooooooooooooooooooooooooooooooooooooooooooooooooooooooo results, the space-time extension of Hamian—Rissanen algorithm is very efficient and accurate. The deviations of the parameter pre-estimates d from the correct values are generally lower than 10%. This algorithm greatly reduces the possibility of (h-‘larquardt’s algoritln‘n) converging to a wrong local optima. and hence squares estimates of the model, Q and 9, are obtained by finding the best ,8 = where T is the number of points in time of the observations, N is the number of sites in each observation, S (Q, 9) is the residuals sum of square errors, Q is a vector containing the estimates for parameters in autoregressive terms, and ('9 is a vector containing the estimates for parameters in moving average terms. 3.4 Diagnostic Checking Stage Once the model type, orders. and parameter estimates are determined, the acquired model can be used to forecast or to analyze the future behavior of the system. How— ever, to prevent 1.1nacceptable forecasting errors, further diagnostic checks must be made to ensure the selected model is appropriate. The methods for checking the adequacy of the models include: residuals" autocorrelation analyses; residuals‘ ran- domness testing, residuals” variance checking [48, 49], and statistical significance test- ing of parameter estimates [48]. In this study, residuals" space-time autocorrelation analyses, significance testing of parameters, measures for the model’s adequacy, and residuals’ variance checking are used for diagnostic checking. Regarding the residuals of the STACF and STPACF, inadequacy of a candidate model can be revealed by the presence of some significant correlation values. If a can— didate model adequately represents the data, residuals should be white(uncorrelated in space and time) noise. In other words, all autocorrelations (in space and in time) at nonzero lags are equal to zero (or approxin‘iately equal to zero, in practice). A candidate model is also checked for any unduly complexity, i.e. having sta- tistically insignificant parameters. Thus, statistical significance of model parameters 41 is tested for this purpose using the hypotheses that a parameter of interest 13,. is zero while keeping other parameters unrestricted. That is, let B“ = [81,132,- ".131; = O,Bk+1, - - -]’. The test for this hypothesis is based on the statistic [TN — we) — Ire-imam — 8(3)] N F1 mam-(é) (3.“) 5(fi) where T is the number of points in time in the dataset, N is the total number of sites in the dataset, S ( B") is the sum of square residuals for parameter B with the restriction Bk 2 0, S (B) is the sum of square residuals for parameter B without restriction, Q is a vector containing the estimated parameters for autoregressive terms in column vector form, Q is a vector containing the estimated parameters for moving average terms in column vector form, K(Q) is the number of elements in Q, K(C:')) is the number of elements in ('9. Any estimated parameter shown to be statistically insignificant. Should be removed from the model, and the number of spatial or temporal orders may be reduced, giving a simpler model that is then considered as the next candidate model in an iterative modeling procedure. The results of autocorrelation analysis may be obscure for model identification. and high order models are not necessarily desirable although the residuals' variances are small. As shown in Box and Jenkins [9], an autoregressive model with invertibility can expressed as an infinite sum of moving average terms (or with very high maximum orders); and a moving average model with causality can be expressed as an infinite sum of autoregressive terms (or with very high maximum orders). In practice. we need to consider another model type (e.g. STMA) if the selection of current model 42 type (e.g. STAR) causes very high orders in space or in time. To sun’imarize the model fitness and the model parsimony, we need to introduce some control criteria. to provide measures for evaluating candidate models, as well as to penalize the excessive use of extra parameters in the fitting of models. I extended the space-time extension of bias-corrected Akaike’s information crite- rion [1], referred to as the AICC, to be. the criterion in determining the model type. orders and checking adequacy of parameters. It can be derived directly from Hurvich and Tsai [32] by using the likelihood function for STARMA [48] S(5) 202 L(B, 02|Z) = (27r02)_%~:e.rp(— ) (3.12) and obtained as Q) + TNIK<) + K(é)] (., 1,) 5:2 TN — K() — me) — 1' ' ' " AICCSTARMA = TN10£‘(27”32) + In the definition of AICC, the last term is the penalty for high order models, which makes the value of AICC higher than that of lower order model. Thus. while de- termining an adequate model. we are searching for a model (given type. orders. and parameters) with a minimum value of AICC. However, because AICC has a tendency of overfitting autoregressions [33, 51], a second criterion, Bayesian Information Criterion (BIC), is suggested to correct this 43 tendency [2. 50]. The BIC is also extended to the space-time case as S . - BICSTARAIA = TNIOgmfléz) + (B) + (I\'(Q) + K(@)) - Iog(TN). (3.14) 0-2 We illustrate the perforn'iance of the model type selection with AICC and BIC, by using 150 Monte Carlo simulated datasets based on STAR, STMA and Mixed processes (50 datasets for each subclass). In this illustration, we used each of the four criteria: AICC only; BIC only; -AICC*BIC; and variance of residuals. The parameters of each simulated model were randomly chosen, with higher weights on the parameters having lower spatial or temporal orders. As expected, using AICC or variance as the only criterion has a bias toward indicating (overparameterized) mixed model (Table 3.1). The AICC is useful when the candidates are only STAR and STMA (excluding mixed model). If we would like to select one out of these three subclasses, BIC per- forms much better than AICC. As is shown in Table 3.1, all of the 50 simulated STA R datasets were correctly classified as STAR, i.e. 100% accuracy. If the given simulated datasets are STMA and STARMA, the accuracies are 86% and 66%, respectively. Furthermore, we found that using the product of AICC and BIC as the criterion can provide even better performance (Table 3.1). The accuracy of model type determina— tion for all subclasses can be as high as 78%. In summary, these methods provide a reliable quai'ititative measure for (‘letermining model choice, whereas previous research relied on subjective qualitative behavior of the STACF and STPACF. With regard to variance checking, it has been shown by Pfeifer [44] that if a 44 Table 3.1. Accuracies of model type selection using variances, AICC, BIC. and — AICC*BIC based on 150 Monte Carlo simulated datasets. Using v(u‘i(mcc.s of residuals Sinmlated datasets Datasets identified as based 011/ STAR STMA Mixed STAR 4% 0% 96% STMA 4% 6% 90% Mixed 8% 2% 90% Using AICC Simulated datasets Datasets identified as based on / STAR STMA Mixed STAR 16% 0% 84% STMA 4% 6% 90% Mixed 8% 2% 90% Using BIC Simulated datasets Datasets identified as based on / STAR STMA Mixed STAR 100% 0% 0% STMA 4% 86% 10% Mixed 18% 16% 66% Using -AICC*BIC Simulated datasets Datasets identified as based on / STAR STMA Mixed STAR 100% 0% 0% STMA 4% 78% 18% Mixed 16% 4% 80% process is pure white noise. then 1 z m ($.10) zza,riancc(/3m(3)I where 610(5)) is the sample STACF as defined in Equation 2.5. That is. if the estimated model appropriately represents the dataset, the variance of the sample STACF of residuals should satisfy Equation 3.15. 3.5 A Refined STARMA Modeling Procedure Combining the criteria and algoritlnns developed. a STARMA modeling procedure was developed and used for analyses in this research, as shown in the. flowchart in Fig- ure 3.3. The intention of the development of this refined modeling procedure has two purposes. The first is to reduce the computing time for optimization in the estimation stage, and to increase reliability during parameter estimation, especially for nonlinear models (STMA and mixed). The second is to overcome the difficulty in model identi- fication resulting from the ambiguous behaviors of some STACF/STPACF. These two advances greatly promote the feasibility and robustness for automating the STARMA modeling procedure computationally. This refined procedure. is similar to the. three-stage iterative procedure in [48]. in which STACF and STPACF are used to identify model types/orders and ordi— nary least square is used for parameter estimation. For our approach. in those cases where the model type cannot be. qualitatively identified by STACF/STPACF. the 46 criteria AICC and BIC are used to iteratively determine the best candidate model. In parameter estimation, the space-time extension of Hannan-Rissanen algorithm is implemented to generate the starting point for optimization process. For the first di— agnostic test, every parameter‘s statistical significance is evaluated, in order to decide whether the current numbers of maximum temporal and spatial orders need to be re- duced, increased, or neither. Finally, the residuals' STACF/STPACF are calculated to ensure that no significant spatial / temporal correlation exists in the residuals. If there are correlations, another iteration is needed for further refining the candidate model. 3.6 Cross Correlations and Space-Time Regres- sion Models Cross correlation is a standard method for estimating the degree to which two series are correlated. It is robust to noise and used widely in many disciplines for fez- ture/signal detection and image/data processing. To further extend the usefulness of the research, cross correlation analyses among the two space—time. variables are conducted. Cross correlation analysis is in‘iplemented in IEAST and also applied to the analysis of the epidemics of West Nile virus. This is especially useful when inves- tigating space-time relation between two variables. It does not necessarily give us the direction of causation flow, but, at the least, we can gain more information of how (and to what degree) variables are similar to each other. 47 Suppose we have two (correlated) space-time variables Z and Y. Z is dependent on an independent variable Y based on the Space-Time Regression model [40] C d z, = Z: w,,w<’>Y,..,, + 5,, (3.10) k=1 1 :0 where Z, is expressed as a weighted sum of Y, lagged in space and time, with added space-time noise 5,. In simple words, the variable Y is the cause. of the variable Z. The cross correlation of Z and Y is defined in [40] and transformed into a format as below E[Z;W(-“)Y,_k] x/Elzzztl - E[Y.)'W<5>Y.1' pz.y(s. k) = (3-17) where s is the spatial lag and kt is the temporal lag between the two series. In the cross correlation function, the two variables (Z and Y) are ordered, that is, pz‘y(8, 1:.) # pyz(8, k). [)z‘y(8, k) gives the cross correlation of the major variable Z and the. secondary variable Y, and shows the correlation values for various spatial and temporal lags (Z leading Y for some positive temporal lag). The differences or lags in space and time of the peak (of the cross correlation function) from zero are the distances that the two variables separate in space and time. In the analysis conducted for the spread of West Nile virus, these lags indicate the time delay and any spatial lags between the human cases and the dead crows findings. The cross correlation is based on an assumption of space-time regression model (Equation 3.16). The linear relationship between variable Z and Y is shown in Figure 3.4. It is evident that the 48 autocorrelation of Y also contributes to the cross correlation pz_y(s, k) because pzy = SpaceTim(2C01‘7'(Z,Y) = SpaceTimcCorr(STR(Y),Y) = Spa(:'cTimeC0rr(2::1 21;, wk1W(’)Y1_k. Y). The cross correlation can be spurious. \Nhen calculating the cross correlation of two variables, this will inflate the cross correlation values and further make it difficult to achieve a parsimonious model. One solution to this problem is to use pic-whitening to eliminate the autocorrelated components in both variables. To calculate [)Z.y(h’. k) with pre—whitening, we can first fit a STARMA model for the input data (Y) to reduce sufficiently the residuals to white noise; then, filter the input data with this model to obtain the white noise residuals, Y". Subsequently, we can filter the response data (Z) with the same STARMA model, and finally, calculate the cross correlation between the filtered response Z" and the filtered input Y“, that is, Pzaw (s, k). 3.7 The Development of IEAST The integrated statistical computing environment, IEAST (Integrated Environment for Analyzing STARMA models), was designed for analyzing and modeling stochastic processes based on general Space-Time AutoRegressive Moving Average (STARMA) models in two-dimensional gridded space. The current version only supports spatial dimension N x N, that is, where the numbers of cells in the X and Y directions are equal. IEAST takes data files (.dat), spatial weight files (.wet), program files 49 (.pgm), and spatial/temporal configurations (in the Setup menu) as inputs. then according to program or user selection, outputs plots, text numerical analysis results, and model equations for STARMA analysis (see Figure 3.5). This system was developed on RedHat Linux 9.0 using GNU OCTAVE v2.1.40 as the programming language. which is a high-level language primarily intended for numerical matrix computations. OCTAVE provides a convenient command line inter- face for numerically solving linear and nonlinear problems, and for performing other numerical Operations. OCTAVE is largely compatible with the commercial software MATLAB® and free of charge. Because of source code compatibility. IEAST can be easily run under Windows/Unix/MacOS without modifications. To provide flexibility, there are two user interfaces provided in IEAST: the menu- driven mode and the interpreter mode (i.e. IEAST programming mode). In the menu-driven mode, users can carry out the modeling or analyzing procedure by selecting among a series of hierarchical menu commands. Users must control all of the flow of modeling procedures by themselves. This is convenient, but inefficient in terms of time spent. With the programming mode, a simple and powerful IEAST programming language is provided. Highly integrated instructions are provided for users to compose STARMA modeling or space-time analysis procedures. Users can design an efficient and autonomous modeling procedure for specific applications by simply combining instructions of core functions and flow controls. It is commonly possible to implement an entire STARMA modeling procedure in less than 20 lines of code. This progrannnability is especially useful when the iterative. procerhuc is necessary for analysis (as in the modeling procedure shown in Figure 3.3). F igure 3.3 is a suggested general modeling procedure for using IEAST. This system is designed for two-dimensional STARMA analysis. However, by carefully designing the spatial weight matrices, this system can be adapted to one- dimensional or even three-dimensional systems. If the computer memory and com- puting power are sufficient, there are no limitations to the spatial dimensions of the dataset that can be explored. However, if the spatial dimension is too small, boundary effects would be significant and should be taken into consideration. After trading-off the computing power and practical needs, the largest spatial and temporal orders of STARMA models in this implementation are limited. IEAST consists of three major parts: menu-driven user interface, interpreter. and core functions. Figure 3.6 illustrates the logical structure of IEAST. The core functions provide a complete set of functionalities for STARMA modeling and space- time analysis, and they are also the basis for the menu—driven interface and interpreter. The core functions are categorized as below: 0 System setup: These functions are for configuring a given space—time system, and include loading space-time data, importing datasets from various formats, generating datasets from simulations, loading/setting up spatial correlation structures. and 2D/ 3D visualization of space-time datasets in time/ frequency domains. 0 Preprocessing: To generalize STARMA models, there are provided various transformations (or preprocessing) for datasets to be transformed into station- ary ones before analysis. These functions include de-seasonalizing, de—trending. differencing, subsequencing/resampling, smoothing, recovering missing data, pre—whitening with STARMA, and the Box-Cox transformation. 0 Correlation Analysis: Space—time auto-correlation (single variable) and cross- correlation (between two variables) analyses, and plotting ability are pr(.)vided. 0 Model Identification: This contains functions for identifying model type and maximum spatial/ temporal orders based on the behaviors of STACF and STPACF. 0 Parameter Estimation: For linear processes (STAR), this function simply pro- vides best linear estimators for model parameters. For nonlinear processes (STMA and mixed), two steps are needed. The first step calculates a start- ing point using the space—time extended Hannan-Rissanen algorithm. Then this point is used as an initial point for the multi-variate nonlinear optimization process for convergence of the estimates of parameters of models (lVIarquai'(lt‘.’s algorithm). In addition, several mechanisms for estimating spatial correlation structures based on a given dataset and a model (with known type and orders) are provided. 0 Diagnostic Analysis: Diagnostic functions provide various criteria for evaluating the adequacy of the obtained model. In our system, the residuals STACF and STPACF, variance analysis, AICC / BIC , and parameter significance analysis are used. The functionalities of menu-driven interface for analysis are distributed in the menu hierarchy. In this interface, users conduct step-by-step set up and analyze. the. dataset through the menus. follow their own modeling procedure. The first two levels of the menu system are shown in Figure 3.7. The current version of IEAST (version 1.30.01) has the following features: 0 A true spatial and temporal analysis software 0 Full configurability - spatial and temporal correlation structures can be specif- ically designed by users 0 Interactive programming environment - a script interpreter with highly inte- grated instructions is provided to simplify programming efforts and to improve efficiency of analysis and modeling 0 Estimation of spatial weighting structure - The spatial weighting structure can be retrieved and estimated from given datasets. o Flexibility of defining spatial structure - a three-step definition for spatial struc- ture makes IEAST very flexible in defining spatial structures, i.e. SOD ——+ 5 RM —+ S WM 0 Improved estimation algorithm - the introduction of space-time extension of Hannan-Rissanen algoritlml greatly reduces the estimation time and possil‘nlity of converging to wrong results 0 Improved diagnostic measures - the space-time AIC / BIC and statistical signif- icance for each parameter provide good criteria for the selection of model type and orders, which is especially useful in the prograimning environment 0 Cross correlation analysis - space—time cross correlation between two variables is provided 0 2D / 3D plotting abilities - both data and analysis results can be plotted in two- or three—dimensional graphs For further information and usage, please refer to the IEAST manual, which can be downloaded from the website http://fried.for.msu.edu/“ieast. Space-Time Autocorrelation (STACF) ‘ I I I I I I r Spatial lag 1 --.x--_ Spatial lag 2 ......,., a 9 ,_; g a). 2 ll ; x- i ‘5 i .C .‘ 0: 5 -o.5 - ~- ~ 2 4 6 8 10 12 14 16 Temporal lag Space-Time Partial Autocorrelation (STPACF) l l r I If f I I Spatial lag o —+—— Spatial lag 1 --.x--- Spatial lag 2 ....*_, Phi(s:T-lag, l) L i i 1 I .5 P b 2 4 6 8 1 0 1 2 14 16 Temporal lag Figure 3.1. STACF and STPACF for a simulated pure STAR dataset with maximum temporal order 2 and maximum spatial order 1: Z, = 0.52,-1 + 0.3W“)Z,_1 + 0.152;-2 + 51. CH C"! Space-Time Autocorrelation (STACF) 1 i T I I I r . I I 1 ' Spatial lag o _.,__ Spatial lag 1 -----.x Spatial lag 2 -......* Spatial lag 3 ----a ------ ES .‘9 rn 0' II 4‘. ES .2 a: -0.5 -- d -1 i i , i i i ,5 4 6 8 1 0 1 2 14 1 6 Temporal lag Space-Time Partial Autocorrelation (STPACF) 1 i i ! I i 1 . r i J .» Spatial lag 0 ._._ Spatial lag 1 ---)(--- Spatial lagz -.......,., Spatial lag 3 we“- 3 9 3 E o. ‘0.5 '- --1 _1 i 1 i 1 i i i 2 4 6 8 10 12 14 16 Temporal lag Figure 3.2. STACF and STPACF for a pure STMA dataset with maximum temporal order 1 and maximum spatial order 1: Z, = e, — (—0.6)€,_1 — (—0.4)W(1)e,_,. Dataset input Pre—processing (dc—season. de-trend,...) Space-time autocorrelations (STACF/STPACF) tod eci e Iterative! tcstin AICC/BIC to dccidct ' model type based on No y g wt oooooooooooooooooooooooooooooooooooooooooooooooooooooooo lterativcly set model type as STAR. STMA, STARMA based on STACF/STPACF 7‘ [ DCCidC S/I‘ orders J [Set to max possible S/I' orders] [ Decide model type [ Pre/Estimation based on STACF/STPACF ‘ [ Pre/Estimation 3 Calculate space-time AICC/BIC Dec or inc Statistical Significance max S/I' orders analysis of parameters ll Need to doc/inc max S/T orders? . , Choose the best model type [Resrdual s STACF/STPACF ] 5 based on AICC/BIC analysrs : ---------------------------------------------------------- Yes Exist significmt Sfl‘ correlations? Figure 3.3. The refined STARMA modeling procedure. The area enclosed by the dotted line is for the case where model type cannot be decided from the STACF/STPACF. AICC/ BIC of every possible model type is evaluated and com— pared, in order to select the best one. Space—Time Regression Model Z=S T R( Y) 1 Random noise Figure 3.4. The linear relationship between space-time variable Z and independent variable Y under the assumption of space-time. regression (STR) model. Dataset file l Weight Matrices Spatial/Temporal Programs (.dat) file (.wet) configurations (.pgm) Plots Text numerical results Model equation 1 Figure 3.5. Inputs and outputs of IEAST. III-“I013": mama?» .Hnmummogi" denies?” Eamon . . 93:890. mgfimmom, mace macaw , ” 255.3 manganese”, . 55mg: , gfiwam m QZCIOOH>5»;me Q Ema—momsob hmmaamcos gh0>=\\wqma Qhwnomoaaomlw ~=8€88~ jmEfiaIaBaJN w 350—“ JN Lennon—man JN\\ _... J m82m~83J Nm m <. manna ( @389 L (2:08 38W L CamacmomsoaL f :52. L (Ems—momsnaLf mz,.2 —.11W<1>zH +.10w<2>z._, —.07W(2’Z,_2 +.02W<2)Z,_3 —.09W<3>z,_-l —.04W(3)Z,_2 —-.02W<3>z,_3 +.04W<4>z,_1 —.11W<4>z,_2 —.03w<4>z,_3 +E,. 71 Table 4.1. The estimates and significance levels of the model 1.)arameters (STAR(MaxTr—B, MaxS=4)). Estimates of parameters Spatial lag Temporal lag 8:0 8:1 8:2 8:3 8:4 T=1 .26 .36 .10 -.09 .04 T=2 .04 -.18 -.07 —.04 -.11 T=3 -.02 -.11 .02 —.02 —.03 Significance levels of parameters Spatial lag Temporal lag 8:0 8:1 S22 S23 S24 T=1 fl Jill $0 .100 .400 T=2 .040 m .040 .300 £10 T=3 .400 .01 .250 .900 .600 Two diagnostic checks were conducted for this analysis: significance test for the. parameters and correlation analysis for the residuals. The significance test. results are shown in Table 4.1. Several parameters at relatively high temporal or spatial lags were statistically insignificant. However, technically the maximum temporal and spatial orders from the candidate STAR(maxT=3, maxS=4) should not be reduced, because highly significant parameters were found at S=4/T=2 (spatial lag 4 and temporal lag 2) and S=1/T=3 (spatial lag 1 and temporal lag 3). The STACF/STPACF of residuals show no significant correlation among residuals (Figure 4.8). Therefore, the candidate model represented the underlying process of the DCD adequately. as it is well-fitted by the data and it is parsimonious. In addition to the 10 x 10 cell configuration, the same area was divided into various arrays of cells up to 20 x 20, for analyses in different configurations. Based on the same spatial correlation structure (i.e. same SOM and SRM) and de-trending methods, these modeling results were consistent. The maximum spatial order and the parameters of the model differed when different arrays of cells were used to model the same data as is expected, because smaller cell size is associated with higher maximum spatial order needed for the model. However, the physical meaning of the modeling results does not change (as shown in Table 4.2 and 4.3). Three analyses using various arrays of cells: 10 x 10, 16 x 16, and 20 x 20 were performed. Only parameters that met the 0.01 significance level were entered into the final models. Based on the significance level, the maximum spatial orders for these three cases can be determined as fourth, sixth, and seventh, respectively. The maximum spatial order can be transformed into distance in the following way. For configuration 10 x 10, the cell size is 3.16 x 2.58 miles. According to the spatial lag definition in Equation 4.1, the spatial lag four is two cells away in X-direction and one cell away in Y-direction from the central cell (lag zero). Thus, the physical distance in the configuration 10 x 10 for spatial lag four can be calculated as the following. physical distance = \/(2 - 3.16)2 + (1 - 2.58)2 = 6.83 miles. After transformations, the equivalent physical distances (measured between cell cen- ters) for these maximum spatial orders were found to be quite uniform (6.83. 6.76. and 6.45 miles, respectively). This shows that the spatially correlated areas from these modeling results are consistently about 6.4 miles in radius, regardless of the. cell size used. Furthermore, additional analyses were conducted over a wide range of still other 73 Table 4.2. The maximum spatial lags and the corresponding distances between cell midpoints for that lag, for three models using different cell numbers (larger numl')ers, smaller areas). Estimated Equivalent Cell numbers max S order real distances 10x10 4 6.83 miles 16x16 6 6.76 miles 20x20 7 6.45 miles configurations, including: various datasets of sub-sampling of various parts of the entire area; and various other cell configurations and SOMs and SRMs. For each of these, the refined iterative modeling procedure and the de-trending method were used. In all cases (results not shown), the model that was converged to was consistent with the model reported. The results show that the incidence of reported dead crows per week fitted a. STAR model. A high degree of space-time autoregression was found, and no evidence of moving average influences. The latter result indicates that although there was both stochasticity and statistical noise in the observed values (the estimated variance is 1.451), these sources of variation were not directly or immediately shared spatially or temporally among cells. Rather, they were shared only through the autm'egressive structure attributable to spatial spreading of WNV. Diagnostically, the envelope of the STPACF cuts to zero after temporal lag three. whereas the STACF tails off with increasing time lag (Figure 4.7). The robustness of this diagnostic check, especially so since the process is STAR, was also indicated in the extensive irwestigations us- ing computer simulations. Moreover, the residuals are those expected from random (white) noise. There do not appear to be any particular reasons to expect STMA 74 spatially-sl'iared stochastic inputs or shared random errors in the DCD. The results also show that the process underlying the DCD has low spatial and temporal orders. Moreover, the autoregressive coefficients show a pattern that gen- erally fits those expected for most biological processes. For example, consider the 10 x 10 array of cells and the spatial order matrix (SOM) of Equation 4.1 and the spatial relation matrix (SRM) of Equation 4.2. As could be expected, the autore- gressive parameters are positive and large for close spatial proximities. i.e. small distances, and drop in value as distance increases (Table 4.1). Exceptions involving fairly large negative parameters for short distances at time lags two and three are discussed below. The results also showed that the temporal order of the process is nearly lV’larkovian, since large and significant parameters are observed only for time lags of one and two weeks. It has been argued that the Markovian property should be common in biological processes [24]. In addition, the prominence of one- and two-week lagged effects is consistent with the short interval between infection and death [35], and with empirical evidence of intense localized epizootics [20, 52]. Using the same SOM and SRM, the best fitted model is STAR(maxT=3, maxS=4), and completely consistent results were found using various other arrays of cells. The maximum spatial lags consistently gave maximum distances of ca. 6.4 miles, in terms of distances between cell-centers. I also examined models with various other SOMs and SRMs, and again obtained completely consistent results. H(,)wever. as noted in the results and further explained below, there was evidence. that the. mod— els could be reduced further in size; the model might be reduced to STAR(maxT=2. maxS=2), consistent with the array of positive autoregressive effects. 75 Table 4.3. The maximum spatial lags for positive autoregressive effects and the corresponding distances between cell midpoints for that lag, for three models using different cell numbers. Estimated Equivalent Cell numbers max S order real distances 10x10 2 4.07 miles 16x16 3 3.96 miles 20x20 4 3.41 miles Focusing on positive autoregressive effects, it appears that the DCD operates over rather smaller distances, with maximum distances of direct effects of locations on one another of ca. four miles or less. However, the distances over which the virus load is primarily moving (spreading) may still be much smaller. There is a remarkable reduction of autoregressive effects for spatial lag 2 (0.10) compared to spatial lag 1 (0.36) for the 10 x 10 array of cells and SOM of Equation 4.1 and SRM of Equation 4.2. The main difference between the two is the length of the boundaries between cells with various spatial lags. This strongly suggests that the length of boundary largely determines the rates of spread of WNV among cells. If that is the case, it would indicate that \NNV directly spreads mostly over very short distances across cell borders, less than a mile or possibly even a few hundred meters. Crows can fly long distances, so the scale over which the DCD show large positive autoregressive effects may seem surprisingly small. However, the mosquito vector is active primarily during hours when the crows are roosting, and crows are territorial and show high fidelity to roosting sites [57]. Thus, crows may largely be being infected and transmitting the virus at their roosting sites, and thus apparently not spatially spreading the disease, although they would amplify the viral load locally. In addition. 76 the distances of flights between where dead crews are found and their roosting sites are superimposed on this process [58]. Although such movements do not spread the disease, they could influence the short term autoregression of DCD. While crows are a very good indicator of viral loads, they may not be largely responsible for spreading the disease spatially. Other animals or even the mosquito vector itself may be more responsible for the local spatial spread per se. Two points should be reinforced. First, since crows are one of the most ii‘ifected animals, they could be very important in amplifying the disease locally (say within a cell). Second, the results do not mean that the crows do not on occasion spread WNV over long distances, and perhaps even start a new local epidemic. It is unlikely that such low frequency events would be detected in the analyses. Nevertheless, the results do suggest that there is high probability that local epidemics could be contained. for example, through a concentrated eradication effort, by spraying mosquito repellant at mosquito breeding ponds, etc. One specific feature of the autoregressive structure that deserves special exami- nation is the rather large negative autoregressive parameters at time lag 2 for some. spatial lags. This was observed in all models, and appears to be a real effect. It may be explained by the depletion of the crow population. Depletion of the crow popula- tion was extreme: during the course of the WNV season in 2002, up to 70 % of crews in Detroit died from 'WNV (K. Signs and J. Patterson, unpublished). Hence, for example, if a cell had unusually high DCD for a given time period (week), especially during the height of the WNV season, a significant portion of the crow population in that cell would have died, leaving fewer to be infected and die one to two weeks 77 later. Thus, the results suggest that there are actually two DCD autoregressive pro- cesses. The first is dependent on the local and nearby W NV loads. and the second is a weaker echo effect caused by depletion. The mixture of the two may make temporal lag one DCD autoregressive parameters all positive, but not all of the temporal lag two parameters. 4.6 Results for the Analysis and Modeling of the Human Case Dataset The analyzing procedure followed for human case data was similar to that for DCD. The same spatial correlation structure (i.e. same SOM and SRM as in Equation 4.1 and 4.2) was used for spatial dimension 10 x 10, but a fourth order polynomial trend surface combined with the temporal trend (weekly average over space.) (as shown in Figure 4.9) and the grand mean were removed from the human case data before STARMA modeling. The autocorrelations for the de-trended human case data. are shown in Figure 4.10. It is not surprising that the autocorrelations fluctuate considerably since the number of cases of human case data is much smaller than that of DCD. However, the behavior of STACF/STPACF still revealed some clues about model type. The curves in the STACF gradually decay to zero as temporal lag increases. On the contrary, the curves in the STPACF drop quickly and cross axis at a short temporal lag, and then fluctuate and gradually diminish at higher temporal lags. This suggests that the pure STAR model is a good candidate for the de-trended data. 78 Regarding the maximum orders of the candidate model, we can start with by using six as the maximum for both spatial and temporal lags, because the STACF showed that there’s no significant correlation for spatial and temporal lags greater than six. Although a STAR model with maximum spatial and temporal order six may be too complex for the data, we start from the candidate model STAR(lVIaxT=6,l\IIaxS=6) and then iteratively refine the model to achieve a parsimonious one. Given the candidate model STAR(MaxT26,l\‘IaxS—-—6), the model parameter esti- mates and the significance level for each parameter can be found in the stages of esti- mation and diagnostic checking, respectively (Table 4.4). As in the analysis of DCD, if significance level is higher than 0.01, the corresponding parameter is classified as insignificant and then is removed from the candidate model. Once all parameters for a temporal lag (or a spatial lag) are removed, the maximum temporal or spatial order may be reduced. In the significance analysis listed in Table 4.4, all significance levels less than equal to 0.01 are underlined. The results show that the maximum temporal order four and spatial order six are sufficient for these significant parameters. Thus. the new candidate model is a reduced STAR(MaxT=4,MaxS=6). Next, based on the new candidate model, the new parameter estimates and corresponding significance levels are calculated as in Table 4.5. Note that the values of retained parameters are very similar to those estimated using STAR(h=IaxT=6,l\/IaxS=6). Based on significance tests, it is not possible to further reduce the maximum orders because the parameter with temporal lag four and spatial lag six is significant (significance level S 0.01). Substituting the parameter estimates into the model equation for STAR(MaxT=4,MaxS=6), a mathematical equation for 79 Table 4.4. Initial estimates and significance levels of the model parameters (STAR(MaxT=6, MaxS=6)) for the human dataset. Estimates of parameters Spatial lag Temporal lag 8:0 8:1 8:2 8:3 S=4 8:5 8:6 T=1 .26 .06 —.10 —.29 —.05 -.30 -.06 .12 .27 .13 -.12 -.11 -.22 -.11 .07 .10 —.15 .05 -.00 .06 -.01 .04 —.17 -.07 -.02 .16 .25 .11 -.01 -.10 -.04 .10 -.06 .11 .06 -.04 .08 .09 .03 -.03 -.19 ~09 Hawk-la acne—cow Significance levels of parameters Spatial lag Temporal lag 820 8:1 822 8:3 S=4 8:5 8:6 T=1 £01 .700 .040 w .500 .100 .150 T 2 m .100 .040 .150 .300 .020 .020 T 3 .900 .400 .040 .150 .900 .400 .900 T24 .900 ,0_1_Q .900 .900 .020 010 ._0_m T 5 T 6 .500 .300 .900 .080 .500 .300 .150 .900 .250 .150 .900 .600 .040 .040 80 Table 4.5. Final estimates and significance levels of the model parameters (STAR(MaxT24, MaxS=6)) for the human dataset. Estimates of parameters Spatial lag Temporal lag 8:0 8:1 8:2 8:3 S=4 S=5 S=6 T=1 .26 .06 -.11 -.30 -.07 -.30 -.06 T=2 .10 .27 .12 -.11 —.09 -.22 —.11 T=3 .05 .07 -.15 .06 .01 .06 -.01 T=4 .03 -.14 -.07 .03 .16 .21 .10 Significance levels of parameters Spatial lag Temporal lag 8:0 8:1 8:2 8:3 S=4 S=5 S=6 T=1 £01 .700 .040 m1 .400 M . 150 T=2 .020 . 100 .060 . 150 .400 .040 .020 T=3 .900 .900 .040 . 150 .700 .400 .900 T=4 .900 .040 .800 .400 .020 .020 ill—0 representing the underlying process of the human epidemics can be obtained as in Equation 4.5. The autocorrelation analysis for the residuals shows no large correlations among spatial lags zero to six and temporal lags one to four. However, there do exist some correlations in higher temporal and spatial lags. The correlations in higher spatial and temporal orders (which are not included in the. model) are not as significant as that in the lower orders (Table 4.5). Therefore, this model (Equation 4.5) can be accepted as the final model for representing the process of human epidemics. In summary, the human epidemic in Detroit Metro area is highly autocorrelated within spatial lag six (equivalent to 12.57-12.86 miles away in physical distance) and four weeks time period. 81 z,= .QGZH +.10z,__2 +.osz,_3 +.03z,_, +.06W<1>z,_1 +.27'sz,_2 +.07W<'”z,_3 —.14W<‘>z,-,, —.11W<2>z,_1 +.12W<2)z,_.2 —.15W<2>z,_3 —-.07W<2>z,._.1 —.30W<3>z,_1 —.11W<3>z,_‘2 +.06w<3>z,_3 +.03w<3>z,_4 (4' —.07W<4>z,_l —.09W<4>z,__2 +.01W<4>z,_3 +.16w<4>z,_4 H —.30W<5)Z,_1 —.06W(6th_1 —-.22W<5>z,.2 —.11W(6)Zt_2 +.06W<5>Z,_3 —.01W(6)Zt_3 +.21W(5)Z,_., +.10W(")Z,_4 +5! There are several points needed to be noted. The area defined for analysis for DCD is different from that for human case data. Nonetheless, this would not affect the final results because the spatial correlation structure used for both datasets is directionless or isotropic. The highly infective area (up to spatial lag 6. or about 12 miles) for human is larger than that for dead crows (up to spatial lag 4, or about 6 miles). This infers that the disease spreading of human cases is, spatially, faster than that of crow cases. It can be explained by the mobility and active area of human is much larger than that of crows. Because of this, the observation error in human case dataset is also larger than that in DCD. From the parameter estimates, the significant parameters for DCD tend to be gathered around lower spatial orders, and tend to be dispersed in a wide area for human case data. In short spatial lags (one or larger) and temporal lags(two or larger), especially for spatial lag one at temporal two or higher, the parameters for 82 human case data is positive and significantly above zero. Contrarily, the parameters for DCD are negative and significantly below zero. Because human cases should not have evident depletion effect as compared to crow cases, this further justifies that the depletion effect among crow is the major reason causing these negative parameters around low spatial orders for DCD. 4.7 Summary of Analyses for the Dead Crow Dataset and the Human Case Dataset As needed for the first step of the STARMA analyses, the space-time autocorrelation function (STACF) and the space-time partial autocorrelation function (STPACF) are shown in Figure 4.7 for the DCD and in Figure 4.10 for the human case data. For both datasets, the STACF values are very large at short spatial lags (distances) and temporal lags, and they tail-off as either temporal or spatial lag increases. The STPACFs of both become and stay near-zero (cut-off) at small spatial and temporal lags. These fundamental aspects of the both dataset are diagnostic in slmwing that both DCD and human cases behaved as a STAR process, and possessed no shared stochastic inputs. In other words, each space—time series behaved as a purely autore- gressive process with added white noise (spatially and temporally random stochastic inputs). Moreover, both datasets have small spatial and temporal orders, the maxi— mum possible values of which are determined by the lags at which the STPACF cuts to zero. 83 For the DCD, based on inspection of the STPACF, the maximum spatial lag is four and the maximum temporal lag is three. Estimates of the space-time autoregression coefficients for various spatial and temporal lags for the STAR(maxT23, maxS=4) structure are shown in Table 4.1. Large and highly significant positive coefficients are only observed for (T21, S20), (T21, S21), (T22, S20). This indicates that most of the spread operates between nearest neighbor cells, and over a two week time period. It is also noteworthy that the coefficient value is much larger for rooks move. (spatial lag S 21) than for bishops move (S22) nearest neighboring cells. Similarly, the STPACF for the human case data indicate that the maximum spatial lag for positive autoregression is six and the maximum temporal lag is four. Estimates of the space-time autoregression coefficients for various spatial and temporal lags for the STAR(maxT24, maxS26) model are shown in Table 4.5. Diagnostic checking for model adequacy is conducted not only by the significance test of parameters, but also by inspection of the STACF/STPACF of the residuals. The graphs of the STACF/STPACF of the residuals for the de-trended DCD, as shown in Figure 4.8 are flat throughout, closely resembling that of the random white noise. This is consistent with the fitted STAR model having captured all of the autoregression in the de-trended data [16]. The same was observed in the STACF and STPACF of the residuals of the fitted model for the human case data. As noted before, I further examined both datasets with other model construc- tions. I varied cell sizes and used sub-areas of the entire dataset, and in all cases, for both datasets, the best fitted models were purely STAR models, and these had fully consistent distances and time scales. I also used different polynomials for spatial 84 de-trending. Again, the results were fully consistent, except that very large-order polynomials generally tended show reductions in autoregressive coefficients. This is expected, since over-de-trending generally removes some autoregression or autocorre- lation [8]. 4.8 Results for the Cross Analysis Between Hu- man Case and Dead Crow Datasets Dead crews are sentinels of WNV viral activity generally and locally, and there, is a burgeoning effort to use dead crow reports as a local signal of an outbreak (Theophilides et a1 [52] with their Knox statistic and Watson et a1 [56] with their geostatistical analysis). However, there were no studies that considered spatial and temporal interactions simultaneously, nor that presented a cross-correlation analysis between dead crow reports and human cases. In addition, there have been no studies that incorporated specific spatial units and project the spread of WNV. Therefore, a few question arose: how predictive in space units and time units are DCD for human cases? What is the predictive power? If the two processes (dead crow epizootic curve and human case epidemic curve) are demonstrably space-time regressive, then a model can be developed in which spatial and temporal lags with specific model parameters are provided, showing the ” strength” of the relationships at defined space-time units. In this study, the cross correlations between the DCD and the human case dataset are investigated using the method mentioned in Section 3.6. The results provide a. 85 foundation for creating a space-time regression in which dead crows data are used as an indicator variable to derive the space-time dynamics of human epidemics. A plot of 1,753 dead crow sightings was successfully classified by latitude and longitude and by date of sighting, spanning 28 weeks. The plot of 385 human cases of the same rectangular area covered as in previous dead crow sightings was classified and spanned 14 weeks. These data points and the area for the cross correlation analysis are shown in Figure 4.12 and Figure 4.13. These data were plotted within a nearly square area of 39.17 miles x 39.23 miles over a map of the 111etropolitan region of Detroit, iV'Iichigan. This area was again divided into a grid of 10 x 10 cells, each cell with 3.92 miles x 3.92 miles dimensions. Because the human case data spanned from Julian week 30 to week 43 and DCD from week 16 to week 43, zeros was applied to human case data for the cells in time periods in which there were no cases reported. Thus, in space and time the DCD and human case data cover the same time period and physical area. Accordingly, the data set was of spatial dimensions 10 x 10 and had 28 time periods. We also assumed that the spatial correlation structure of the data. was isotropic and invariant with time. A higher order (up to ten) of spatial correlation structure. is used for the cross analysis (the SOM and SRM are shown in Equation 4.6 and 4.7). The method in section 3.6 was used for calculating cross correlations, and no spatial and temporal de—trending was used. 86 SOM 2 .083 .050 SR!” 2 .050 .050 .083 .083 .050 .083 .083 .083 .083 .083 .050 .083 .050 .083 .125 .063 .063 .063 .125 .083 .050 01 Cr! CB 00 CD CI! (30003 .083 .083 .125 .063 .125 .250 .125 .063 .125 .083 .083 p-i O>OOO C31 0060142006:- 10 .050 .083 .063 .125 .250 .250 .250 .125 .063 .083 .050 i—IMJACHKIQD iii-[0 CI! 9 DOV-JO CH ON] .050 .083 .063 .250 .000 .250 .250 .063 .083 .050 [0.330% “>53 CI! 9 C" C) 00 430033 .050 .083 .063 .125 .250 .250 .250 .125 .063 .083 .050 10 290003010101me CDOONNNWO .083 .083 .125 .063 .125 .250 .125 .063 .125 .083 .083 .050 .083 .125 .063 .063 .063 .125 .083 .083 .050 .083 .083 .083 .083 .083 .050 .083 .083 .050 .050 .083 E 4.7) As shown in the analysis of results for the DCD and the human cases, the two space.- time process both behave as a STAR process and possessed no shared stochastic or statistical noise inputs. The space-time cross-correlations between the DCD and the 87 human case data are shown in Figure 4.14. The spatial lags span from 0 to 10, and the temporal lags from -12 to +12. The temporal lags in the cross correlation represent the number of weeks in which human cases reported lagged after dead crow cases reported. The cross correlation coefficients were very large (as high as 0.7). Across all spatial lags the peaks were such that DCD preceded the human cases by three weeks. and these cross correlation curves are perfectly aligned. Although this result does not provide a proof of causation flow between dead crows and human cases, the result do support the use of dead crow data as a space-time indicator variable for forecasting human epidemics. Moreover, the correlations drop smoothly and symmetrically for temporal lags deviating from the three week time lag and spatial lags increasing from 0 to 10. Generally the correlations also decrease smoothly with spatial distance. An exception is the slightly greater values for spatial lag one than spatial lag zero. This may infer that crows died in adjacent cells is even more correlated with the incidence of human cases than the crows that died in the same cell. However. this outcome may simply result from sampling effects. At least, we know that the cross correlations between human cases and dead crows for spatial lag one are as large. as those for spatial zero. When temporal lag decreases to -8 or below, the correlations between these two datasets are negligible (less than 0.1). This result indi :ates that the crow epidemic that occurred eight weeks earlier no longer correlates with current human epidemic in space and time. When spatial lag increases up to 10 (equivalent to 21.1 miles), the cross correlations are reduced to as low as 0.2. As mentioned in section 3.6, the autocorrelation of the DCD can spuriously con- tribute to cross correlations. To eliminate this effect, the cross correlation is conducted 88 with pre—whitening. Both the DCD and human case data are ‘filtered' with the esti- mated crow model equation (Equation 4.4), and then calculate the cross correlation between the two filtered space—time series as usual. The result is shown in Figure 4.15. Generally, the cross correlation is lower, but there is no significant difference between Figure 4.14 and Figure 4.15 in terms of pattern. This provides compelling evidence that the high cross correlation between the human and the dead crow data is causal and not spurious. 4.9 Discussion The STARMA framework of statistical space-time modeling appears to hold consid— erable promise for characterizing the space—time autoregressive structure and corre- lations for many biological processes, including many variables associated with the, spread of zoonotic infectious diseases. The application of STARMA to dead crow data (DCD) as indicators of the environmental load of the W NV in metropolitan Detroit showed that DCD follows a STAR process having low spatial and temporal orders. This outcome is consistent over various definitions of spatial lag structures. The results indicate that the incidence of dead crow does not spread spatially very far very fast. Indeed, a focus on positive autoregressive effects and the contrast be.- tween rooks—move and bishops-move nearest-neighboring cells suggests that most of the spread is determined by the length of shared boundaries. Biological C(msiderations suggest that although crows are a likely cause of local amplification of WNV. they are not spatially spreading the virus very far very fast. It is possible that other hosts 89 or even the mosquito vector itself may be more responsible for the spatial spread. In addition, there appear to be significant. effects of depletion of the crow population on the space-time dead crow incidence. The autocorrelation results for modeling the human case data appear more er- ratic than that for DCD (as shown in the STACF/STPACF for human case data in Figure 4.10), however it follows the same. pattern as the DCD. It should also be kept in mind that the number of human cases is much smaller than that for DCD. Such paucity in the number of samples can introduce statistical noise during analysis. In space-time modeling, the condition of abundant samples in time and in space is critical. It is important to discuss the de-trending over time and space. The fact that we de-trended temporally over the season means that the autoregressive parameters are not rates of spread, but rather relative rates given the overall increase or decrease in WNV. The former might be recovered by adding multipliers to the autoregressive (the variance of the errors term 5:, might also change over time) structures, depending on the exact week. More interesting is the rather striking spatial non-stationarity. Since these have been removed by polynomial regression, the analyses reveals nothing about why they occurred. Finally. in the cross analysis, the two datasets revealed very high correlation which reaches the maximum values for each spatial lag at temporal lag -3. That is, the num- ber of crows dying in neighboring cells three weeks earlier has the highest ('rorrelation with the number of human infections. This result supports the assumption in the epi— demiological literature that crow epidemic is leading human epidemic by three weeks 90 and it strengthens the basis for using dead crow data as an indicator for forecasting human epidemics. 91 Iirivt Macomb : O OaMand Ann Arbor Figure 4.1. Dead crows reported in the Detroit Metro area. The rectangular area enclosed by thick line is the area retained after truncation. 92 Week 1 Week 2 Week 3 Week 4 Week 5 Week 10 Week 12 Week 18 Week 19 Week 20 Week 24 Week 25 20 15 10 5 0 Week 26 Week 27 Week 28 Figure 4.2. Dead crow dataset in the Detroit Metro area with 10 x 10 cells in space and spanning 28 weeks in time. 93 Figure 4.3. Human cases reported in the Detroit Metro area. The rectangular area enclosed by thick line is the area retained after truncation. 94 Week 1 Week 2 Week 3 Week 6 Week 7 Week 8 Week 11 Week 12 Week 13 Figure 4.4. Human cases in the Detroit Metro area with 10x10 cells in space and spanning 13 weeks in time. Figure 4.5. The spatial trends (averaging over time) of dead crow (upper) and human case data. The areas cover by these two datasets are not aligned geographically. 96 Trend Surface using PloyRegr - 4th order Z(x.y.t) OPPPP‘fff b thm mam I O Temporal Trend of Z (Average over space) 4.5 | l l TemporaI'Trend —o— 3.5 - .. 2.5 I- ~ Z(t) 0.5 — _ time Figure 4.6. The spatial and temporal trends removed from the dead crow dataset. 97 Space-Time Autocorrelation (STACF) l l f F l r T Spatial Iago —+—— Spatial Iag1 ---X--- Spatial Iagz ------x J Spatial lags ~~~~~~ aw Spatial lag4 -1.-- Spatial lags mo.- Spatial lagG We M I \ 3 ; l . , . _ 0,s:T-lag) Rho(l.k - i i g 4 5 6 7 8 Temporal lag 9 1O Space-Time Partial Autocorrelation (STPACF) I l l T l l r Spatial lag o _._ Spatial lag 1 --.x--- Spatial lag 2 My: .- Spatial lag 3 ...... 3.... Spatial '39 4 ---'._. Spatial lag 5 --.o..- Spatial '89 6 -- "Q- .. Phi(s:T-lag, I) i L i Temporal lag 5 6 10 Figure 4.7. Space-time autocorrelation analysis of the de-trendod (lead crow dataset. 98 Residuals' Space-Time Autocorrelation (R-STACF) R-Rho(l,k=0,s:T-lag) r T l T r I I Spatial lag o .__.__ Spatial '39 1 ---x..-- Spatial lag 2 nun... Eflxmanaga ...... {J ...... Spatlai '39 4 -~—.._. Spatial lag 5 --_o.._ .....$9§Ii§i.i§9§..17."? .., ‘ -O.5 b “d " ‘ ' i i ; r 2 4 6 3 1O 12 14 15 Temporal lag Residualsv space-Time Partial Autocorrelation (Fl-STPACF) 1 l l I I I I . I S E E 5 Spatiallago __.__. Spatial lag1 ---x--- Spatial lagZ ------x Spatlallaga ...... a. ...... Spatial lag4 _.-..-. l 2 t i : fflXMaHag5 ~4am 0.5 L Spal'allagsg‘ g? 3 . E I q i ‘r -O.5 — _ ”q '1 ' l ' i i i i 2 4 6 3 10 12 14 16 Temporal lag Figure 4.8. Space-time autocorrelation analysis for the modeling residuals of dead crow dataset 99 Z(X.y.t) 2.5 2 1.5 1 0.5 0 -0.5 0.9 Trend Surface using PloyRegr - 4th order Temporal Trend of 2 (Average over space) 0.8 - 0.7 - 0.5 - Z(t) 0.3 - 0.2 - 0.1 — l l I T time 14 Figure 4.9. The spatial and temporal trends removed from the human case dataset. 100 Space-Time Autocorrelation (STACF) 1 I i i i i 7 Spatial iagc') _._ Spatial lag 1 ---x--- Spatial I39 2 "'*--. Spatial I39 3 ...... 3.... Spatial lag 4 -----—. Spatial lag 5 “9.- ....§Pa.tl§l,l§9.§.'.T.‘.,‘,.’..'T...._ Rho(l,k=0,s:T-lag) _ i— i l 1 2 3 4 5 6 7 8 9 1 0 Temporal lag I —A . ~- . I : : - , P - Space-Time Partial Autocorrelation (STPACF) 1 i i i I 1. i I Spatial Iago —-i— Spatial lag1 --—x--- Spatial lagz m-x- Spatial lag 3 ------ am Spatial lag4 ---n--- 5 3 3 3 3 ; Spatial lag 5 ---o-- 0 5 __l .. .. . . .. Spatiallage ‘- Phi(s:T-Iag, I) _ - . . .. . .r r . ....... ,. .. .........r...r. .......... - ., i L . s s 1 2 3 4 5 6 7 8 9 10 Temporal lag Figure 4.10. Space-time autocorrelation analysis for the de-trended human case dataset. 101 Residuals‘ Space—Time Autocorrelation (Ft—STACF) I 1 l l’ 1 l l ' . l ' Spatial lag 0 —i— Spatial Iagi ---x- Spatial lag 2 x Spatial lag 3 “or Spatial lag 4 — 1 - , Spatial lag 5 — —o - 05 _,. . , , Spatial lags o A a '0. 2 3'. ° .: : E’ i .C ‘F a: .05 — ............................ _ -1 1 I l 1 l I 1 1 2 3 4 5 6 7 e 9 10 Temporal lag Residuals' Space-Time Partial Autocorrelation (R-STPACF) 1 I I I r I T ' . I Spatial lag 0 ——+—— Spatial lagi -----x- Spatial Iag 2 max-- Spatial lags ma . Spatial lag 4 - -I — Spatial lag 5 — ~o — 05 e » ,, . , , . . . _ Spatial lag 833a”- 35’ L", E ‘¥ 0: -0.5 - fl .1 I 1 L l I l l l 1 2 3 4 5 6 7 8 9 10 Temporal lag Figure 4.11. Space—time autocorrelation analysis for the modeling residuals of human case dataset 102 Figure 4.12. Dead crows reported and area for cross analysis. 103 Figure 4.13. Human cases reported and area for cross analysis. 104 Rho(l,k=0,s;'r-|ag) @2258 2509.226 03mm 028.26: Exam-383 J. n . . ._ d magi. in C III: mega .8 d ...x... momma .8 m 55.x mama»: .8 w ...... m. ...... magi. So 5. iili. wage. .3 m 3-9- mumza _m© m .. :0. wage. .mou 5.»--. i 25$. .8 m -5»:-. mama». .8 m IT mama. .8 no 1-4--- .- ............................................................................................................................................................................................................ i Lo o .3338. .8 3 Figure 4.14. Space-time cross correlation between the human case data and the dead crow data. 105 Rho(l,k=0,s:T-lag) mxgmaama mumom-._..3m 08mm 026.30: “magma-x03 d a woman. in 0 lil momzmz .8 d -lx.-- mum—am. .8 m -...:*. mom—am. in w ...... m. ...... mafia. .8 a 5-9.-.. momma in m 5-9.- mama. .mom .9 . mum—gm. So u 5-»..- a @0922 50 m 5-9.--. magi. in m |.4Il momma in do ..ifl p m a _ a -8 -m o m a .3388. in Figure 4.15. Space-time cross correlation between the filtered human case data and the filtered dead crow data. 106 CHAPTER 5 Conclusion In this thesis, the application has demonstrated the effectiveness of STARMA mod— eling in analyzing and modeling space-time data. specifically in disease spreading. At the same time. it also provides a demonstration of using the refined modeling framework and related statistical tools/ algorithms. If the progrannning language in IEAST is used. we can even automate the whole process of analysis and modeling- and increase efficiency and flexibility of cmnputations. However, there are some points or limitations to be kept in mind when modeling space-time data. First, the number of data points in space and time need to be ’sufficiently’ large, or the random distribution (point pattern processes) of the given data can be quite different than the basic assumption of normal distribution in our STARMA model. In this research. several transformation methods (transforming to near normal distribution) have been applied to the DCD and human data to try to reduce the influence of the potential non-normality. Incidentally, after comparing the analysis results of transformed and non-transformed data, any differences are subtle. 107 This may result from the nature of autocorrelations. that is. the relative values or positions (in time or in space) are more important than the absolute values. Second. the STARMA model assumes that the value of the space—time variable for a given time and location is ’linearly’ correlated with its past values in time and in neighboring locations. Although in reality there is no purely linear process. most of them are. either close to or able to be transformed into linear process. If we know some process is highly non-linear. some methods have to be used to transform it into a linear one, or STARMA models are not appropriate for modeling such processes. Third. as mentioned in the basic assumptions, the correlation structure (in space and in time) of the process underlying the given data should be tune—invariant. In fact, no model can represent a process with constantly changing characteristics. Fourth. as implied in the context, STARMA modeling using IEAST does not allow missing observations. For cells without samples or observations, its value simply is regarded as zero (e.g. no cases for the WNV modeling). If missing observations is an important issue. the dataset needs to be interpolated or extrapolated in advance with some other temporal or spatial methods, such as Krigging. smoothing, etc. For future research, there are two directions for follow-up. First, anisotropic (directional) correlation structure can be an option while modeling. In my research. the spatial correlation structure (or weight matrices) was assumed and assigned to be uniform or directionless. Not every natural process is isotropic. It is especially interesting that an unknown spatial correlation structure can be estimated for a. given space-time data. For example. with some epidemics, there may exist a significant tendency in some specific direction. Once the spatial correlation structure is retrieved 108 from the given data. it can be very useful in forecasting the spatial spreading tendency. In IEAST. there is a functionality for estimating spatial correlation structure based on a given dataset. However. for most cases. the nmnber of parameters to be estimated for a correlation structure is too large and the computation needed intensive. involving too much computation to be applied to real data to get useful results. The work of developing algorithms or statistical tools to simplify this calculation is necessary. However. in general the model for forecasting human epidemic from dead crow data is very useful for prediction and forecast. 109 [ll [2] l3] [4] [10] [11] BIBLIOGRAPHY H. Akaike. Information theory and an extension of the maximum likelihood principle. In Petrov. B., F.. and C., editors. 2nd International Symposium on Information Theory. pages 267—281. Akademiai Kiado, Budapest. 1973. H. Akaike. Time series analysis and control through parametric models. In D. F. Findley. editor, Applied Time Series Analysis. pages 1-23. New York, 1978. Academic Press. L. A. Aroian. Time series in m—dimensions: autoregressive models. Communi- cations in Statistics - Simulation and Computation. BQ (5)2491’513, 1980. L. A. Aroian. Time series in m-dimensions definition, problems and prospects. Communications in Statistics - Simulation and Computation, B9 (5):453 465. 1980. R. J. Bennett. The representation and identification of spatio—temporal systems: an example of population diffusion in north-west england. Transactions of the Institute of British Geographers, 66:73 "94, 1975. R. J. Bennett. Spatial Time Series. Pion Ltd, London, 1979. J. S. Besag. Spatial interaction and the statistical analysis of lattice systems. Journal of the Royal Statistical Society Series B. 36:197—242. 1974. J .-P. Bocquet-Appel and R. R. Sokal. Spatial autocorrelation analysis of trend residuals in biological data. Systematic Zoology, 38, 4:331~341. 1989. G. E. P. Box and G. M. Jenkins. Time Series Analysis: Forecasting and Control. Holden-Day, San Francisco, 1970. Centers for Disease Control and Prevention. Outbreak of West Nile-like viral encephalitis. Morbidity Mortality Weekly Report. 48:845-849. 1999. Centers for Disease Control and Prevention. Provisional surveillance summary of the West Nile virus epidemic - United States, January-November 2002. Morbidity Mortality Weekly Report, 51:1129--1133, 2002. 110 [12] [13] [14] [15] [16} [17] [18] [19] [20] [21] [22] [23] [24] A. D. Cliff, P. Haggett. J. K. Ord. K. A. Bassett. and R. B. Davies. Elements of Spatial Structure: A Quantitative Approach. Cambridge University Press, New York, 1975. A. D. Cliff and J. K. Ord. Spatial autocorrelation. London, 1973. A. D. Cliff and J. K. Ord. Model building and the analysis of spatial pattern in human geography. Journal of the Royal Statistical Society Series B, 37:297 348- 1975. A. D. Cliff and J. K. Ord. Space-time modelling with an application to regional forecasting. Transactions of the Institute of British Geographers. 64:119-128. 1975. A. D. Cliff and J. K. Ord. Spatial processes: Models and applications. Pion Ltd.. London,1981. N. Cressie. Comment on ’An approach to statistical spatial-temporal modeling of meterological fields’ by MS. Handcock and J .R. Wallis. Journal of the American Statistical Association, 89:379-382. 1994. Division of Vector-Borne Diseases, National Center for Infectious Diseases. Cen- ters for Disease Control and Prevention. West Nile virus statistics. surveillance, and control. 2004. J. Durbin. The fitting of time series models. Review of the International Statis- tical Institute, 28:233—243, 1960. M. Eidson, N. Komar, F. Sorhage, R. Melson, T. Talbot, and F. Mostashari. Crow deaths as a sentinel surveillance system for West Nile virus in the north- eastern United States. Emerging Infectious Diseases, 7:615~620, 2001. M. Eidson, N. Komar, F. Sorhage, R. Melson, T. Talbot, and F. Mostashari. Crow deaths as a sentinel surveillance system for West Nile virus in the north- eastern united states. Emerging Infectious Diseases, 7:615—620. 2001. M. Eidson, L. Kramer, W. Stone, Y. Hagiwara. K. Schmit. and New York State West Nile Virus Surveillance Team. Dead bird surveillance as an early warning system for West Nile virus. Emerging Infectious Diseases, 726314335, 2001. M. Eidson, M. Miller, and L. Kramer. Dead crow densities and human cases of West Nile virus, New York State. Emerging Infectious Diseases, 7:662~664. 2001. B. K. Epperson- Spatial and space-time correlations in ecological models. Eco- logical Modelling, 132:63< 76. 2000. 111 [25] [26] [27] [28l [29] l30l [31] [3?] [33] [34] [35] [36] [37] B. K. Epperson. Geographical Genetics. Princeton University Press. New Jersey. 2003. S. Fitzgerald, J. Patterson, M. Kiupel. H. Simmons. S. Grimes, C. Sarver. R. Ful- ton. B. Steficek, T. Cooley, J. Massey. and J Sikarskie. Clinical and pathologic features of West Nile virus infection in native North American owls (Family Strigidae). Avian Diseases. 47(3):602 610, 2003. R. Giacomini and C. W. J. Granger. Aggregation of space—time processes. Jour- nal of Econometrics, 118:7 26, 2001. S. C. Guptill, K. G. Julian. G. L. Campbell, S. D. Price. and A. A. Marfin. Early-season avian deaths from West Nile virus as warnings of human infection. Emerging Infectious Diseases, 9:483—484, 2003. E. J. Hannan and J. Rissanen. Recursive estimation of mixed auto—regressive moving average order. Biometrika, 69:81-94. 1982. P. M. Hooper and G. J. D. Hewing. Some properties of space-time processes. Geographical Analysis, 13:203-223. 1981. H. C. Huang and N. Cressie. Spatio—temporal prediction of snow water equivalent using the kalman filter. Computational Statistics and Data Analysis. 222159 5175. 1996. C. M. Hurvich and C.-L. Tsai. Regression and time series model selection in small samples. Biometrika, 76:297~307. 1989. R. H. Jones. Fitting autoregressions. Journal of the American Statistical Asso- ciation, 70 (351)2590’592, 1975. N. Komar. West Nile viral encephalitis. Revue Scientifique et Technique, 192166 - 176, 2000. N. Komar, S. Langevin, S. Hinten, N. Nemeth, Edwards 13.. D. Hettler. B. Davis. R. Bowen, and Michel Bunning- Experimental infection of North American birds with the New York 1999 strain of West Nile viral. Emerging Infectious Diseases. 9:311—322, 2003. H. Lutkerpohl. Forecasting Aggregated Vector ARMA Processes. Springer—Verlag. Berlin, 1987. H. Lutkerpohl. Introduction to Multiple Time Series Analysis. Springer-V 5rlag. Berlin, 1993. 112 [38] [391 [40] [41] [42] [43] [44] [46] [47] [48] [49] [50] K. V. Mardia, C. R. Goodall. E. Redfern. and F. Alonso. The kriged kalman filter. Test, 7:217—285, 1998. D. W. Marquardt. An algorithm for least squares estimation of nonlinear param— eters. Journal of the Society of Industrial and Applied Mathematics. 11:431-441. 1963. R. L. Martin and J. E. Oeppen. The identification of regional forecasting models using space-time correlation functions. Transactions of the Institute of British Geographers, 66295—118. 1975. R. G. McLean, S. R. Ubico, D. E. Docherty, W. R. Hansen, L. Sielo. and T. S. McNamera. West Nile virus transmission and ecology in birds. Annals of the New York Academy of Sciences. 951:54—57, 2001. C. A. Oprian. V. S. Taneja, D. A. V088, and L. A. Aroian. General considerations and interrelationships between MA and AR models, time series in n1 dimensions. the ARMA model. Communications in Statistics - Simulation and Computation. B9 (5)2515—532. 1980. L. Petersen and J. Roehrig. West Nile virus: a reemerging global pathogen. Emerging Infectious Diseases. 7(4):611—614, 2001. P. E. Pfeifer. Spatial-Dynamic Modeling (Unpublished Ph.D. dissertation). Geor- gia Institute of Technology, Atlanta, Georgia. 1979. P. E. Pfeifer and S. J. Deutsch. A comparison of estimation procedures for the parameters of the STAR model. Communications in Statistics - Simulation and Computation, 89 (3):255—270, 1980. P. E. Pfeifer and S. J. Deutsch. Independence and sphericity test for the residuals of space-time ARMA models. Communications in Statistics - Simulation and Computation, 89 (5):533—549, 1980. P. E. Pfeifer and S. J. Deutsch. Stationarity and invertibility regions for low order STARMA models. Communications in Statistics - Simulation and Computation. 39 (5):551—562. 1980. P. E. Pfeifer and S. J. Deutsch. A three-stage iterative procedure for space-time modeling. Technometrics, 22 (1):35-—47, 1980. P. E. Pfeifer and S. J. Deutsch. Variance of the sample space-time autocorrelatimi function. Journal of the Royal Statistical Society Series B, 43 (1):28~33. 1981. C. Schwarz. Estimating the dimension of a model. The Annals of Statistics- 6 (2):461464, 1978. 113 [51] [52] [53] [54] [56] [57] [59] R. Shibata. Selection of the order of an autoregressive model by Akaike’s infor- mation criterion. Biometrika. 63:117—126. 1976. C. N. Theophilides, S. C. Ahearn, S. Grady, and M. Merlino. Identifying West Nile virus risk areas: The dynamic continuous-area space-time system. American. Journal of Epidemiology, 157:843—854, 2003. M. J. Turell, M. O’Guinn, and J. Oliver. Potential for New York mosquitoes to transmit West Nile virus. American Journal of Tropical Medicine and Hygiene. 62:413-414. 2000. M. J. Turell, M. R. Sardelis, and D. J. Dohm. Potential North American vectors of West Nile virus. Annals of the New York Academy of Sciences, 9512317324- 2001. D. A. Voss, C. A. Oprian, and L. A. Aroian. Moving average models- time series in 111-dimensions. Communications in Statistics - Simulation and Computation. 139 (5):467—489, 1980. J. Watson, R. Jones, K. Gibbs, and W. Paul. Dead crow reports and location of human West Nile virus cases, Chicago, 2002. Emerging Infectious Diseases, 10(5):938—940, 2004. S. A. Yaremych, R. J. Novak, A. Raim, P. Mankin, and R. E. Warner. Home range and habitat use by American Crows in relation to West Nile virus transmission. The Wilson Bulletin (in press), 2004a. S. A. Yaremych, R. E. Warner. P. Mankin, J. Brown, and R. J. Novak. West Nile virus and high death rate in American Crows. Emerging Infectious Diseases, 10(4):709~711, 2004b. S. A. Yaremych, R. E. Warner, M. T. Wyngaerde, A. M. Ringia, R. Lampman. and R. J. Novak. West Nile virus detection in American Crows. Emerging Infectious Diseases. 9:1319-1321, 2003. 114 lllllllllllllllllllllllllllllllllllll 3 1293 02736 1041