This work presents some common issues in the statistical analysis of time series of environmental area. The discussion and the presentation of solutions is raised by the study of a time series of the oxygen concentration variable in a water quality monitoring site in the river Vouga hydrological basin in Portugal. Issues such as trends, seasonality, temporal correlation and detection of change points are addressed

Costa, Marco

Monteiro, Magda

Marco Costa

Magda Monteiro

Crossref

A mixed-effect state space model to environmental data

Repositório Institucional da Universidade de Aveiro

A mixed-effect state space model to environmental dataMarco Costa and Magda Monteiro  Citation: AIP Conference Proceedings 1648, 110002 (2015); doi: 10.1063/1.4912409 View online: http://dx.doi.org/10.1063/1.4912409 View Table of Contents: http://scitation.aip.org/content/aip/proceeding/aipcp/1648?ver=pdfcov Published by the AIP Publishing  Articles you may be interested in Performance of mixed effects for clustered binary data models AIP Conf. Proc. 1643, 80 (2015); 10.1063/1.4907428  Reverberation models as an aid to interpret data and extract environmental information J. Acoust. Soc. Am. 136, 2297 (2014); 10.1121/1.4900309  Lee-Carter state space modeling: Application to the Malaysia mortality data AIP Conf. Proc. 1602, 1002 (2014); 10.1063/1.4882606  Investigation of the use of lower frequency acoustic data in helicopter environmental noise modeling. J. Acoust. Soc. Am. 127, 1835 (2010); 10.1121/1.3384280  Environmental Effects on Assignment and Geometry of the Triplet State of Benzaldehyde J. Chem. Phys. 57, 1809 (1972); 10.1063/1.1678489   This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to  IP:193.137.168.16 On: Wed, 11 Mar 2015 16:34:19AMixed-effect State Space Model to Environmental DataMarco Costa∗,† and Magda Monteiro∗,∗∗∗Escola Superior de Tecnologia e Gestão de Águeda, Universidade de Aveiro, 3754–909 Águeda, Portugal†Centro de Matemática e Aplicações Fundamentais, Universidade de Lisboa, 1649-003 Lisboa, Portugal∗∗Centro de Investigação e Desenvolvimento em Matemática e Aplicações, Universidade de Aveiro, 3810-193Aveiro, PortugalAbstract. This work presents some common issues in the statistical analysis of time series of environmental area. Thediscussion and the presentation of solutions is raised by the study of a time series of the oxygen concentration variable ina water quality monitoring site in the river Vouga hydrological basin in Portugal. Issues such as trends, seasonality, temporalcorrelation and detection of change points are addressed.Keywords: Time series analysis, Kalman filter, state space model, environmental data.PACS: 02.50.Ey, 05.45.TpINTRODUCTION AND DESCRIPTION OF THE DATAThe time series analysis of environmental data has been an important tool on the assessment and monitoring of naturalprocesses. However, the statistical analysis of environmental data raises a number of issues related to the characteristicsof this type of data. Oftentimes, the hydrometeorological and environmental data are monthly observations witha strong seasonal character. Some analyzes can circumvent this problem considering the annual averages, if itis appropriate ([9]). The Gaussian assumption for the data distribution is not suitable in many applications. Thisquestion is relevant due to many statistical procedures assume this assumption. When data distribution are skewed,the observations can be transformed, usually, using the logarithmic function. However, this procedure has problemsin the interpretation of the models and their parameters. A temporal correlation structure on environmental dataare very common. For instance, [6] applies state space models in order to accommodate the time dependence ofmonthly water quality variables, [1] presents inference propers of regression models which residuals follow anyautoregressive stationary process. One of the most common properties is the non-stationarity of a time series. The non-stationarity can assume many forms. For instance, the existence of a trend or heterogeneity in variance, the changeof a property in a known instant or in a unknown point (usually designed as a change-point). An environmentaltime series may have some of these properties simultaneously implying that its analysis should incorporate varioustechniques in an appropriate way. This work presents the modeling of a water quality variable time series whichpresents a set of properties. This example allows to expose statistics issues on the time series modeling and itpresents some solutions for them. This paper examines the monthly data of the Dissolved Oxygen (mg/l) in theCarvoeiro water monitoring site in the hydrological basin of the river Vouga, Portugal. Data were collected from theSNIRH (Sistema Nacional de Recursos Hídricos), the national information system for water resources of Portugal, athttp://snirh.apambiente.pt/. The modeling of water quality variables will be performed using both linearregression ([5, 11]) and state space models ([3, 6]). Linear models are usually applied for their simplicity and wellknown statistical properties. State space models associated with Kalman filter allow incorporating some componentsin the model in a dynamic point of view and the Kalman filter produces forecasts with the respective confidenceintervals.REGRESSION MODELINGThe exploratory analysis of data shows that there are a clear structural change in the trend during in the year of 2000(see Figure 2). Instead of considering two separate analysis, we intend to find a global model that accommodates thisproperty. Thus, it was followed the methodology adopted in [2]. Considering the DO variable as Yt , t = 1, ...,283, afterfitting a regression model with a quadratic trend up to a certain time t0 to all possible month/year and a constant to theProceedings of the International Conference on Numerical Analysis and Applied Mathematics 2014 (ICNAAM-2014)AIP Conf. Proc. 1648, 110002-1–110002-4; doi: 10.1063/1.4912409© 2015 AIP Publishing LLC 978-0-7354-1287-3/$30.00110002-1 This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to  IP:193.137.168.16 On: Wed, 11 Mar 2015 16:34:19		  	          	          	  		  	          	          	  FIGURE 1. Sample autocorrelation function and sample partial autocorrelation function of the residuals of the regressionanalysis.same series from the next month/year, we choose the one with a smaller residual sum of squares. This procedure ledto t0 = 139, i.e., October 2000 and the estimated model of the trend wasTt ={ −0.0014t2+0.01986t+6.8991 i f t ≤ 1399.128 i f t > 139. (1)From residuals the seasonal coefficients were estimated as suggested by [9], i.e., the easiest way to handle theseasonality is to subtract from each month the overall average of this month,sr+12k, with r = 1, ...,12 and k = 0,1, ....As will be seen this procedure is sufficient for this stage of modeling. On the one hand, from Figure 1, which representsboth the sample autocorrelation function and the partial autocorrelation function, it is clear that the residuals seriesfollows an AR(1) process, that is, there is a temporal correlation structure. On the other hand, we concluded thatsubseries variances of each month of the year are heterogeneous and they tend to have greater variability associated tothe highest seasonal coefficients (sample linear correlation coefficient equal to 0.38).A MIXED-EFFECT STATE-SPACE MODELIn order to incorporate the temporal correlation and to accommodate the heterogeneity of the variances of the subseries of each month, we adopte a mixed-effect state-space model. The model is defined asYt = (Tt +St)βt + et (2)βt = µ+φ(βt−1−µ)+ εt (3)where the trend component Tt has a structure defined in (1), St is a periodic function with the twelve averages asabove, the process {βt} is a stationary AR(1) with Gaussian errors εt with mean µ and et is the observation equationerror, assumed to be a Gaussian white noise process. The model has two main components: a regression structurewhich incorporate both trend and seasonality and an unobservable process {βt}, the state. It is assumed that stateprocess will calibrate the deterministic structure Tt + St . Another important feature of the model is that by its ownformulation it allows the variances heterogeneity previously identified. As the state process {βt} is unknown, it mustbe predicted by the Kalman filter recursions. The Kalman filter provides optimal unbiased linear one-step-aheadpredictions of the unknown states and their corresponding mean square error (MSE), which are denoted by β̂t|t−1and Pt|t−1 respectively. Furthermore, the Kalman filter provides optimal filtered predictions β̂t|t and their MSE Pt|t .If the errors are further assumed to have a conditional Normal distribution, then β̂t|t−1 is the conditional mean ofβt . Model’s parameters Θ = (µ ,φ ,σ2e ,σ2ε ) are obtained through the maximum likelihood estimation ([7]). Table 1presents the maximum likelihood estimates of Θ. In this paper we adopted a decomposition approach type to estimatethe model’s parameters, that is, estimating parameters step-by-step in each phase of the modeling procedure. However,it is possible to estimate the model (2)-(3) through the EM-algorithm method (see for instance [10]). The forecastsone step-ahead with the respective empirical confidence interval at 95% were obtained by Kalman filter equation,that is, Ŷt|t−1±1.96√MSEŶt|t−1 . Figure 2 shows the observed DO concentration, the forecast one step-ahead with therespective confidence interval at 95%. The percentage of observations outside of the respective empirical confidenceinterval is 6.36% (18/283). The residuals analysis shows that there are some outliers; the biggest residual is respected110002-2 This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to  IP:193.137.168.16 On: Wed, 11 Mar 2015 16:34:19                     			FIGURE 2. Dissolved oxygen concentration and the forecast one step-ahead with the respective confidence intervals at 95%.TABLE 1. Gaussian maximum likelihoodestimates of the state space model.µ̂ φ̂ σ̂2ε σ̂2e1.00077 0.43001 0.01790 0.24292to January 1992. The p-value of the Kolmogorov-Smirnov test to the normality is equal to 2.8% leading to rejectionof null hypothesis. However, the Gaussian distribution is not rejected (K-S p-value=20.0%) when the largest outliersare not considered.CHANGE POINT DETECTIONThe modeling procedure exposed before shows a change in the trend component of the DO concentration. This factis clear from the exploratory analysis. However, other structural changes are not so evident and they need a statisticalanalysis in order to evaluate their existence. The model (2)-(3) established a calibration procedure considering Tt +Sta intrinsic component of the variable which, in a dynamic way, is calibrate by the state process {βt}. Thus, it isreasonable to investigate changes in this structure which may be relevant in the water monitoring process. For thispurpose, the filtered predictions of the state process, β̂t|t can be analyzed relatively to the existence of a change (ormore) in their mean as the best prediction of an AR(1) process. Since the sample size is considerable, the results on theGaussian distribution discussed before was devalued. We used test statistics that are called maximum type and referredin [4]. The null hypothesis claims that X1, X2, ..., Xn are distributed according to N(µ,σ2). The alternative claims thatthere exists a time point k ∈ {1,2, ...,n− 1} such that X1, X2, ..., Xk are distributed according to N(µ1,σ2) and Xk+1,..., Xn are distributed according to N(µ2,σ2) with µ1 ̸= µ2. Supposing σ2 is unknown, the problem described can betreated based on two-sample t tests. The test statistic T (n) is the maximum of the absolute values of two sample t-teststatisticsT (n) = max1≤k<n|Tk|= max1≤k<n√(n− k)kn∣∣∣β k|k−β ∗k|k∣∣∣ 1skwhere β̂ k|k = 1k ∑ki=1 β̂i|i, β̂∗k|k = 1n−k ∑ni=k+1 β̂i|i and sk =√√√√ 1n−2[∑ki=1(β̂i|i− β̂ k|k)2+∑ni=k+1(β̂i|i− β̂∗k|k)2].The test statistic T (n) were computed considering the sequences of βt|t . Based on [9], the interpolated value ton = 283 is 3.207. However, this critical value assumes that the observations are uncorrelated. So, it is necessary tocorrect this critical value according to [8]. Thus, in order to obtain 5% critical values, the value 3.207 is multiplied by[(1+ φ̂)(1− φ̂)−1]1/2, where φ̂ is obtained in Table 1. Thus, 5% critical value is 5.079. Table 2 presents the resultsfrom a binary segmentation procedure (proposed initially by [12]) in order to detect multiple change point. Thisprocedure ends when the observed value of the test statistics is less than the respective critical value (the critical valueis obtained by interpolation according to the subsample size). This procedure identified two change points, April 1990and December 1991. Figure 3 represents the change points and the three mean levels of the calibration factors process.For chronological order, the mean levels of the process {βt} are 1.190, 0.835 and 1.005. This means that in the first110002-3 This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to  IP:193.137.168.16 On: Wed, 11 Mar 2015 16:34:19TABLE 2. Results from the binary segmentation procedure.ni T (ni)obs month/yearcritical valuenot correctedcritical valuecorrected283 5.1893 Apr-90 3.2066 5.0790270 5.7600 Dec-91 3.2040 5.0749!#%!#$%$%%%%%%%%% %!	%!%"%"%#%$%$%%%%	 !!""##$%	%FIGURE 3. Representation of the change points and the three mean levels of the {β̂t}.period of time the DO variable was, on average, 19% higher than those predicted by the regression component. In thesecond period of time the observations were, on average, 16.5% lower and in last period of time the observation were,on average, very close to the regression model predictions. From the environmental point of view, the analysis of thecalibration factor in each time can be relevant in order to identify unanticipated changes.ACKNOWLEDGMENTSMarco Costa was partially supported by Fundação para a Ciência e a Tecnologia, PEst OE/ MAT/ UI0209/ 2011.Magda Monteiro was partially supported by Portuguese funds through the CIDMA - Center for Research and Develop-ment in Mathematics and Applications, and the Portuguese Foundation for Science and Technology ("FCT–Fundaçãopara a Ciência e a Tecnologia"), within project PEst-OE/MAT/UI4106/2014.REFERENCES1. T. Alpuim, and A. El-Shaarawi, On the efficiency of regression analysis with AR(p) errors, Journal of Applied Statistics 35,717–737 (2008).2. T. Alpuim, and A. El-Shaarawi, Modeling monthly temperature data in Lisbon and Prague, Environmetrics 20, 835–852 (2009).3. M. Costa, and A.M. Gonçalves, Clustering and forecasting of dissolved oxygen concentration on a river basin, Stoch. Env. Res.Risk A. 25, 151–163 (2011).4. M. Costa, and A.M. Gonçalves, Application of Change-Point Detection to a Structural Component of Water Quality Variables,In Numerical Analysis and Applied Mathematics ICNAAM 2011, edited by S. Theodore et al., AIP Conference Proceedings1389, American Institute of Physics, New York, 2011, pp. 1565–1568.5. A.M.. Gonçalves, and T. Alpuim, Water Quality Monitoring using Cluster Analysis and Linear Models, Environmetrics 22,933–945 (2011).6. A.M. Gonçalves, and M. Costa, Predicting seasonal and hydro-meteorological impact in environmental variables modelling viaKalman filtering, Stoch. Env. Res. Risk A. 27, 1021–1038 (2013).7. A. Harvey, Forecasting, Structural Time Series Models and the Kalman Filter, Cambridge University Press, Cambridge.8. D. Jarušková, Change-point detection meteorological measurement, Mon Weather Rev 124, 1535–1543 (1996).9. D. Jarusˇková, Some problems with application of change-point detection methods to environmental data, Environmetrics 8,469–483 (1997).10. P. Kokic, S. Crimp, and M. Howden, Forecasting climate variables using a mixed-effect state-space model, Environmetrics 22,409–419 (2011).11. J.A. Renwick, A.B. Mullan, and A. Porteous, Statistical downscaling of New Zealand climate,Weather and Climate 29, 24–44(2009).12. L.J. Vostrikova, Detecting "disorder" in multidimensional random processes, Soviet Mathematics Doklady 24, 55–59 (2009).110002-4 This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to  IP:193.137.168.16 On: Wed, 11 Mar 2015 16:34:19

A Mixed-effect State Space Model to Environmental Data

http://ria.ua.pt/bitstream/10773/13619/1/A%20mixed-effect%20state%20space%20model%20to%20environmental%20data%20Costa%20Monteiro%20ICNAAM2015.pdf

A Mixed-effect State Space Model to Environmental Data

Abstract

Similar works

Full text

Available Versions

Crossref

Repositório Institucional da Universidade de Aveiro